Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

implicit BCs and Laplacian operator #262

Open
wants to merge 20 commits into
base: stack/implicitOperators
Choose a base branch
from
Open

implicit BCs and Laplacian operator #262

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
a529edb
Fix operator & return a scalar instead of the vector
HenningScheufler Feb 15, 2025
b6d159f
geometricScheme compute delta and nonOrthoDelta
HenningScheufler Feb 15, 2025
602d5e6
ENH: added facenormalgradient
HenningScheufler Feb 15, 2025
45d7786
WIP added LaplacianOperator
HenningScheufler Feb 15, 2025
e698a60
FIX: 1DMesh face ordering
HenningScheufler Feb 15, 2025
a437e47
add next functionality to be able to consume the tokens
HenningScheufler Feb 15, 2025
7ad1ed2
operators consume tokens similar to a stream
HenningScheufler Feb 15, 2025
1bd7e16
fix: faceNormalGradient
HenningScheufler Feb 15, 2025
7f2259f
implementLaplacian
HenningScheufler Feb 15, 2025
8697b27
laplacian uses fixedValue BCs
HenningScheufler Feb 16, 2025
8b43753
added deltaCoeffs to surfaceNormalGradient
HenningScheufler Feb 16, 2025
0c134bd
FIX: valueFraction return the wrong type
HenningScheufler Feb 16, 2025
dda140d
BC are handled as mixed BC as default to avoid the dispatch per BC
HenningScheufler Feb 16, 2025
fd19c94
Merge branch 'stack/implicitOperators' into feat/implicitBCs
HenningScheufler Feb 16, 2025
8b1a5f5
added implicit Laplacian operator
HenningScheufler Feb 16, 2025
d05ac35
added default constructors and outOfRange check
HenningScheufler Feb 16, 2025
3928a2b
FIX: laplacian operator check of div
HenningScheufler Feb 16, 2025
e7e1154
Fix: scaling of operators sign inside the dsl was not accounted for
HenningScheufler Feb 22, 2025
1abf469
FIX: solverdict always stores doubles
HenningScheufler Feb 22, 2025
2f0eca4
FIX: solutionDict in backwardEuler no longer hardcoded
HenningScheufler Feb 22, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions include/NeoFOAM/dsl/expression.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -195,6 +195,15 @@ template<typename leftOperator, typename rightOperator>
return lhs;
}

template<typename leftOperator, typename rightOperator>
[[nodiscard]] inline Expression operator-(leftOperator lhs, rightOperator rhs)
{
Expression expr(lhs.exec());
expr.addOperator(lhs);
expr.addOperator(Coeff(-1) * rhs);
return expr;
}

[[nodiscard]] inline Expression operator-(const SpatialOperator& lhs, const SpatialOperator& rhs)
{
Expression expr(lhs.exec());
Expand Down
45 changes: 45 additions & 0 deletions include/NeoFOAM/dsl/operator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,4 +16,49 @@ class Operator
};
};


/* @class OperatorMixin
* @brief A mixin class to simplify implementations of concrete operators
* in NeoFOAMs dsl
*
* @ingroup dsl
*/
template<typename FieldType>
class OperatorMixin
{

public:

OperatorMixin(const Executor exec, const Coeff& coeffs, FieldType& field, Operator::Type type)
: exec_(exec), coeffs_(coeffs), field_(field), type_(type) {};


Operator::Type getType() const { return type_; }

virtual ~OperatorMixin() = default;

virtual const Executor& exec() const final { return exec_; }

Coeff& getCoefficient() { return coeffs_; }

const Coeff& getCoefficient() const { return coeffs_; }

FieldType& getField() { return field_; }

const FieldType& getField() const { return field_; }

/* @brief Given an input this function reads required coeffs */
void build([[maybe_unused]] const Input& input) {}

protected:

const Executor exec_; //!< Executor associated with the field. (CPU, GPU, openMP, etc.)

Coeff coeffs_;

FieldType& field_;

Operator::Type type_;
};

} // namespace NeoFOAM::dsl
44 changes: 2 additions & 42 deletions include/NeoFOAM/dsl/spatialOperator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,8 @@ class SpatialOperator

SpatialOperator(SpatialOperator&& eqnOperator);

SpatialOperator& operator=(const SpatialOperator& eqnOperator);

void explicitOperation(Field<scalar>& source);

void implicitOperation(la::LinearSystem<scalar, localIdx>& ls);
Expand Down Expand Up @@ -213,46 +215,4 @@ SpatialOperator operator*([[maybe_unused]] CoeffFunction coeffFunc, const Spatia
return result;
}

/* @class OperatorMixin
* @brief A mixin class to simplify implementations of concrete operators
* in NeoFOAMs dsl
*
* @ingroup dsl
*/
template<typename FieldType>
class OperatorMixin
{

public:

OperatorMixin(const Executor exec, FieldType& field, Operator::Type type)
: exec_(exec), coeffs_(), field_(field), type_(type) {};

Operator::Type getType() const { return type_; }

virtual ~OperatorMixin() = default;

virtual const Executor& exec() const final { return exec_; }

Coeff& getCoefficient() { return coeffs_; }

const Coeff& getCoefficient() const { return coeffs_; }

FieldType& getField() { return field_; }

const FieldType& getField() const { return field_; }

/* @brief Given an input this function reads required coeffs */
void build([[maybe_unused]] const Input& input) {}

protected:

const Executor exec_; //!< Executor associated with the field. (CPU, GPU, openMP, etc.)

Coeff coeffs_;

FieldType& field_;

Operator::Type type_;
};
} // namespace NeoFOAM::dsl
2 changes: 1 addition & 1 deletion include/NeoFOAM/finiteVolume/cellCentred/operators/ddtOperator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ class DdtOperator : public dsl::OperatorMixin<VolumeField<scalar>>
public:

DdtOperator(dsl::Operator::Type termType, VolumeField<scalar>& field)
: dsl::OperatorMixin<VolumeField<scalar>>(field.exec(), field, termType),
: dsl::OperatorMixin<VolumeField<scalar>>(field.exec(), dsl::Coeff(1.0), field, termType),
sparsityPattern_(SparsityPattern::readOrCreate(field.mesh())) {

};
Expand Down
62 changes: 42 additions & 20 deletions include/NeoFOAM/finiteVolume/cellCentred/operators/divOperator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,29 +48,44 @@ class DivOperatorFactory :
// does not work and we cant template the entire class because the static create function
// cannot access keyExistsOrError and table anymore.
virtual void
div(VolumeField<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi
) = 0;
div(VolumeField<scalar>& divPhi,
const SurfaceField<scalar>& faceFlux,
VolumeField<scalar>& phi,
const dsl::Coeff operatorScaling) = 0;

virtual void
div(la::LinearSystem<scalar, localIdx>& ls,
const SurfaceField<scalar>& faceFlux,
VolumeField<scalar>& phi) = 0;
VolumeField<scalar>& phi,
const dsl::Coeff operatorScaling) = 0;

virtual void
div(Field<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi) = 0;
div(Field<scalar>& divPhi,
const SurfaceField<scalar>& faceFlux,
VolumeField<scalar>& phi,
const dsl::Coeff operatorScaling) = 0;

virtual void
div(VolumeField<Vector>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi
) = 0;
div(VolumeField<Vector>& divPhi,
const SurfaceField<scalar>& faceFlux,
VolumeField<Vector>& phi,
const dsl::Coeff operatorScaling) = 0;

virtual void
div(Field<Vector>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi) = 0;
div(Field<Vector>& divPhi,
const SurfaceField<scalar>& faceFlux,
VolumeField<Vector>& phi,
const dsl::Coeff operatorScaling) = 0;

virtual VolumeField<scalar>
div(const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi) = 0;
div(const SurfaceField<scalar>& faceFlux,
VolumeField<scalar>& phi,
const dsl::Coeff operatorScaling) = 0;

virtual VolumeField<Vector>
div(const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi) = 0;
div(const SurfaceField<scalar>& faceFlux,
VolumeField<Vector>& phi,
const dsl::Coeff operatorScaling) = 0;

// Pure virtual function for cloning
virtual std::unique_ptr<DivOperatorFactory> clone() const = 0;
Expand All @@ -90,7 +105,9 @@ class DivOperator : public dsl::OperatorMixin<VolumeField<ValueType>>

// copy constructor
DivOperator(const DivOperator& divOp)
: dsl::OperatorMixin<VolumeField<ValueType>>(divOp.exec_, divOp.field_, divOp.type_),
: dsl::OperatorMixin<VolumeField<ValueType>>(
divOp.exec_, divOp.coeffs_, divOp.field_, divOp.type_
),
faceFlux_(divOp.faceFlux_),
divOperatorStrategy_(
divOp.divOperatorStrategy_ ? divOp.divOperatorStrategy_->clone() : nullptr
Expand All @@ -102,7 +119,7 @@ class DivOperator : public dsl::OperatorMixin<VolumeField<ValueType>>
VolumeField<ValueType>& phi,
Input input
)
: dsl::OperatorMixin<VolumeField<ValueType>>(phi.exec(), phi, termType),
: dsl::OperatorMixin<VolumeField<ValueType>>(phi.exec(), dsl::Coeff(1.0), phi, termType),
faceFlux_(faceFlux),
divOperatorStrategy_(DivOperatorFactory::create(phi.exec(), phi.mesh(), input)) {};

Expand All @@ -112,16 +129,16 @@ class DivOperator : public dsl::OperatorMixin<VolumeField<ValueType>>
VolumeField<ValueType>& phi,
std::unique_ptr<DivOperatorFactory> divOperatorStrategy
)
: dsl::OperatorMixin<VolumeField<scalar>>(phi.exec(), phi, termType), faceFlux_(faceFlux),
divOperatorStrategy_(std::move(divOperatorStrategy)) {};
: dsl::OperatorMixin<VolumeField<scalar>>(phi.exec(), dsl::Coeff(1.0), phi, termType),
faceFlux_(faceFlux), divOperatorStrategy_(std::move(divOperatorStrategy)) {};

DivOperator(
dsl::Operator::Type termType,
const SurfaceField<scalar>& faceFlux,
VolumeField<ValueType>& phi
)
: dsl::OperatorMixin<VolumeField<scalar>>(phi.exec(), phi, termType), faceFlux_(faceFlux),
divOperatorStrategy_(nullptr) {};
: dsl::OperatorMixin<VolumeField<scalar>>(phi.exec(), dsl::Coeff(1.0), phi, termType),
faceFlux_(faceFlux), divOperatorStrategy_(nullptr) {};


void explicitOperation(Field<scalar>& source)
Expand All @@ -131,13 +148,15 @@ class DivOperator : public dsl::OperatorMixin<VolumeField<ValueType>>
NF_ERROR_EXIT("DivOperatorStrategy not initialized");
}
NeoFOAM::Field<NeoFOAM::scalar> tmpsource(source.exec(), source.size(), 0.0);
divOperatorStrategy_->div(tmpsource, faceFlux_, this->getField());
const auto operatorScaling = this->getCoefficient();
divOperatorStrategy_->div(tmpsource, faceFlux_, this->getField(), operatorScaling);
source += tmpsource;
}

void div(Field<ValueType>& divPhi)
{
divOperatorStrategy_->div(divPhi, faceFlux_, this->getField());
const auto operatorScaling = this->getCoefficient();
divOperatorStrategy_->div(divPhi, faceFlux_, this->getField(), operatorScaling);
}

la::LinearSystem<scalar, localIdx> createEmptyLinearSystem() const
Expand All @@ -155,18 +174,21 @@ class DivOperator : public dsl::OperatorMixin<VolumeField<ValueType>>
{
NF_ERROR_EXIT("DivOperatorStrategy not initialized");
}
divOperatorStrategy_->div(ls, faceFlux_, this->getField());
const auto operatorScaling = this->getCoefficient();
divOperatorStrategy_->div(ls, faceFlux_, this->getField(), operatorScaling);
}


void div(la::LinearSystem<scalar, localIdx>& ls)
{
divOperatorStrategy_->div(ls, faceFlux_, this->getField());
const auto operatorScaling = this->getCoefficient();
divOperatorStrategy_->div(ls, faceFlux_, this->getField(), operatorScaling);
};

void div(VolumeField<scalar>& divPhi)
{
divOperatorStrategy_->div(divPhi, faceFlux_, this->getField());
const auto operatorScaling = this->getCoefficient();
divOperatorStrategy_->div(divPhi, faceFlux_, this->getField(), operatorScaling);
}


Expand Down
35 changes: 24 additions & 11 deletions include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenDiv.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,31 +29,44 @@ class GaussGreenDiv : public DivOperatorFactory::Register<GaussGreenDiv>


virtual void
div(VolumeField<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess what we can do to reduce the amount of variants here is to move the dispatching to dsl/operators.hpp and support only one concrete implementation per ValueType ie. div(span<scalar> out, ), div(span<Vector> out, ) and one for Tensor ...

) override;
div(VolumeField<scalar>& divPhi,
const SurfaceField<scalar>& faceFlux,
VolumeField<scalar>& phi,
const dsl::Coeff operatorScaling) override;

virtual void
div(la::LinearSystem<scalar, localIdx>& ls,
const SurfaceField<scalar>& faceFlux,
VolumeField<scalar>& phi) override;
VolumeField<scalar>& phi,
const dsl::Coeff operatorScaling) override;

virtual void
div(Field<scalar>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi
) override;
div(Field<scalar>& divPhi,
const SurfaceField<scalar>& faceFlux,
VolumeField<scalar>& phi,
const dsl::Coeff operatorScaling) override;

virtual VolumeField<scalar>
div(const SurfaceField<scalar>& faceFlux, VolumeField<scalar>& phi) override;
div(const SurfaceField<scalar>& faceFlux,
VolumeField<scalar>& phi,
const dsl::Coeff operatorScaling) override;

virtual void
div(VolumeField<Vector>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi
) override;
div(VolumeField<Vector>& divPhi,
const SurfaceField<scalar>& faceFlux,
VolumeField<Vector>& phi,
const dsl::Coeff operatorScaling) override;

virtual void
div(Field<Vector>& divPhi, const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi
) override;
div(Field<Vector>& divPhi,
const SurfaceField<scalar>& faceFlux,
VolumeField<Vector>& phi,
const dsl::Coeff operatorScaling) override;

virtual VolumeField<Vector>
div(const SurfaceField<scalar>& faceFlux, VolumeField<Vector>& phi) override;
div(const SurfaceField<scalar>& faceFlux,
VolumeField<Vector>& phi,
const dsl::Coeff operatorScaling) override;

std::unique_ptr<DivOperatorFactory> clone() const override;

Expand Down
13 changes: 10 additions & 3 deletions include/NeoFOAM/finiteVolume/cellCentred/operators/gaussGreenLaplacian.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,17 +29,24 @@ class GaussGreenLaplacian : public LaplacianOperatorFactory::Register<GaussGreen
la::LinearSystem<scalar, localIdx> createEmptyLinearSystem() const override;

virtual void laplacian(
VolumeField<scalar>& lapPhi, const SurfaceField<scalar>& gamma, VolumeField<scalar>& phi
VolumeField<scalar>& lapPhi,
const SurfaceField<scalar>& gamma,
VolumeField<scalar>& phi,
const dsl::Coeff operatorScaling
) override;

virtual void laplacian(
Field<scalar>& lapPhi, const SurfaceField<scalar>& gamma, VolumeField<scalar>& phi
Field<scalar>& lapPhi,
const SurfaceField<scalar>& gamma,
VolumeField<scalar>& phi,
const dsl::Coeff operatorScaling
) override;

virtual void laplacian(
la::LinearSystem<scalar, localIdx>& ls,
const SurfaceField<scalar>& gamma,
VolumeField<scalar>& phi
VolumeField<scalar>& phi,
const dsl::Coeff operatorScaling
) override;

std::unique_ptr<LaplacianOperatorFactory> clone() const override;
Expand Down
Loading
Loading