diff --git a/kernel/include/Coefficients/AxiCylindricalCoefficient.hpp b/kernel/include/Coefficients/AxiCylindricalCoefficient.hpp new file mode 100644 index 00000000..9c19afc2 --- /dev/null +++ b/kernel/include/Coefficients/AxiCylindricalCoefficient.hpp @@ -0,0 +1,43 @@ +/** + * @file AxiCylindricalCoefficient.hpp + * @author Clément Introïni (clement.introini@cea.fr) + * @brief Built a MFEM Coefficient for AxiCylindricalCoefficient + * @version 0.1 + * @date 2026-06-16 + * + * @copyright CEA (C) 2026 + * + * This file is part of SLOTH. + * + * SLOTH is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * SLOTH is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + * + */ +#pragma once +#include +#include +#include +#include + +#include "Coefficients/Coefficients.hpp" +#include "mfem.hpp" // NOLINT [no include the directory when naming mfem include file] + +class AxiCylindricalCoefficient : public mfem::Coefficient { + private: + mutable mfem::Vector transip; + + public: + AxiCylindricalCoefficient(); + + mfem::real_t Eval(mfem::ElementTransformation& T, const mfem::IntegrationPoint& ip) override; +}; \ No newline at end of file diff --git a/kernel/include/Coefficients/Coefficients.hpp b/kernel/include/Coefficients/Coefficients.hpp index f311ace2..eb1a70a2 100644 --- a/kernel/include/Coefficients/Coefficients.hpp +++ b/kernel/include/Coefficients/Coefficients.hpp @@ -62,8 +62,9 @@ class Coefficients { void add(Coefficient coef); std::vector getCoefficients() const; - size_t size() noexcept; - Coefficient operator[](size_t i); + size_t size() const noexcept; + Coefficient& operator[](size_t i); + const Coefficient& operator[](size_t i) const; std::vector get_types() const; diff --git a/kernel/include/Coefficients/ListCoefficients.hpp b/kernel/include/Coefficients/ListCoefficients.hpp index 1e67edd3..774d06fb 100644 --- a/kernel/include/Coefficients/ListCoefficients.hpp +++ b/kernel/include/Coefficients/ListCoefficients.hpp @@ -25,6 +25,7 @@ */ #pragma once +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Coefficients/CommonCoefficients.hpp" #include "Coefficients/ProductFunctionCoefficients.hpp" #include "Coefficients/SumFunctionCoefficients.hpp" diff --git a/kernel/include/Coefficients/SlothBaseCoefficient.hpp b/kernel/include/Coefficients/SlothBaseCoefficient.hpp index 58c74029..2d06a7e8 100644 --- a/kernel/include/Coefficients/SlothBaseCoefficient.hpp +++ b/kernel/include/Coefficients/SlothBaseCoefficient.hpp @@ -117,5 +117,5 @@ class SlothBaseCoefficient { void set_time(double time); void set_bdr_index_coef(std::vector ids); - std::vector get_bdr_index_coef(); + std::vector get_bdr_index_coef() const; }; diff --git a/kernel/include/Integrators/AllenCahnNLFormIntegrator.hpp b/kernel/include/Integrators/AllenCahnNLFormIntegrator.hpp index 087b5894..d97dc869 100644 --- a/kernel/include/Integrators/AllenCahnNLFormIntegrator.hpp +++ b/kernel/include/Integrators/AllenCahnNLFormIntegrator.hpp @@ -72,7 +72,7 @@ class AllenCahnNLFormIntegrator : public SlothNLFormIntegrator { public: void init() override; - AllenCahnNLFormIntegrator(const std::vector& u_old, + AllenCahnNLFormIntegrator(Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/AllenCahnNLFormIntegrator.tpp b/kernel/include/Integrators/AllenCahnNLFormIntegrator.tpp index 546f219a..ee9286d3 100644 --- a/kernel/include/Integrators/AllenCahnNLFormIntegrator.tpp +++ b/kernel/include/Integrators/AllenCahnNLFormIntegrator.tpp @@ -33,6 +33,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Coefficients/SlothBaseCoefficient.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" #include "MAToolsProfiling/MATimersAPI.hxx" @@ -60,10 +61,10 @@ */ template AllenCahnNLFormIntegrator::AllenCahnNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : SlothNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : SlothNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->integrator_name_ = "AllenCahn"; this->check_variables_consistency(); @@ -167,15 +168,21 @@ void AllenCahnNLFormIntegrator::AssembleElementVector( std::span(vaux_gf_at_ip)); const double lamb = this->compute_coefficient(lambda[blk], std::span(u_values), std::span(vaux_gf_at_ip)); - gradU *= coef_mobi * lamb * ip.weight * Tr.Weight(); - ; + + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + + gradU *= coef_mobi * lamb * weight_coef; + gradPsi.AddMult(gradU, *elvect[blk]); // Given u, compute (w'(u), psi), psi is shape function double ww = coef_mobi * this->compute_gradient_energy_coefficient( double_well_energy[blk], blk, std::span(u_values), std::span(vaux_gf_at_ip)); - ww *= ip.weight * Tr.Weight(); + ww *= weight_coef; add(*elvect[blk], ww, Psi, *elvect[blk]); } } @@ -246,7 +253,12 @@ void AllenCahnNLFormIntegrator::AssembleElementGrad( const double lamb = this->compute_coefficient(lambda[blk], std::span(u_values), std::span(vaux_gf_at_ip)); - const double mob_lamb = coef_mobi * lamb * ip.weight * Tr.Weight(); + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + + const double mob_lamb = coef_mobi * lamb * weight_coef; AddMult_a_AAt(mob_lamb, gradPsi, *elmats(blk, blk)); // Compute w'(u)*(du,psi), psi is shape function ( // w''(u)) @@ -254,7 +266,7 @@ void AllenCahnNLFormIntegrator::AssembleElementGrad( coef_mobi * this->compute_hessian_coefficient(double_well_energy[blk], blk, blk, std::span(u_values), std::span(vaux_gf_at_ip)); - fun_val *= ip.weight * Tr.Weight(); + fun_val *= weight_coef; AddMult_a_VVt(fun_val, Psi, *elmats(blk, blk)); // w'(u)*(du, psi) } @@ -286,7 +298,12 @@ void AllenCahnNLFormIntegrator::AssembleElementGrad( coef_mobi * this->compute_hessian_coefficient(double_well_energy[jblk], blk, jblk, std::span(u_values), std::span(vaux_gf_at_ip)); - fun_val *= ip.weight * Tr.Weight(); + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + + fun_val *= weight_coef; AddMult_a_VVt(fun_val, Psi, *elmats(blk, jblk)); } } @@ -316,7 +333,13 @@ void AllenCahnNLFormIntegrator::AssembleElementGrad( coef_mobi * this->compute_hessian_coefficient(double_well_energy[jblk], blk, jblk, std::span(u_values), std::span(vaux_gf_at_ip)); - fun_val *= ip.weight * Tr.Weight(); + + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + + fun_val *= weight_coef; AddMult_a_VVt(fun_val, Psi, *elmats(blk, jblk)); } } diff --git a/kernel/include/Integrators/BlockAllenCahnNLFormIntegrator.hpp b/kernel/include/Integrators/BlockAllenCahnNLFormIntegrator.hpp index 36f4e480..5015ed7a 100644 --- a/kernel/include/Integrators/BlockAllenCahnNLFormIntegrator.hpp +++ b/kernel/include/Integrators/BlockAllenCahnNLFormIntegrator.hpp @@ -65,7 +65,7 @@ class BlockAllenCahnNLFormIntegrator : public SlothNLFormIntegrator { public: void init() override; - BlockAllenCahnNLFormIntegrator(const std::vector& u_old, + BlockAllenCahnNLFormIntegrator(Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/BlockAllenCahnNLFormIntegrator.tpp b/kernel/include/Integrators/BlockAllenCahnNLFormIntegrator.tpp index 9d7b913c..d7464ca8 100644 --- a/kernel/include/Integrators/BlockAllenCahnNLFormIntegrator.tpp +++ b/kernel/include/Integrators/BlockAllenCahnNLFormIntegrator.tpp @@ -32,6 +32,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Coefficients/SlothBaseCoefficient.hpp" #include "Integrators/SlothGridFunction.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" @@ -59,10 +60,10 @@ */ template BlockAllenCahnNLFormIntegrator::BlockAllenCahnNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : SlothNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : SlothNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->check_variables_consistency(); } @@ -163,18 +164,23 @@ void BlockAllenCahnNLFormIntegrator::AssembleElementVector( u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip); } - const double xx = ip.weight * Tr.Weight(); + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + el[blk]->CalcPhysDShape(Tr, gradPsi); gradPsi.MultTranspose(*elfun[blk], gradU); - gradU *= xx * this->compute_coefficient(lambda[blk], std::span(u_values), - std::span(vaux_gf_at_ip)); + gradU *= + weight_coef * this->compute_coefficient(lambda[blk], std::span(u_values), + std::span(vaux_gf_at_ip)); gradPsi.AddMult(gradU, *elvect[blk]); const double ww = - xx * (eta + this->compute_gradient_energy_coefficient( - double_well_energy[blk], blk, std::span(u_values), - std::span(vaux_gf_at_ip))); + weight_coef * (eta + this->compute_gradient_energy_coefficient( + double_well_energy[blk], blk, std::span(u_values), + std::span(vaux_gf_at_ip))); add(*elvect[blk], ww, Psi, *elvect[blk]); } @@ -212,10 +218,15 @@ void BlockAllenCahnNLFormIntegrator::AssembleElementVector( u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip); } + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + const double coef_mob = this->compute_coefficient(mobility[blk], std::span(u_values), std::span(vaux_gf_at_ip)) * - ip.weight * Tr.Weight(); + weight_coef; const double ww = -coef_mob * eta; @@ -281,18 +292,22 @@ void BlockAllenCahnNLFormIntegrator::AssembleElementGrad( u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip); } - const double xx = ip.weight * Tr.Weight(); + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + el[blk]->CalcPhysDShape(Tr, gradPsi); const double coef_lambda = this->compute_coefficient( lambda[blk], std::span(u_values), std::span(vaux_gf_at_ip)); - AddMult_a_AAt(xx * coef_lambda, gradPsi, *elmats(blk, blk)); + AddMult_a_AAt(weight_coef * coef_lambda, gradPsi, *elmats(blk, blk)); double fun_val = - xx * this->compute_hessian_coefficient(double_well_energy[blk], blk, blk, - std::span(u_values), - std::span(vaux_gf_at_ip)); + weight_coef * this->compute_hessian_coefficient(double_well_energy[blk], blk, blk, + std::span(u_values), + std::span(vaux_gf_at_ip)); AddMult_a_VVt(fun_val, Psi, *elmats(blk, blk)); } } @@ -318,8 +333,11 @@ void BlockAllenCahnNLFormIntegrator::AssembleElementGrad( el[blk]->CalcShape(ip, Psi); Tr.SetIntPoint(&ip); - double ww = ip.weight * Tr.Weight(); - AddMult_a_VVt(ww, Psi, *elmats(blk, off_blk)); + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + AddMult_a_VVt(weight_coef, Psi, *elmats(blk, off_blk)); } } // Block 1 0 dR(eta)dphi=d(M eta )/dphi @@ -368,10 +386,16 @@ void BlockAllenCahnNLFormIntegrator::AssembleElementGrad( u_values[off_blk] = (*elfun[off_blk]) * Psi; u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip); } + + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + const double coef_mob = -this->compute_coefficient(mobility[blk], std::span(u_values), std::span(vaux_gf_at_ip)) * - ip.weight * Tr.Weight(); + weight_coef; AddMult_a_VVt(coef_mob, Psi, *elmats(blk, blk)); } diff --git a/kernel/include/Integrators/CahnHilliardNLFormIntegrator.hpp b/kernel/include/Integrators/CahnHilliardNLFormIntegrator.hpp index b069a0b4..1265d7a3 100644 --- a/kernel/include/Integrators/CahnHilliardNLFormIntegrator.hpp +++ b/kernel/include/Integrators/CahnHilliardNLFormIntegrator.hpp @@ -66,7 +66,7 @@ class CahnHilliardNLFormIntegrator : public SlothNLFormIntegrator { public: void init() override; - CahnHilliardNLFormIntegrator(const std::vector& u_old, + CahnHilliardNLFormIntegrator(Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/CahnHilliardNLFormIntegrator.tpp b/kernel/include/Integrators/CahnHilliardNLFormIntegrator.tpp index e19b0b57..85bab459 100644 --- a/kernel/include/Integrators/CahnHilliardNLFormIntegrator.tpp +++ b/kernel/include/Integrators/CahnHilliardNLFormIntegrator.tpp @@ -31,6 +31,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Coefficients/SlothBaseCoefficient.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" #include "MAToolsProfiling/MATimersAPI.hxx" @@ -57,10 +58,10 @@ */ template CahnHilliardNLFormIntegrator::CahnHilliardNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : SlothNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : SlothNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->integrator_name_ = "CahnHilliard"; this->check_variables_consistency(); } @@ -194,19 +195,23 @@ void CahnHilliardNLFormIntegrator::AssembleElementVector( u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip); } - const double xx = ip.weight * Tr.Weight(); + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } el[blk]->CalcPhysDShape(Tr, gradPsi); gradPsi.MultTranspose(*elfun[blk], gradU); - gradU *= -xx * this->compute_coefficient(lambda[blk], std::span(u_values), - std::span(vaux_gf_at_ip)); + gradU *= + -weight_coef * this->compute_coefficient(lambda[blk], std::span(u_values), + std::span(vaux_gf_at_ip)); gradPsi.AddMult(gradU, *elvect[blk]); // Given u, compute (w'(u), psi), psi is shape function const double ww = - xx * (mu - this->compute_gradient_energy_coefficient( - double_well_energy[blk], blk, std::span(u_values), - std::span(vaux_gf_at_ip))); + weight_coef * (mu - this->compute_gradient_energy_coefficient( + double_well_energy[blk], blk, std::span(u_values), + std::span(vaux_gf_at_ip))); add(*elvect[blk], ww, Psi, *elvect[blk]); } @@ -243,11 +248,16 @@ void CahnHilliardNLFormIntegrator::AssembleElementVector( u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip); } + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + el[blk]->CalcPhysDShape(Tr, gradPsi); gradPsi.MultTranspose(*elfun[blk], gradU); gradU *= this->compute_coefficient(mobility[blk], std::span(u_values), std::span(vaux_gf_at_ip)) * - ip.weight * Tr.Weight(); + weight_coef; gradPsi.AddMult(gradU, *elvect[blk]); } } @@ -312,19 +322,22 @@ void CahnHilliardNLFormIntegrator::AssembleElementGrad( u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip); } - const double xx = -ip.weight * Tr.Weight(); + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } el[blk]->CalcPhysDShape(Tr, gradPsi); const double coef_lambda = this->compute_coefficient( lambda[blk], std::span(u_values), std::span(vaux_gf_at_ip)); - AddMult_a_AAt(xx * coef_lambda, gradPsi, *elmats(blk, blk)); + AddMult_a_AAt(-weight_coef * coef_lambda, gradPsi, *elmats(blk, blk)); double fun_val = - xx * this->compute_hessian_coefficient(double_well_energy[blk], blk, blk, - std::span(u_values), - std::span(vaux_gf_at_ip)); + -weight_coef * this->compute_hessian_coefficient(double_well_energy[blk], blk, blk, + std::span(u_values), + std::span(vaux_gf_at_ip)); AddMult_a_VVt(fun_val, Psi, *elmats(blk, blk)); } } @@ -350,8 +363,12 @@ void CahnHilliardNLFormIntegrator::AssembleElementGrad( el[blk]->CalcShape(ip, Psi); Tr.SetIntPoint(&ip); - double ww = ip.weight * Tr.Weight(); - AddMult_a_VVt(ww, Psi, *elmats(blk, off_blk)); + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + + AddMult_a_VVt(weight_coef, Psi, *elmats(blk, off_blk)); } } // Block 1 0 dR(mu)dphi=d(div M grad mu)/dphi @@ -401,10 +418,15 @@ void CahnHilliardNLFormIntegrator::AssembleElementGrad( u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip); } + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + const double coef_mob = this->compute_coefficient(mobility[blk], std::span(u_values), std::span(vaux_gf_at_ip)) * - ip.weight * Tr.Weight(); + weight_coef; el[blk]->CalcPhysDShape(Tr, gradPsi); AddMult_a_AAt(coef_mob, gradPsi, *elmats(blk, blk)); diff --git a/kernel/include/Integrators/DiffusionFluxNLFormIntegrator.hpp b/kernel/include/Integrators/DiffusionFluxNLFormIntegrator.hpp index 792f4045..0b3065e0 100644 --- a/kernel/include/Integrators/DiffusionFluxNLFormIntegrator.hpp +++ b/kernel/include/Integrators/DiffusionFluxNLFormIntegrator.hpp @@ -72,7 +72,7 @@ class DiffusionFluxNLFormIntegrator : public SlothNLFormIntegrator { void init() override; public: - DiffusionFluxNLFormIntegrator(const std::vector& u_old, + DiffusionFluxNLFormIntegrator(Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/DiffusionFluxNLFormIntegrator.tpp b/kernel/include/Integrators/DiffusionFluxNLFormIntegrator.tpp index 83f81919..b21f80a3 100644 --- a/kernel/include/Integrators/DiffusionFluxNLFormIntegrator.tpp +++ b/kernel/include/Integrators/DiffusionFluxNLFormIntegrator.tpp @@ -32,6 +32,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Coefficients/SlothBaseCoefficient.hpp" #include "Integrators/SlothGridFunction.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" @@ -80,10 +81,10 @@ void DiffusionFluxNLFormIntegrator::get_coefficients() { */ template DiffusionFluxNLFormIntegrator::DiffusionFluxNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : SlothNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : SlothNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->expected_list_.push_back(GlossaryType::Diffusivity); } @@ -186,8 +187,11 @@ void DiffusionFluxNLFormIntegrator::AssembleElementVector( // Diffusion flux (see child classes) this->add_diffusion_flux(Tr, nElement, ip, dim); - - this->Flux_ *= ip.weight * Tr.Weight(); + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + this->Flux_ *= weight_coef; this->gradPsi.AddMult(this->Flux_, *elvect[blk], 1.0); } @@ -248,10 +252,16 @@ void DiffusionFluxNLFormIntegrator::AssembleElementGrad( } Tr.SetIntPoint(&ip); + + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + const double coeff_diffu = this->compute_coefficient(stab_diffusion[blk], std::span(u_values), std::span(vaux_gf_at_ip)) * - ip.weight * Tr.Weight(); + weight_coef; el[blk]->CalcPhysDShape(Tr, gradPsi); AddMult_a_AAt(coeff_diffu, gradPsi, *elmat(blk, blk)); } diff --git a/kernel/include/Integrators/DiffusionNLFormIntegrator.hpp b/kernel/include/Integrators/DiffusionNLFormIntegrator.hpp index e84b0e4a..1adbcd3c 100644 --- a/kernel/include/Integrators/DiffusionNLFormIntegrator.hpp +++ b/kernel/include/Integrators/DiffusionNLFormIntegrator.hpp @@ -59,7 +59,7 @@ class DiffusionNLFormIntegrator : public SlothNLFormIntegrator { void init() override; public: - DiffusionNLFormIntegrator(const std::vector u_old, + DiffusionNLFormIntegrator(Geometry geometry, const std::vector u_old, const std::vector aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/DiffusionNLFormIntegrator.tpp b/kernel/include/Integrators/DiffusionNLFormIntegrator.tpp index 4ed7035a..8601ee8e 100644 --- a/kernel/include/Integrators/DiffusionNLFormIntegrator.tpp +++ b/kernel/include/Integrators/DiffusionNLFormIntegrator.tpp @@ -32,6 +32,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Coefficients/SlothBaseCoefficient.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" #include "Parameters/Parameter.hpp" @@ -60,10 +61,10 @@ */ template DiffusionNLFormIntegrator::DiffusionNLFormIntegrator( - const std::vector u_old, + Geometry geometry, const std::vector u_old, const std::vector aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : SlothNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) {} + : SlothNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) {} /** * @brief Initialize the diffusion integrator. @@ -155,7 +156,12 @@ void DiffusionNLFormIntegrator::AssembleElementVector( gradPsi.MultTranspose(*elfun[blk], gradU); double diffu = this->compute_coefficient(diffusion[blk], std::span(u_values), std::span(vaux_gf_at_ip)); - const double coeff_diffu = diffu * ip.weight * Tr.Weight(); + + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + const double coeff_diffu = diffu * weight_coef; gradU *= coeff_diffu; gradPsi.AddMult(gradU, *elvect[blk]); } @@ -227,11 +233,17 @@ void DiffusionNLFormIntegrator::AssembleElementGrad( std::span(vaux_gf_at_ip)); el[blk]->CalcPhysDShape(Tr, gradPsi); - AddMult_a_AAt(diffu * ip.weight * Tr.Weight(), gradPsi, *elmat(blk, blk)); + + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + + AddMult_a_AAt(diffu * weight_coef, gradPsi, *elmat(blk, blk)); gradPsi.MultTranspose(*elfun[blk], gradU); gradPsi.AddMult(gradU, vec); - AddMult_a_VWt(grad_diffu * ip.weight * Tr.Weight(), Psi, vec, *elmat(blk, blk)); + AddMult_a_VWt(grad_diffu * weight_coef, Psi, vec, *elmat(blk, blk)); } } } diff --git a/kernel/include/Integrators/FickNLFormIntegrator.hpp b/kernel/include/Integrators/FickNLFormIntegrator.hpp index 75c78f36..299346e5 100644 --- a/kernel/include/Integrators/FickNLFormIntegrator.hpp +++ b/kernel/include/Integrators/FickNLFormIntegrator.hpp @@ -56,7 +56,7 @@ class FickNLFormIntegrator : public DiffusionNLFormIntegrator { void get_coefficients() override; public: - FickNLFormIntegrator(const std::vector u_old, + FickNLFormIntegrator(Geometry geometry, const std::vector u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); virtual ~FickNLFormIntegrator() = default; diff --git a/kernel/include/Integrators/FickNLFormIntegrator.tpp b/kernel/include/Integrators/FickNLFormIntegrator.tpp index 23bcf8b1..f91120de 100644 --- a/kernel/include/Integrators/FickNLFormIntegrator.tpp +++ b/kernel/include/Integrators/FickNLFormIntegrator.tpp @@ -30,6 +30,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Integrators/DiffusionNLFormIntegrator.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" #include "Parameters/Parameter.hpp" @@ -57,12 +58,13 @@ * */ template -FickNLFormIntegrator::FickNLFormIntegrator(const std::vector u_old, +FickNLFormIntegrator::FickNLFormIntegrator(Geometry geometry, + const std::vector u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : DiffusionNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : DiffusionNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->integrator_name_ = "Fick"; this->expected_list_.push_back(GlossaryType::Diffusivity); diff --git a/kernel/include/Integrators/FourierNLFormIntegrator.hpp b/kernel/include/Integrators/FourierNLFormIntegrator.hpp index c06eb922..17a7ecad 100644 --- a/kernel/include/Integrators/FourierNLFormIntegrator.hpp +++ b/kernel/include/Integrators/FourierNLFormIntegrator.hpp @@ -56,7 +56,7 @@ class FourierNLFormIntegrator : public DiffusionNLFormIntegrator { void get_coefficients() override; public: - FourierNLFormIntegrator(const std::vector& u_old, + FourierNLFormIntegrator(Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/FourierNLFormIntegrator.tpp b/kernel/include/Integrators/FourierNLFormIntegrator.tpp index 53bf38ce..cf4c386e 100644 --- a/kernel/include/Integrators/FourierNLFormIntegrator.tpp +++ b/kernel/include/Integrators/FourierNLFormIntegrator.tpp @@ -30,6 +30,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Integrators/DiffusionNLFormIntegrator.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" #include "Parameters/Parameter.hpp" @@ -58,10 +59,10 @@ */ template FourierNLFormIntegrator::FourierNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : DiffusionNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : DiffusionNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->integrator_name_ = "Fourier"; this->expected_list_.push_back(GlossaryType::Conductivity); } diff --git a/kernel/include/Integrators/HeatTimeNLFormIntegrator.hpp b/kernel/include/Integrators/HeatTimeNLFormIntegrator.hpp index d8d6059b..55eea2c0 100644 --- a/kernel/include/Integrators/HeatTimeNLFormIntegrator.hpp +++ b/kernel/include/Integrators/HeatTimeNLFormIntegrator.hpp @@ -48,7 +48,7 @@ class HeatTimeNLFormIntegrator : public TimeNLFormIntegrator { void get_coefficients() override; public: - HeatTimeNLFormIntegrator(const std::vector u_old, + HeatTimeNLFormIntegrator(Geometry geometry, const std::vector u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/HeatTimeNLFormIntegrator.tpp b/kernel/include/Integrators/HeatTimeNLFormIntegrator.tpp index 8c565ffc..c365ca1c 100644 --- a/kernel/include/Integrators/HeatTimeNLFormIntegrator.tpp +++ b/kernel/include/Integrators/HeatTimeNLFormIntegrator.tpp @@ -29,6 +29,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" #include "Integrators/TimeNLFormIntegrator.hpp" #include "MAToolsProfiling/MATimersAPI.hxx" @@ -55,10 +56,10 @@ */ template HeatTimeNLFormIntegrator::HeatTimeNLFormIntegrator( - const std::vector u_old, + Geometry geometry, const std::vector u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : TimeNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : TimeNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->integrator_name_ = "HeatTime"; this->expected_list_.push_back(GlossaryType::Concentration); diff --git a/kernel/include/Integrators/LatentHeatNLFormIntegrator.hpp b/kernel/include/Integrators/LatentHeatNLFormIntegrator.hpp index 5abf3995..f3ea1e6f 100644 --- a/kernel/include/Integrators/LatentHeatNLFormIntegrator.hpp +++ b/kernel/include/Integrators/LatentHeatNLFormIntegrator.hpp @@ -74,7 +74,7 @@ class LatentHeatNLFormIntegrator : public SlothNLFormIntegrator { virtual void check_variables_consistency(); public: - LatentHeatNLFormIntegrator(const std::vector& u_old, + LatentHeatNLFormIntegrator(Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/LatentHeatNLFormIntegrator.tpp b/kernel/include/Integrators/LatentHeatNLFormIntegrator.tpp index 97225f84..a2ee98ab 100644 --- a/kernel/include/Integrators/LatentHeatNLFormIntegrator.tpp +++ b/kernel/include/Integrators/LatentHeatNLFormIntegrator.tpp @@ -33,6 +33,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" #include "MAToolsProfiling/MATimersAPI.hxx" #include "Parameters/Parameter.hpp" @@ -58,10 +59,10 @@ */ template LatentHeatNLFormIntegrator::LatentHeatNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : SlothNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : SlothNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->integrator_name_ = "LatentHeat"; this->latent_time_step_ = this->params_.template get_param_value("latent_time_step"); this->check_variables_consistency(); @@ -144,7 +145,7 @@ void LatentHeatNLFormIntegrator::AssembleElementVector( const mfem::Array& elvect) { int num_blocks = el.Size(); std::vector u_values(2 * num_blocks); - std::vector vaux_gf_at_ip(this->vaux_gf_.size()); + std::vector vaux_gf_at_ip(2 * this->nb_vaux_); for (int blk = 0; blk < num_blocks; ++blk) { // Catch_Time_Section("LatentHeatNLFormIntegrator:AssembleElementVector"); int nd = el[blk]->GetDof(); @@ -164,9 +165,9 @@ void LatentHeatNLFormIntegrator::AssembleElementVector( el[blk]->CalcShape(ip, Psi); // Tr.SetIntPoint(&ip); - for (size_t k = 0; k < this->vaux_gf_.size(); ++k) { + for (size_t k = 0; k < this->nb_vaux_; ++k) { vaux_gf_at_ip[k] = this->vaux_gf_[k].GetValue(Tr, ip); - vaux_gf_at_ip[k + num_blocks] = this->vaux_old_gf_[k].GetValue(Tr, ip); + vaux_gf_at_ip[k + this->nb_vaux_] = this->vaux_old_gf_[k].GetValue(Tr, ip); } // Get values for (int off_blk = 0; off_blk < num_blocks; ++off_blk) { @@ -174,10 +175,15 @@ void LatentHeatNLFormIntegrator::AssembleElementVector( u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip); } + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + const double latent_heat = this->get_latent_heat_at_ip(blk, std::span(u_values), std::span(vaux_gf_at_ip)) * - ip.weight * Tr.Weight(); + weight_coef; add(*elvect[blk], latent_heat, Psi, *elvect[blk]); } } @@ -259,8 +265,9 @@ template double LatentHeatNLFormIntegrator::get_latent_heat_at_ip( [[maybe_unused]] unsigned int blk, [[maybe_unused]] const std::span& values, [[maybe_unused]] const std::span& aux_values) { - std::span local_auxvalues(aux_values.begin(), aux_values.begin() + this->nb_blk_); - std::span local_auxvalues_n(aux_values.begin() + this->nb_blk_, aux_values.end()); + std::span local_auxvalues(aux_values.begin(), aux_values.begin() + this->nb_vaux_); + std::span local_auxvalues_n(aux_values.begin() + this->nb_vaux_, aux_values.end()); + const double mobility_value = this->get_mob_at_ip(blk, values, local_auxvalues); const double phi = local_auxvalues[this->phase_field_index_]; diff --git a/kernel/include/Integrators/MassDiffusionFluxNLFormIntegrator.hpp b/kernel/include/Integrators/MassDiffusionFluxNLFormIntegrator.hpp index e6adb99e..1f616f7d 100644 --- a/kernel/include/Integrators/MassDiffusionFluxNLFormIntegrator.hpp +++ b/kernel/include/Integrators/MassDiffusionFluxNLFormIntegrator.hpp @@ -82,7 +82,8 @@ class MassDiffusionFluxNLFormIntegrator : public DiffusionFluxNLFormIntegrator& u_old, + MassDiffusionFluxNLFormIntegrator(Geometry geometry, + const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/MassDiffusionFluxNLFormIntegrator.tpp b/kernel/include/Integrators/MassDiffusionFluxNLFormIntegrator.tpp index 67641ee1..aa788499 100644 --- a/kernel/include/Integrators/MassDiffusionFluxNLFormIntegrator.tpp +++ b/kernel/include/Integrators/MassDiffusionFluxNLFormIntegrator.tpp @@ -33,6 +33,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Integrators/SlothGridFunction.hpp" #include "Options/Options.hpp" #include "Parameters/Parameters.hpp" @@ -90,10 +91,10 @@ void MassDiffusionFluxNLFormIntegrator::get_parameters() { */ template MassDiffusionFluxNLFormIntegrator::MassDiffusionFluxNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : DiffusionFluxNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : DiffusionFluxNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->get_parameters(); this->check_variables_consistency(); } diff --git a/kernel/include/Integrators/MeltingBaseNLFormIntegrator.hpp b/kernel/include/Integrators/MeltingBaseNLFormIntegrator.hpp index 4088ed21..77406ca6 100644 --- a/kernel/include/Integrators/MeltingBaseNLFormIntegrator.hpp +++ b/kernel/include/Integrators/MeltingBaseNLFormIntegrator.hpp @@ -73,7 +73,7 @@ class MeltingBaseNLFormIntegrator : public SlothNLFormIntegrator { virtual void check_variables_consistency(); public: - MeltingBaseNLFormIntegrator(const std::vector& u_old, + MeltingBaseNLFormIntegrator(Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/MeltingBaseNLFormIntegrator.tpp b/kernel/include/Integrators/MeltingBaseNLFormIntegrator.tpp index d322c61f..0274c1a9 100644 --- a/kernel/include/Integrators/MeltingBaseNLFormIntegrator.tpp +++ b/kernel/include/Integrators/MeltingBaseNLFormIntegrator.tpp @@ -33,6 +33,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" #include "MAToolsProfiling/MATimersAPI.hxx" #include "Parameters/Parameter.hpp" @@ -58,10 +59,10 @@ */ template MeltingBaseNLFormIntegrator::MeltingBaseNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : SlothNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) {} + : SlothNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) {} /** * @brief Initialize the integrator. @@ -162,7 +163,13 @@ void MeltingBaseNLFormIntegrator::AssembleElementVector( std::span(u_values), std::span(vaux_gf_at_ip))) + seed; - ww *= ip.weight * Tr.Weight(); + + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + + ww *= weight_coef; add(*elvect[blk], ww, Psi, *elvect[blk]); } } @@ -233,7 +240,12 @@ void MeltingBaseNLFormIntegrator::AssembleElementGrad( interpolation_potential[blk], blk, blk, std::span(u_values), std::span(vaux_gf_at_ip)); // this->energy_derivatives(2, Tr, ip)(u); - fun_val *= ip.weight * Tr.Weight(); + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + + fun_val *= weight_coef; AddMult_a_VVt(fun_val, Psi, *elmats(blk, blk)); // w'(u)*(du, psi) } diff --git a/kernel/include/Integrators/MeltingCalphadNLFormIntegrator.hpp b/kernel/include/Integrators/MeltingCalphadNLFormIntegrator.hpp index cad16fba..1a80dde8 100644 --- a/kernel/include/Integrators/MeltingCalphadNLFormIntegrator.hpp +++ b/kernel/include/Integrators/MeltingCalphadNLFormIntegrator.hpp @@ -54,7 +54,8 @@ class MeltingCalphadNLFormIntegrator : public MeltingBaseNLFormIntegrator std::string secondary_phase_; int dgm_primary_phase_index_; int dgm_secondary_phase_index_; - int nucleus_index_; + int primary_nucleus_index_; + int secondary_nucleus_index_; double melting_temperature_; double melting_enthalpy_; @@ -72,7 +73,7 @@ class MeltingCalphadNLFormIntegrator : public MeltingBaseNLFormIntegrator const std::span& aux_values) override; public: - MeltingCalphadNLFormIntegrator(const std::vector& u_old, + MeltingCalphadNLFormIntegrator(Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/MeltingCalphadNLFormIntegrator.tpp b/kernel/include/Integrators/MeltingCalphadNLFormIntegrator.tpp index 6b259133..f3b7e856 100644 --- a/kernel/include/Integrators/MeltingCalphadNLFormIntegrator.tpp +++ b/kernel/include/Integrators/MeltingCalphadNLFormIntegrator.tpp @@ -31,6 +31,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Integrators/MeltingBaseNLFormIntegrator.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" #include "MAToolsProfiling/MATimersAPI.hxx" @@ -57,10 +58,10 @@ */ template MeltingCalphadNLFormIntegrator::MeltingCalphadNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : MeltingBaseNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : MeltingBaseNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->integrator_name_ = "MeltingCalphad"; this->get_parameters(); @@ -144,8 +145,10 @@ void MeltingCalphadNLFormIntegrator::check_driving_forces() { */ template void MeltingCalphadNLFormIntegrator::check_nucleus() { - bool nucleus_found = false; - this->nucleus_index_ = -1; + bool secondary_nucleus_found = false; + bool primary_nucleus_found = false; + this->secondary_nucleus_index_ = -1; + this->primary_nucleus_index_ = -1; for (std::size_t i = 0; i < this->aux_infos_.size(); ++i) { const auto& variable_info = this->aux_infos_[i]; MFEM_VERIFY(!variable_info.empty(), "Empty variable_info encountered."); @@ -164,11 +167,16 @@ void MeltingCalphadNLFormIntegrator::check_nucleus() { ": the name of the phase and the symbol 'nucleus'"); if (variable_info[0] == this->secondary_phase_) { - this->nucleus_index_ = i; - nucleus_found = true; + this->secondary_nucleus_index_ = i; + secondary_nucleus_found = true; + } + if (variable_info[0] == this->primary_phase_) { + this->primary_nucleus_index_ = i; + primary_nucleus_found = true; } } - MFEM_VERIFY(nucleus_found, "Nucleus for secondary phase must be set."); + MFEM_VERIFY(primary_nucleus_found, "Nucleus for primary phase must be set."); + MFEM_VERIFY(secondary_nucleus_found, "Nucleus for secondary phase must be set."); } /** @@ -203,7 +211,7 @@ double MeltingCalphadNLFormIntegrator::get_phase_change_at_ip( } /** - * @brief Compute the value of a seed of the secondary phase at integration point + * @brief Compute the value of a seed of the secondary or primary phases at integration point * * @tparam VARS Template parameter defining the variables used in the integrator. * @@ -221,10 +229,8 @@ double MeltingCalphadNLFormIntegrator::get_seed_at_ip( [[maybe_unused]] const std::span& aux_values) { // Nucleus must be equal to zero except when phase transition starts - // std::span u(values.begin(), values.begin() + this->nb_blk_); - // std::span un(values.begin() + this->nb_blk_, values.end()); - - const double seed = -aux_values[this->nucleus_index_]; + const double seed = + -aux_values[this->secondary_nucleus_index_] - aux_values[this->primary_nucleus_index_]; return seed; } diff --git a/kernel/include/Integrators/MeltingConstantNLFormIntegrator.hpp b/kernel/include/Integrators/MeltingConstantNLFormIntegrator.hpp index 3acc9395..0b558c0e 100644 --- a/kernel/include/Integrators/MeltingConstantNLFormIntegrator.hpp +++ b/kernel/include/Integrators/MeltingConstantNLFormIntegrator.hpp @@ -59,7 +59,8 @@ class MeltingConstantNLFormIntegrator : public MeltingBaseNLFormIntegrator const std::span& aux_values) override; public: - MeltingConstantNLFormIntegrator(const std::vector& u_old, + MeltingConstantNLFormIntegrator(Geometry geometry, + const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/MeltingConstantNLFormIntegrator.tpp b/kernel/include/Integrators/MeltingConstantNLFormIntegrator.tpp index a9fa1a9b..f2f2fab8 100644 --- a/kernel/include/Integrators/MeltingConstantNLFormIntegrator.tpp +++ b/kernel/include/Integrators/MeltingConstantNLFormIntegrator.tpp @@ -31,6 +31,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Integrators/MeltingBaseNLFormIntegrator.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" #include "MAToolsProfiling/MATimersAPI.hxx" @@ -57,10 +58,10 @@ */ template MeltingConstantNLFormIntegrator::MeltingConstantNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : MeltingBaseNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : MeltingBaseNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->integrator_name_ = "MeltingConstant"; this->get_parameters(); diff --git a/kernel/include/Integrators/MeltingTemperatureNLFormIntegrator.hpp b/kernel/include/Integrators/MeltingTemperatureNLFormIntegrator.hpp index f911417c..97ff8b0d 100644 --- a/kernel/include/Integrators/MeltingTemperatureNLFormIntegrator.hpp +++ b/kernel/include/Integrators/MeltingTemperatureNLFormIntegrator.hpp @@ -61,7 +61,8 @@ class MeltingTemperatureNLFormIntegrator : public MeltingBaseNLFormIntegrator& aux_values) override; public: - MeltingTemperatureNLFormIntegrator(const std::vector& u_old, + MeltingTemperatureNLFormIntegrator(Geometry geometry, + const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/MeltingTemperatureNLFormIntegrator.tpp b/kernel/include/Integrators/MeltingTemperatureNLFormIntegrator.tpp index dc21f3fb..fcae7344 100644 --- a/kernel/include/Integrators/MeltingTemperatureNLFormIntegrator.tpp +++ b/kernel/include/Integrators/MeltingTemperatureNLFormIntegrator.tpp @@ -57,10 +57,10 @@ */ template MeltingTemperatureNLFormIntegrator::MeltingTemperatureNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : MeltingBaseNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : MeltingBaseNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->integrator_name_ = "MeltingTemperature"; this->get_parameters(); this->check_variables_consistency(); diff --git a/kernel/include/Integrators/NeumannNLFormIntegrator.hpp b/kernel/include/Integrators/NeumannNLFormIntegrator.hpp index c2018fb8..1e3a96e1 100644 --- a/kernel/include/Integrators/NeumannNLFormIntegrator.hpp +++ b/kernel/include/Integrators/NeumannNLFormIntegrator.hpp @@ -48,7 +48,7 @@ class NeumannNLFormIntegrator : public RobinNLFormIntegrator { void get_coefficients() override; public: - NeumannNLFormIntegrator(const std::vector& u_old, + NeumannNLFormIntegrator(Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients, const unsigned int block, diff --git a/kernel/include/Integrators/NeumannNLFormIntegrator.tpp b/kernel/include/Integrators/NeumannNLFormIntegrator.tpp index 43f73875..790a1209 100644 --- a/kernel/include/Integrators/NeumannNLFormIntegrator.tpp +++ b/kernel/include/Integrators/NeumannNLFormIntegrator.tpp @@ -30,6 +30,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Integrators/NeumannNLFormIntegrator.hpp" #include "Integrators/RobinNLFormIntegrator.hpp" #include "Integrators/SlothGridFunction.hpp" @@ -78,10 +79,11 @@ void NeumannNLFormIntegrator::get_coefficients() { */ template NeumannNLFormIntegrator::NeumannNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients, const unsigned int block, const unsigned int bdr_id) - : RobinNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients, block, bdr_id) { + : RobinNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients, block, + bdr_id) { this->integrator_name_ = "Neumann"; } diff --git a/kernel/include/Integrators/RobinNLFormIntegrator.hpp b/kernel/include/Integrators/RobinNLFormIntegrator.hpp index 72d70c71..367290d1 100644 --- a/kernel/include/Integrators/RobinNLFormIntegrator.hpp +++ b/kernel/include/Integrators/RobinNLFormIntegrator.hpp @@ -61,7 +61,7 @@ class RobinNLFormIntegrator : public SlothNLFormIntegrator { void init() override; public: - RobinNLFormIntegrator(const std::vector& u_old, + RobinNLFormIntegrator(Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients, const unsigned int block, const unsigned int bdr_id); diff --git a/kernel/include/Integrators/RobinNLFormIntegrator.tpp b/kernel/include/Integrators/RobinNLFormIntegrator.tpp index 9d601736..3b21c756 100644 --- a/kernel/include/Integrators/RobinNLFormIntegrator.tpp +++ b/kernel/include/Integrators/RobinNLFormIntegrator.tpp @@ -32,6 +32,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Integrators/RobinNLFormIntegrator.hpp" #include "Integrators/SlothGridFunction.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" @@ -97,11 +98,11 @@ void RobinNLFormIntegrator::get_coefficients() { */ template RobinNLFormIntegrator::RobinNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients, const unsigned int block, const unsigned int bdr_id) - : SlothNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients), + : SlothNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients), blk_(block), bdr_id_(bdr_id), robin_a(Coefficient(Glossary::Default, 0.0)), @@ -179,14 +180,18 @@ void RobinNLFormIntegrator::AssembleElementVector( u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip); } + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } const double val = - Tr.Weight() * (this->compute_coefficient(robin_a, std::span(u_values), + weight_coef * (this->compute_coefficient(robin_a, std::span(u_values), std::span(vaux_gf_at_ip)) * u_values[blk] - this->compute_coefficient(robin_b, std::span(u_values), std::span(vaux_gf_at_ip))); - add(*elvect[blk], ip.weight * val, Psi, *elvect[blk]); + add(*elvect[blk], val, Psi, *elvect[blk]); } } } @@ -245,13 +250,18 @@ void RobinNLFormIntegrator::AssembleElementGrad( } Tr.SetIntPoint(&ip); + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + const double coeff_robin = (this->compute_coefficient(robin_a, std::span(u_values), std::span(vaux_gf_at_ip)) + this->compute_gradient_coefficient(robin_a, blk, std::span(u_values), std::span(vaux_gf_at_ip)) * u_values[blk]) * - ip.weight * Tr.Weight(); + weight_coef; AddMult_a_VVt(coeff_robin, Psi, *elmat(blk, blk)); } } diff --git a/kernel/include/Integrators/SlothNLFormIntegrator.hpp b/kernel/include/Integrators/SlothNLFormIntegrator.hpp index f9eec42c..36bde61f 100644 --- a/kernel/include/Integrators/SlothNLFormIntegrator.hpp +++ b/kernel/include/Integrators/SlothNLFormIntegrator.hpp @@ -39,6 +39,7 @@ #include "Coefficients/Coefficient.hpp" #include "Coefficients/Coefficients.hpp" #include "MAToolsProfiling/MATimersAPI.hxx" +#include "Options/ProblemsOptions.hpp" #include "Parameters/Parameters.hpp" #include "Utils/Utils.hpp" #include "mfem.hpp" // NOLINT [no include the directory when naming mfem include file] @@ -56,6 +57,8 @@ class SlothNLFormIntegrator : public mfem::BlockNonlinearFormIntegrator { std::vector> vect_aux_infos_; protected: + Geometry geometry_ = Geometry::Cartesian; + std::string integrator_name_ = ""; virtual void AssembleElementVector(const mfem::Array& el, mfem::ElementTransformation& Tr, @@ -77,6 +80,7 @@ class SlothNLFormIntegrator : public mfem::BlockNonlinearFormIntegrator { Parameters params_; std::vector coefficients_; unsigned int nb_blk_; + unsigned int nb_vaux_; std::vector get_aux_gf(); std::vector> get_aux_infos(); @@ -98,9 +102,11 @@ class SlothNLFormIntegrator : public mfem::BlockNonlinearFormIntegrator { const std::span& values, const std::span& aux_values); + bool isAxisymmetric() const noexcept; + public: virtual void init() = 0; - SlothNLFormIntegrator(const std::vector u_old, + SlothNLFormIntegrator(Geometry geometry, const std::vector u_old, const std::vector aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); virtual ~SlothNLFormIntegrator() = default; diff --git a/kernel/include/Integrators/SlothNLFormIntegrator.tpp b/kernel/include/Integrators/SlothNLFormIntegrator.tpp index 0abc07ed..74551875 100644 --- a/kernel/include/Integrators/SlothNLFormIntegrator.tpp +++ b/kernel/include/Integrators/SlothNLFormIntegrator.tpp @@ -35,9 +35,11 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Coefficients/Coefficient.hpp" #include "Coefficients/Coefficients.hpp" #include "MAToolsProfiling/MATimersAPI.hxx" +#include "Options/ProblemsOptions.hpp" #include "Parameters/Parameters.hpp" #include "Utils/Utils.hpp" #include "mfem.hpp" // NOLINT [no include the directory when naming mfem include file] @@ -50,14 +52,21 @@ * @param auxvars */ template -SlothNLFormIntegrator::SlothNLFormIntegrator(const std::vector u_old, +SlothNLFormIntegrator::SlothNLFormIntegrator(Geometry geometry, + const std::vector u_old, const std::vector aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : u_old_(u_old), aux_old_gf_(aux_old), params_(params), coefficients_(coefficients) { + : geometry_(geometry), + u_old_(u_old), + aux_old_gf_(aux_old), + params_(params), + coefficients_(coefficients) { this->nb_blk_ = this->u_old_.size(); + this->manage_auxiliary_variables(auxvars); + this->aux_gf_ = this->get_aux_gf(); this->aux_infos_ = this->get_aux_infos(); } @@ -80,6 +89,8 @@ void SlothNLFormIntegrator::manage_auxiliary_variables(std::vector this->vect_aux_infos_.emplace_back(std::move(var_info)); } } + + this->nb_vaux_ = this->vect_aux_gf_.size(); } /** @@ -122,7 +133,7 @@ std::vector> SlothNLFormIntegrator::get_aux_infos */ template void SlothNLFormIntegrator::check_coefficient_types(std::list expected_types) { - for (auto coefficients : this->coefficients_) { + for (const auto& coefficients : this->coefficients_) { auto vect_types = coefficients.get_types(); std::list TestedGlossaryType; TestedGlossaryType.assign(vect_types.begin(), vect_types.end()); @@ -160,10 +171,10 @@ std::optional SlothNLFormIntegrator::get_coefficient(const in GlossaryType type, unsigned int id, std::optional bdr_id) { - Coefficients coefficients = this->coefficients_[blk]; + const Coefficients& coefficients = this->coefficients_[blk]; for (unsigned int i = 0; i < coefficients.size(); i++) { - auto coef = coefficients[i]; + const auto& coef = coefficients[i]; if (coef.get_type() == type && coef.get_id() == id) { if (bdr_id.has_value()) { auto bdr_index = coef.get_bdr_index_coef(); @@ -319,3 +330,16 @@ double SlothNLFormIntegrator::compute_hessian_coefficient( } } } + +/** + * @brief Indicates whether the Problem uses an axisymmetric geometry. + * + * @tparam VAR Type representing the problem variables. + * + * @retval true The problem geometry is axisymmetric. + * @retval false The problem geometry is not axisymmetric. + */ +template +bool SlothNLFormIntegrator::isAxisymmetric() const noexcept { + return this->geometry_ == Geometry::Axisymmetric; +} \ No newline at end of file diff --git a/kernel/include/Integrators/ThermalDiffusionFluxNLFormIntegrator.hpp b/kernel/include/Integrators/ThermalDiffusionFluxNLFormIntegrator.hpp index 31f61728..062bf95d 100644 --- a/kernel/include/Integrators/ThermalDiffusionFluxNLFormIntegrator.hpp +++ b/kernel/include/Integrators/ThermalDiffusionFluxNLFormIntegrator.hpp @@ -58,7 +58,8 @@ class ThermalDiffusionFluxNLFormIntegrator : public DiffusionFluxNLFormIntegrato const mfem::IntegrationPoint& ip) override; public: - ThermalDiffusionFluxNLFormIntegrator(const std::vector& u_old, + ThermalDiffusionFluxNLFormIntegrator(Geometry geometry, + const std::vector& u_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); }; diff --git a/kernel/include/Integrators/ThermalDiffusionFluxNLFormIntegrator.tpp b/kernel/include/Integrators/ThermalDiffusionFluxNLFormIntegrator.tpp index 4bf32809..69aa8aec 100644 --- a/kernel/include/Integrators/ThermalDiffusionFluxNLFormIntegrator.tpp +++ b/kernel/include/Integrators/ThermalDiffusionFluxNLFormIntegrator.tpp @@ -28,6 +28,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Integrators/SlothGridFunction.hpp" #include "Options/Options.hpp" #include "Parameters/Parameters.hpp" @@ -62,9 +63,9 @@ void ThermalDiffusionFluxNLFormIntegrator::get_parameters() { */ template ThermalDiffusionFluxNLFormIntegrator::ThermalDiffusionFluxNLFormIntegrator( - const std::vector& u_old, const Parameters& params, + Geometry geometry, const std::vector& u_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : DiffusionFluxNLFormIntegrator(u_old, params, auxvars, coefficients) { + : DiffusionFluxNLFormIntegrator(geometry, u_old, params, auxvars, coefficients) { this->check_variables_consistency(); } diff --git a/kernel/include/Integrators/TimeCHNLFormIntegrator.hpp b/kernel/include/Integrators/TimeCHNLFormIntegrator.hpp index 60fdfc71..a73f442c 100644 --- a/kernel/include/Integrators/TimeCHNLFormIntegrator.hpp +++ b/kernel/include/Integrators/TimeCHNLFormIntegrator.hpp @@ -61,7 +61,7 @@ class TimeCHNLFormIntegrator : public SlothNLFormIntegrator { public: void init() override; - TimeCHNLFormIntegrator(const std::vector& u_old, + TimeCHNLFormIntegrator(Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/TimeCHNLFormIntegrator.tpp b/kernel/include/Integrators/TimeCHNLFormIntegrator.tpp index 8c90f45f..51e93ee7 100644 --- a/kernel/include/Integrators/TimeCHNLFormIntegrator.tpp +++ b/kernel/include/Integrators/TimeCHNLFormIntegrator.tpp @@ -31,6 +31,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" #include "Integrators/TimeCHNLFormIntegrator.hpp" #include "MAToolsProfiling/MATimersAPI.hxx" @@ -57,10 +58,10 @@ */ template TimeCHNLFormIntegrator::TimeCHNLFormIntegrator( - const std::vector& u_old, + Geometry geometry, const std::vector& u_old, const std::vector& aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : SlothNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : SlothNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->integrator_name_ = "SplitTimeDerivative"; this->check_variables_consistency(); @@ -224,10 +225,14 @@ void TimeCHNLFormIntegrator::AssembleElementVector( u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip); } + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } double coef_a = this->compute_coefficient(this->coefficient_A[blk], std::span(u_values), std::span(vaux_gf_at_ip)); - const double ww = coef_a * phi * ip.weight * Tr.Weight(); + const double ww = coef_a * phi * weight_coef; add(*elvect[blk], ww, Psi, *elvect[blk]); } } @@ -322,10 +327,15 @@ void TimeCHNLFormIntegrator::AssembleElementGrad( u_values[off_blk + num_blocks] = this->u_old_[off_blk].GetValue(Tr, ip); } + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + double coef_a = this->compute_coefficient(this->coefficient_A[blk], std::span(u_values), std::span(vaux_gf_at_ip)); - double w = coef_a * Tr.Weight() * ip.weight; + double w = coef_a * weight_coef; AddMult_a_VVt(w, Psi, *elmats(blk, off_blk)); } } diff --git a/kernel/include/Integrators/TimeNLFormIntegrator.hpp b/kernel/include/Integrators/TimeNLFormIntegrator.hpp index dc70bc15..72d97734 100644 --- a/kernel/include/Integrators/TimeNLFormIntegrator.hpp +++ b/kernel/include/Integrators/TimeNLFormIntegrator.hpp @@ -62,7 +62,7 @@ class TimeNLFormIntegrator : public SlothNLFormIntegrator { public: void init() override; - TimeNLFormIntegrator(const std::vector u_old, + TimeNLFormIntegrator(Geometry geometry, const std::vector u_old, const std::vector aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients); diff --git a/kernel/include/Integrators/TimeNLFormIntegrator.tpp b/kernel/include/Integrators/TimeNLFormIntegrator.tpp index c7be183e..05210569 100644 --- a/kernel/include/Integrators/TimeNLFormIntegrator.tpp +++ b/kernel/include/Integrators/TimeNLFormIntegrator.tpp @@ -31,6 +31,7 @@ #include #include +#include "Coefficients/AxiCylindricalCoefficient.hpp" #include "Integrators/SlothNLFormIntegrator.hpp" #include "MAToolsProfiling/MATimersAPI.hxx" #include "Parameters/Parameters.hpp" @@ -55,12 +56,13 @@ * */ template -TimeNLFormIntegrator::TimeNLFormIntegrator(const std::vector u_old, +TimeNLFormIntegrator::TimeNLFormIntegrator(Geometry geometry, + const std::vector u_old, const std::vector aux_old, const Parameters& params, std::vector auxvars, const std::vector& coefficients) - : SlothNLFormIntegrator(u_old, aux_old, params, auxvars, coefficients) { + : SlothNLFormIntegrator(geometry, u_old, aux_old, params, auxvars, coefficients) { this->integrator_name_ = "TimeDerivative"; this->check_variables_consistency(); @@ -229,7 +231,12 @@ void TimeNLFormIntegrator::AssembleElementVector( double coef_b = this->compute_coefficient(this->coefficient_B[blk], std::span(u_values), std::span(vaux_gf_at_ip)); - const double ww = coef_a * coef_b * u * ip.weight * Tr.Weight(); + + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + const double ww = coef_a * coef_b * u * weight_coef; add(*elvect[blk], ww, Psi, *elvect[blk]); } } @@ -295,8 +302,11 @@ void TimeNLFormIntegrator::AssembleElementGrad( double coef_b = this->compute_coefficient(this->coefficient_B[blk], std::span(u_values), std::span(vaux_gf_at_ip)); - - double fun_val = coef_a * coef_b * ip.weight * Tr.Weight(); + double weight_coef = ip.weight * Tr.Weight(); + if (this->isAxisymmetric()) { + weight_coef *= AxiCylindricalCoefficient().Eval(Tr, ip); + } + double fun_val = coef_a * coef_b * weight_coef; AddMult_a_VVt(fun_val, Psi, *elmats(blk, blk)); } } diff --git a/kernel/include/Operators/OperatorBase.hpp b/kernel/include/Operators/OperatorBase.hpp index 8fbb824c..36cbd8d0 100644 --- a/kernel/include/Operators/OperatorBase.hpp +++ b/kernel/include/Operators/OperatorBase.hpp @@ -69,6 +69,7 @@ class OperatorBase : public mfem::Operator { void set_default_solver(); protected: + Geometry geometry_ = Geometry::Cartesian; Parameters nl_solver_params_; std::optional energy_coefficient_; std::optional grad_energy_coefficient_; @@ -94,7 +95,7 @@ class OperatorBase : public mfem::Operator { /// Boundary conditions std::vector*> bcs_; - mfem::Array array_bdr_; + std::vector> array_bdr_; std::vector> src_func_; @@ -173,6 +174,7 @@ class OperatorBase : public mfem::Operator { void overload_preconditioner(VSolverType PRECOND); void overload_preconditioner(VSolverType PRECOND, const Parameters& p_params); + void setGeometry(Geometry geometry); // Virtual methods virtual void initialize(const double& initial_time, Variables& vars, std::vector*> auxvars); diff --git a/kernel/include/Operators/OperatorBase.tpp b/kernel/include/Operators/OperatorBase.tpp index 2528cf35..db5f6a8d 100644 --- a/kernel/include/Operators/OperatorBase.tpp +++ b/kernel/include/Operators/OperatorBase.tpp @@ -406,25 +406,39 @@ void OperatorBase::build_rhs_nonlinear_form(const std::vectorblock_trueOffsets_.Size() - 1; mfem::Array Robin_bdr; mfem::Array Neumann_bdr; + + auto make_boundary_marker = [](int nb_bdr, int j) { + mfem::Array marker(nb_bdr); + marker = 0; + marker[j] = 1; + return marker; + }; + // Loop over variables for (int i = 0; i < fes_size; ++i) { Robin_bdr = this->bcs_[i]->get_marker_array("Robin"); + const int nb_bdr = Robin_bdr.Size(); Neumann_bdr = this->bcs_[i]->get_marker_array("Neumann"); - const int nb_bdr = Robin_bdr.Size(); - this->array_bdr_.SetSize(nb_bdr); - Coefficients coefficients = this->coefficients_[i]; + // All finite element spaces are assumed to live on the same mesh. + // Therefore all boundary marker arrays have the same size. + if (this->array_bdr_.empty()) { + this->array_bdr_.resize(nb_bdr); + } + MFEM_VERIFY(this->array_bdr_.size() == nb_bdr, "Inconsistent number of boundary attributes."); + + const Coefficients& coefficients = this->coefficients_[i]; const int coef_size = coefficients.size(); // Loop over boundaries - for (auto j = 0; j < nb_bdr; j++) { + for (int j = 0; j < nb_bdr; j++) { // Add Neumann integrator if (Neumann_bdr[j] > 0) { // Check if a coefficient is given for this bc // (else Homogeneous Neumann) bool has_neumann_coeff = false; for (int l = 0; l < coef_size; l++) { - auto coef = coefficients[l]; + const auto& coef = coefficients[l]; if (coef.get_type() == GlossaryType::Neumann) { auto bdr_ids = coef.get_bdr_index_coef(); if (std::find(bdr_ids.begin(), bdr_ids.end(), j) != bdr_ids.end()) { @@ -436,18 +450,16 @@ void OperatorBase::build_rhs_nonlinear_form(const std::vectorarray_bdr_ = 0; - this->array_bdr_[j] = 1; + this->array_bdr_[j] = make_boundary_marker(nb_bdr, j); auto integrator_ptr = this->set_bdr_nlfi_ptr("Neumann", u_vect, i, j); - this->RHS->AddBoundaryIntegrator(integrator_ptr, this->array_bdr_); + this->RHS->AddBoundaryIntegrator(integrator_ptr, this->array_bdr_[j]); } } // Add Robin integrator if (Robin_bdr[j] > 0) { - this->array_bdr_ = 0; - this->array_bdr_[j] = 1; + this->array_bdr_[j] = make_boundary_marker(nb_bdr, j); auto integrator_ptr = this->set_bdr_nlfi_ptr("Robin", u_vect, i, j); - this->RHS->AddBoundaryIntegrator(integrator_ptr, this->array_bdr_); + this->RHS->AddBoundaryIntegrator(integrator_ptr, this->array_bdr_[j]); } } } @@ -554,7 +566,7 @@ void OperatorBase::ComputeIntegral(const int& it, const double& t, const fe = this->fes_[id_var]->GetFE(i); const mfem::IntegrationRule* ir; - int intorder = 2 * fe->GetOrder() + 3; // <---------- + int intorder = 2 * fe->GetOrder() + 3; ir = &(mfem::IntRules.Get(fe->GetGeomType(), intorder)); mfem::real_t int_elem = 0.0; @@ -911,43 +923,43 @@ SlothNLFormIntegrator>* OperatorBase::get_rhs_integrat switch (Integrators::from(integrator)) { case Integrators::MassFlux: { return new MassDiffusionFluxNLFormIntegrator>( - vun, vauxn, all_params, this->auxvariables_, this->coefficients_); + this->geometry_, vun, vauxn, all_params, this->auxvariables_, this->coefficients_); } case Integrators::Fick: { - return new FickNLFormIntegrator>(vun, vauxn, all_params, + return new FickNLFormIntegrator>(this->geometry_, vun, vauxn, all_params, this->auxvariables_, this->coefficients_); } case Integrators::Fourier: { return new FourierNLFormIntegrator>( - vun, vauxn, all_params, this->auxvariables_, this->coefficients_); + this->geometry_, vun, vauxn, all_params, this->auxvariables_, this->coefficients_); } case Integrators::CahnHilliard: { return new CahnHilliardNLFormIntegrator>( - vun, vauxn, all_params, this->auxvariables_, this->coefficients_); + this->geometry_, vun, vauxn, all_params, this->auxvariables_, this->coefficients_); } case Integrators::AllenCahn: { return new AllenCahnNLFormIntegrator>( - vun, vauxn, all_params, this->auxvariables_, this->coefficients_); + this->geometry_, vun, vauxn, all_params, this->auxvariables_, this->coefficients_); } case Integrators::SplitAllenCahn: { return new BlockAllenCahnNLFormIntegrator>( - vun, vauxn, all_params, this->auxvariables_, this->coefficients_); + this->geometry_, vun, vauxn, all_params, this->auxvariables_, this->coefficients_); } case Integrators::MeltingTemperature: { return new MeltingTemperatureNLFormIntegrator>( - vun, vauxn, all_params, this->auxvariables_, this->coefficients_); + this->geometry_, vun, vauxn, all_params, this->auxvariables_, this->coefficients_); } case Integrators::MeltingCalphad: { return new MeltingCalphadNLFormIntegrator>( - vun, vauxn, all_params, this->auxvariables_, this->coefficients_); + this->geometry_, vun, vauxn, all_params, this->auxvariables_, this->coefficients_); } case Integrators::MeltingConstant: { return new MeltingConstantNLFormIntegrator>( - vun, vauxn, all_params, this->auxvariables_, this->coefficients_); + this->geometry_, vun, vauxn, all_params, this->auxvariables_, this->coefficients_); } case Integrators::LatentHeat: { return new LatentHeatNLFormIntegrator>( - vun, vauxn, all_params, this->auxvariables_, this->coefficients_); + this->geometry_, vun, vauxn, all_params, this->auxvariables_, this->coefficients_); } default: mfem::mfem_error("RHS Integrators not found. Please check your data."); @@ -974,13 +986,15 @@ SlothNLFormIntegrator>* OperatorBase::get_bdr_integrat const unsigned int block, const unsigned int bdr_id) { switch (Integrators::from(integrator)) { case Integrators::Neumann: { - return new NeumannNLFormIntegrator>( - vun, vauxn, all_params, this->auxvariables_, this->coefficients_, block, bdr_id); + return new NeumannNLFormIntegrator>(this->geometry_, vun, vauxn, all_params, + this->auxvariables_, + this->coefficients_, block, bdr_id); break; } case Integrators::Robin: { - return new RobinNLFormIntegrator>( - vun, vauxn, all_params, this->auxvariables_, this->coefficients_, block, bdr_id); + return new RobinNLFormIntegrator>(this->geometry_, vun, vauxn, all_params, + this->auxvariables_, this->coefficients_, + block, bdr_id); break; } default: @@ -1090,3 +1104,15 @@ void OperatorBase::set_time_coefficients(double time) { (*this->grad_energy_coefficient_).set_time(time); } } + +/** + * @brief Define the geometry of the problem (Cartesian, Axisymmetric) + * + * @tparam T Finite Element collection (mfem object) + * @tparam DIM Spatial dimension + * @param geometry + */ +template +void OperatorBase::setGeometry(Geometry geometry) { + this->geometry_ = geometry; +} diff --git a/kernel/include/Operators/TransientOperator.tpp b/kernel/include/Operators/TransientOperator.tpp index 8e655bcf..213e2e1e 100644 --- a/kernel/include/Operators/TransientOperator.tpp +++ b/kernel/include/Operators/TransientOperator.tpp @@ -660,18 +660,18 @@ SlothNLFormIntegrator>* TransientOperator::get_lhs_int const std::vector& vauxn, const Parameters& all_params) { switch (Integrators::from(integrator)) { case Integrators::TimeDerivative: { - return new TimeNLFormIntegrator>(vun, vauxn, all_params, + return new TimeNLFormIntegrator>(this->geometry_, vun, vauxn, all_params, this->auxvariables_, this->coefficients_); break; } case Integrators::HeatTimeDerivative: { return new HeatTimeNLFormIntegrator>( - vun, vauxn, all_params, this->auxvariables_, this->coefficients_); + this->geometry_, vun, vauxn, all_params, this->auxvariables_, this->coefficients_); break; } case Integrators::SplitTimeDerivative: { return new TimeCHNLFormIntegrator>( - vun, vauxn, all_params, this->auxvariables_, this->coefficients_); + this->geometry_, vun, vauxn, all_params, this->auxvariables_, this->coefficients_); break; } default: diff --git a/kernel/include/Options/ProblemsOptions.hpp b/kernel/include/Options/ProblemsOptions.hpp index 2cdc1a31..9521d4a8 100644 --- a/kernel/include/Options/ProblemsOptions.hpp +++ b/kernel/include/Options/ProblemsOptions.hpp @@ -37,3 +37,5 @@ struct Problems { enum value { Diffusion, AllenCahn, Calphad }; static value from(const std::string&); }; + +enum class Geometry { Cartesian, Axisymmetric }; diff --git a/kernel/include/Problems/Problem.tpp b/kernel/include/Problems/Problem.tpp index 68118b52..2afca411 100644 --- a/kernel/include/Problems/Problem.tpp +++ b/kernel/include/Problems/Problem.tpp @@ -394,6 +394,7 @@ void Problem::do_time_step( this->set_time_coefficients(next_time); + this->oper_.setGeometry(this->geometry_); this->oper_.solve(vect_unk, next_time, current_time, current_time_step, iter); // Store the solution into a temporary mfem::Vector that will be used during updating stage, if diff --git a/kernel/include/Problems/ProblemBase.hpp b/kernel/include/Problems/ProblemBase.hpp index 652b5350..4f4b3e36 100644 --- a/kernel/include/Problems/ProblemBase.hpp +++ b/kernel/include/Problems/ProblemBase.hpp @@ -38,6 +38,7 @@ #include "Coefficients/Coefficients.hpp" #include "Convergence/PhysicalConvergence.hpp" +#include "Options/ProblemsOptions.hpp" #include "Parameters/Parameter.hpp" #include "PostProcessing/postprocessing.hpp" #include "Variables/Variable.hpp" @@ -64,6 +65,7 @@ class ProblemBase { void check_convergence(const std::vector>& unks); protected: + Geometry geometry_ = Geometry::Cartesian; std::string name_{"Unnamed problem"}; VAR& variables_; std::vector auxvariables_; @@ -98,6 +100,7 @@ class ProblemBase { std::vector> get_convergence(); + void setGeometry(Geometry geometry); ///////////////////////////////////////////////////// virtual void initialize([[maybe_unused]] const double& initial_time) {} diff --git a/kernel/include/Problems/ProblemBase.tpp b/kernel/include/Problems/ProblemBase.tpp index a46e8e13..ab04151b 100644 --- a/kernel/include/Problems/ProblemBase.tpp +++ b/kernel/include/Problems/ProblemBase.tpp @@ -38,6 +38,7 @@ #include "Coefficients/Coefficients.hpp" #include "Convergence/PhysicalConvergence.hpp" +#include "Options/ProblemsOptions.hpp" #include "Parameters/Parameter.hpp" #include "PostProcessing/postprocessing.hpp" #include "Variables/Variable.hpp" @@ -345,3 +346,15 @@ template std::vector> ProblemBase::get_convergence() { return this->var_convergence_; } + +/** + * @brief Define the geometry of the problem (Cartesian, Axisymmetric) + * + * @tparam VAR Type representing the problem Variables. + * @tparam PST Type representing the post-processing. + * @param geometry + */ +template +void ProblemBase::setGeometry(Geometry geometry) { + this->geometry_ = geometry; +} diff --git a/kernel/src/Coefficients/AxiCylindricalCoefficient.cpp b/kernel/src/Coefficients/AxiCylindricalCoefficient.cpp new file mode 100644 index 00000000..15abecb7 --- /dev/null +++ b/kernel/src/Coefficients/AxiCylindricalCoefficient.cpp @@ -0,0 +1,52 @@ +/** + * @file AxiCylindricalCoefficient.cpp + * @author Clément Introïni (clement.introini@cea.fr) + * @brief Built a MFEM Coefficient for AxiCylindricalCoefficient + * @version 0.1 + * @date 2026-06-16 + * + * @copyright CEA (C) 2026 + * + * This file is part of SLOTH. + * + * SLOTH is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * SLOTH is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see . + * + */ +#include "Coefficients/AxiCylindricalCoefficient.hpp" + +#include +#include +#include +#include + +#include "mfem.hpp" // NOLINT [no include the directory when naming mfem include file] + +/** + * @brief Construct a new AxiCylindricalCoefficient::AxiCylindricalCoefficient object + * + */ +AxiCylindricalCoefficient::AxiCylindricalCoefficient() : transip(3) {} + +/** + * @brief Eval AxiCylindricalCoefficient at the integration point + * + * @param T + * @param ip + * @return real_t + */ +mfem::real_t AxiCylindricalCoefficient::Eval(mfem::ElementTransformation& T, + const mfem::IntegrationPoint& ip) { + T.Transform(ip, transip); + return abs(transip[0]); +} diff --git a/kernel/src/Coefficients/Coefficients.cpp b/kernel/src/Coefficients/Coefficients.cpp index f61005a5..da4eaddb 100644 --- a/kernel/src/Coefficients/Coefficients.cpp +++ b/kernel/src/Coefficients/Coefficients.cpp @@ -60,7 +60,7 @@ std::vector Coefficients::getCoefficients() const { return this->ve * * @return Number of stored coefficients. */ -size_t Coefficients::size() noexcept { return this->vect_coefficients_.size(); } +size_t Coefficients::size() const noexcept { return this->vect_coefficients_.size(); } /** * @brief Returns the i-th coefficient. @@ -70,7 +70,20 @@ size_t Coefficients::size() noexcept { return this->vect_coefficients_.size(); } * * @throws std::out_of_range if i is out of bounds. */ -Coefficient Coefficients::operator[](size_t i) { +Coefficient& Coefficients::operator[](size_t i) { + if (i >= vect_coefficients_.size()) throw std::out_of_range("Index out of range"); + return this->vect_coefficients_[i]; +} + +/** + * @brief Returns the i-th coefficient. + * + * @param i Index of the coefficient. + * @return Reference to the i-th coefficient. + * + * @throws std::out_of_range if i is out of bounds. + */ +const Coefficient& Coefficients::operator[](size_t i) const { if (i >= vect_coefficients_.size()) throw std::out_of_range("Index out of range"); return this->vect_coefficients_[i]; } diff --git a/kernel/src/Coefficients/SlothBaseCoefficient.cpp b/kernel/src/Coefficients/SlothBaseCoefficient.cpp index c535ebe2..13166b4a 100644 --- a/kernel/src/Coefficients/SlothBaseCoefficient.cpp +++ b/kernel/src/Coefficients/SlothBaseCoefficient.cpp @@ -362,4 +362,4 @@ void SlothBaseCoefficient::set_bdr_index_coef(std::vector ids) { this->bdr_ * * @return Vector of boundary ids where the coefficient is applied. */ -std::vector SlothBaseCoefficient::get_bdr_index_coef() { return this->bdr_index_; } +std::vector SlothBaseCoefficient::get_bdr_index_coef() const { return this->bdr_index_; } diff --git a/tests/HeatTransfer/1D/CMakeLists.txt b/tests/HeatTransfer/1D/CMakeLists.txt index b1da7574..7e029feb 100644 --- a/tests/HeatTransfer/1D/CMakeLists.txt +++ b/tests/HeatTransfer/1D/CMakeLists.txt @@ -1,2 +1,3 @@ add_subdirectory(test1) add_subdirectory(test2) +add_subdirectory(test3) diff --git a/tests/HeatTransfer/1D/test3/CMakeLists.txt b/tests/HeatTransfer/1D/test3/CMakeLists.txt new file mode 100644 index 00000000..3b209c68 --- /dev/null +++ b/tests/HeatTransfer/1D/test3/CMakeLists.txt @@ -0,0 +1,6 @@ +create_test("HeatTransfer1Dtest3" "HeatTransfer1Dtest3" FALSE "2D EXE HEAT" 1) +configure_file(${CMAKE_SOURCE_DIR}/tests/tools/convergence_study.py ${CMAKE_CURRENT_BINARY_DIR}/convergence_study.py COPYONLY) +create_col_comparison_convergence("CompareHeatTransfer1Dtest3Convergence" "convergence_output_ref.csv" "convergence_output.csv" -1 absolute 1e-16 FALSE "HeatTransfer1Dtest3" "1D HEAT" 0.00465 100000) + + + diff --git a/tests/HeatTransfer/1D/test3/main.cpp b/tests/HeatTransfer/1D/test3/main.cpp new file mode 100644 index 00000000..9e780a6d --- /dev/null +++ b/tests/HeatTransfer/1D/test3/main.cpp @@ -0,0 +1,181 @@ +/** + * @file main.cpp + * @author ci230846 (clement.introini@cea.fr) + * @brief 2D axisymmetric heat transfer problem + * @version 0.1 + * @date 2024-09-3 + * + * @copyright Copyright (c) 2024 + * + */ +#include +#include +#include +#include +#include +#include + +#include "Sloth/sloth.hpp" +#include "Sloth/tests.hpp" +///--------------- +/// Main program +///--------------- +int main(int argc, char* argv[]) { + //--------------------------------------- + // Initialize MPI and HYPRE + //--------------------------------------- + setVerbosity(Verbosity::Debug); + mfem::Mpi::Init(argc, argv); + mfem::Hypre::Init(); + // + //--------------------------------------- + // Profiling + Profiling::getInstance().enable(); + //--------------------------------------- + ///////////////////////// + const int DIM = 1; + using FECollection = Test::FECollection; + using VARS = Test::VARS; + using VAR = Test::VAR; + using PST = Test::PST; + using SPA = Test::SPA; + using BCS = Test::BCS; + // + using OPE = TransientOperator; + using PB = Problem; + + // ########################################### + // ########################################### + // Spatial Discretization // + // ########################################### + // ########################################### + // ############################## + // Meshing // + // ############################## + const std::string& mesh_type = "InlineLineWithSegments"; // type of mesh + const int refinement_level = 0; // number of levels of uniform refinement + + const double pellet_radius = 0.00465; + std::vector vect_NN{30, 60, 90, 120}; + for (const auto NN : vect_NN) { + const int order = 1; // finite element order + const int n_elements = NN; + + SPA spatial(mesh_type, order, refinement_level, std::make_tuple(n_elements, pellet_radius)); + + // ############################## + // Boundary conditions // + // ############################## + const double T_initial = 750.0; + auto boundaries = {Boundary("left", 0, "Neumann"), + Boundary("right", 1, "Dirichlet", T_initial)}; + auto Tbcs = BCS(&spatial, boundaries); + // ########################################### + // ########################################### + // Physical models // + // ########################################### + // ########################################### + // #################### + // parameters // + // #################### + // Heat + const auto& rho(1.e3); + const auto& cp(10.); + const auto& cond(2.); + Coefficient density(Glossary::Concentration, rho); + Coefficient heat_capacity(Glossary::Cp, cp); + Coefficient conductivity(Glossary::Conductivity, cond); + + // ############################ + // variables IC + SRC // + // ########################### + // Heat + auto pl = 26.e3; + + auto src_func = std::function( + [pl, pellet_radius]([[maybe_unused]] const mfem::Vector& vcoord, + [[maybe_unused]] double time) { + const double rr = vcoord[0]; + const double func = pl * rr / (M_PI * pellet_radius * pellet_radius); + + return func; + }); + + auto T_analytical = std::function( + [pl, pellet_radius, T_initial, cond]([[maybe_unused]] const mfem::Vector& vcoord, + [[maybe_unused]] double time) { + const double rr = vcoord[0] * vcoord[0]; + const double func = T_initial + pl * (pellet_radius * pellet_radius - rr) / + (4.0 * M_PI * pellet_radius * pellet_radius * cond); + + return func; + }); + + auto heat_vars = VARS(VAR(&spatial, Tbcs, "T", Glossary::Temperature, 2, T_initial, + AnalyticalFunctions(T_analytical))); + // ########################################### + // ########################################### + // Post-processing // + // ########################################### + // ########################################### + + const std::string& main_folder_path = + "Saves_order_" + std::to_string(order) + "_Nx" + std::to_string(NN); + const auto& level_of_detail = 1; + const auto& frequency = 1; + // Heat + const std::string& calculation_path = "Problem1"; + auto p_pst = Parameters( + Parameter("main_folder_path", main_folder_path), + Parameter("calculation_path", calculation_path), Parameter("frequency", frequency), + Parameter("level_of_detail", level_of_detail), Parameter("enable_compute_energies", false)); + auto pst = PST(&spatial, p_pst); + + // #################### + // Problems // + // #################### + + // Heat: + Coefficients coef_pb(density, heat_capacity, conductivity); + std::vector > src_term; + src_term.emplace_back(AnalyticalFunctions(src_func)); + std::vector spatials{&spatial}; + OPE oper(spatials, {"Fourier"}, TimeScheme::EulerImplicit, "HeatTimeDerivative", src_term); + + oper.overload_nl_solver( + NLSolverType::NEWTON, + Parameters(Parameter("description", "Newton solver "), Parameter("print_level", 1), + Parameter("rel_tol", 1.e-9), Parameter("abs_tol", 1.e-9))); + + auto T_cvg = PhysicalConvergence(ConvergenceType::ABSOLUTE_MAX, 1.e-16); + auto CVG = Convergence(T_cvg); + + PB Heat_pb("Heat", oper, heat_vars, {coef_pb}, CVG, pst); + Heat_pb.setGeometry(Geometry::Axisymmetric); + // Coupling 1 + auto cc = Coupling("Heat transfer", Heat_pb); + // ########################################### + // ########################################### + // Time-integration // + // ########################################### + // ########################################### + const auto& t_initial = 0.0; + const auto& t_final = 10.; + const auto& dt = 0.01; + auto time_params = Parameters(Parameter("initial_time", t_initial), + Parameter("final_time", t_final), Parameter("time_step", dt)); + auto time = TimeDiscretization(time_params, cc); + + time.solve(); + //--------------------------------------- + // Profiling stop + //--------------------------------------- + Profiling::getInstance().print(); + } + //--------------------------------------- + // Finalize MPI + //--------------------------------------- + MPI_Finalize(); + //--------------------------------------- + return 0; +} diff --git a/tests/HeatTransfer/1D/test3/ref/convergence_output_ref.csv b/tests/HeatTransfer/1D/test3/ref/convergence_output_ref.csv new file mode 100644 index 00000000..c9e27e17 --- /dev/null +++ b/tests/HeatTransfer/1D/test3/ref/convergence_output_ref.csv @@ -0,0 +1,4 @@ +order 1,3.875e-05,0.000891857 +order 1,5.1666666666666664e-05,0.00158401 +order 1,7.75e-05,0.00355722 +order 1,0.000155,0.0141473 diff --git a/tests/HeatTransfer/2D/CMakeLists.txt b/tests/HeatTransfer/2D/CMakeLists.txt index b1da7574..03456f5e 100644 --- a/tests/HeatTransfer/2D/CMakeLists.txt +++ b/tests/HeatTransfer/2D/CMakeLists.txt @@ -1,2 +1,4 @@ add_subdirectory(test1) add_subdirectory(test2) +add_subdirectory(test3) +add_subdirectory(test4) diff --git a/tests/HeatTransfer/2D/test3/CMakeLists.txt b/tests/HeatTransfer/2D/test3/CMakeLists.txt new file mode 100644 index 00000000..d5c62fea --- /dev/null +++ b/tests/HeatTransfer/2D/test3/CMakeLists.txt @@ -0,0 +1,6 @@ +create_test("HeatTransfer2Dtest3" "HeatTransfer2Dtest3" FALSE "2D EXE HEAT" 1) +configure_file(${CMAKE_SOURCE_DIR}/tests/tools/convergence_study.py ${CMAKE_CURRENT_BINARY_DIR}/convergence_study.py COPYONLY) +create_col_comparison_convergence("CompareHeatTransfer2Dtest3Convergence" "convergence_output_ref.csv" "convergence_output.csv" -1 absolute 1e-16 FALSE "HeatTransfer2Dtest3" "2D HEAT" 1 100000) + + + diff --git a/tests/HeatTransfer/2D/test3/main.cpp b/tests/HeatTransfer/2D/test3/main.cpp new file mode 100644 index 00000000..3a25129b --- /dev/null +++ b/tests/HeatTransfer/2D/test3/main.cpp @@ -0,0 +1,186 @@ +/** + * @file main.cpp + * @author ci230846 (clement.introini@cea.fr) + * @brief 2D axisymmetric heat transfer problem + * @version 0.1 + * @date 2024-09-3 + * + * @copyright Copyright (c) 2024 + * + */ +#include +#include +#include +#include +#include +#include + +#include "Sloth/sloth.hpp" +#include "Sloth/tests.hpp" +///--------------- +/// Main program +///--------------- +int main(int argc, char* argv[]) { + //--------------------------------------- + // Initialize MPI and HYPRE + //--------------------------------------- + setVerbosity(Verbosity::Debug); + mfem::Mpi::Init(argc, argv); + mfem::Hypre::Init(); + // + //--------------------------------------- + // Profiling + Profiling::getInstance().enable(); + //--------------------------------------- + ///////////////////////// + const int DIM = 2; + using FECollection = Test::FECollection; + using VARS = Test::VARS; + using VAR = Test::VAR; + using PST = Test::PST; + using SPA = Test::SPA; + using BCS = Test::BCS; + // + using OPE = TransientOperator; + using PB = Problem; + + // ########################################### + // ########################################### + // Spatial Discretization // + // ########################################### + // ########################################### + // ############################## + // Meshing // + // ############################## + const std::string& mesh_type = "InlineSquareWithQuadrangles"; // type of mesh + const int order_fe = 1; // finite element order + const int refinement_level = 1; // number of levels of uniform refinement + const double L = 1.0; + std::vector vect_NN{30, 60, 120, 240}; + for (const auto NN : vect_NN) { + const int order = 1; // finite element order + const int n_elements = NN; + + SPA spatial(mesh_type, order_fe, refinement_level, + std::make_tuple(n_elements, n_elements, L, L)); + + // ############################## + // Boundary conditions // + // ############################## + auto boundaries = {Boundary("lower", 0, "Neumann"), Boundary("right", 1, "Neumann"), + Boundary("upper", 2, "Neumann"), Boundary("left", 3, "Neumann")}; + auto Xboundaries = {Boundary("lower", 0, "Neumann"), Boundary("right", 1, "Neumann"), + Boundary("upper", 2, "Neumann"), Boundary("left", 3, "Neumann")}; + auto Tbcs = BCS(&spatial, boundaries); + // ########################################### + // ########################################### + // Physical models // + // ########################################### + // ########################################### + // #################### + // parameters // + // #################### + // Heat + const auto& rho(1.); + const auto& cp(1.); + const auto& cond(1.); + Coefficient density(Glossary::Concentration, rho); + Coefficient heat_capacity(Glossary::Cp, cp); + Coefficient conductivity(Glossary::Conductivity, cond); + + // ############################ + // variables IC + SRC // + // ########################### + // Heat + double qtop = 1.0; + auto T_analytical = std::function( + [qtop]([[maybe_unused]] const mfem::Vector& vcoord, [[maybe_unused]] double time) { + const double r = vcoord[0]; + const double z = vcoord[1]; + + return std::exp(-time) * std::cos(M_PI * r * r) * std::cos(M_PI * z); + }); + + auto heat_vars = + VARS(VAR(&spatial, Tbcs, "T", Glossary::Temperature, 2, + AnalyticalFunctions(T_analytical), AnalyticalFunctions(T_analytical))); + // ########################################### + // ########################################### + // Post-processing // + // ########################################### + // ########################################### + + const std::string& main_folder_path = + "Saves_order_" + std::to_string(order) + "_Nx" + std::to_string(NN); + const auto& level_of_detail = 1; + const auto& frequency = 1; + // Heat + const std::string& calculation_path = "Problem1"; + auto p_pst = Parameters( + Parameter("main_folder_path", main_folder_path), + Parameter("calculation_path", calculation_path), Parameter("frequency", frequency), + Parameter("level_of_detail", level_of_detail), Parameter("enable_compute_energies", false)); + auto pst = PST(&spatial, p_pst); + + // #################### + // Problems // + // #################### + + // Heat: + Coefficients coef_pb(density, heat_capacity, conductivity); //, robin_b); + std::vector spatials{&spatial}; + + auto src_func = std::function( + [rho, cp, cond]([[maybe_unused]] const mfem::Vector& vcoord, [[maybe_unused]] double time) { + const double r = vcoord[0]; + const double z = vcoord[1]; + + const double crr = std::cos(M_PI * r * r); + const double srr = std::sin(M_PI * r * r); + const double cz = std::cos(M_PI * z); + + return r * std::exp(-time) * + (-rho * cp * crr * cz + + cond * (4.0 * M_PI * srr * cz + 4.0 * M_PI * M_PI * r * r * crr * cz + + M_PI * M_PI * crr * cz)); + }); + + std::vector > src_term; + src_term.emplace_back(AnalyticalFunctions(src_func)); + + OPE oper(spatials, {"Fourier"}, TimeScheme::EulerImplicit, "HeatTimeDerivative", src_term); + + oper.overload_nl_solver( + NLSolverType::NEWTON, + Parameters(Parameter("description", "Newton solver "), Parameter("print_level", 1), + Parameter("rel_tol", 1.e-9), Parameter("abs_tol", 1.e-9))); + + PB Heat_pb("Heat", oper, heat_vars, {coef_pb}, pst); + Heat_pb.setGeometry(Geometry::Axisymmetric); + // Coupling 1 + auto cc = Coupling("Heat transfer", Heat_pb); + // ########################################### + // ########################################### + // Time-integration // + // ########################################### + // ########################################### + const auto& t_initial = 0.0; + const auto& t_final = 0.005; + const auto& dt = 0.001; + auto time_params = Parameters(Parameter("initial_time", t_initial), + Parameter("final_time", t_final), Parameter("time_step", dt)); + auto time = TimeDiscretization(time_params, cc); + + time.solve(); + //--------------------------------------- + // Profiling stop + //--------------------------------------- + Profiling::getInstance().print(); + } + //--------------------------------------- + // Finalize MPI + //--------------------------------------- + MPI_Finalize(); + //--------------------------------------- + return 0; +} diff --git a/tests/HeatTransfer/2D/test3/ref/convergence_output_ref.csv b/tests/HeatTransfer/2D/test3/ref/convergence_output_ref.csv new file mode 100644 index 00000000..8616df0d --- /dev/null +++ b/tests/HeatTransfer/2D/test3/ref/convergence_output_ref.csv @@ -0,0 +1,4 @@ +order 1,0.004166666666666667,3.9514e-06 +order 1,0.008333333333333333,1.85752e-05 +order 1,0.016666666666666666,7.73116e-05 +order 1,0.03333333333333333,0.000312207 diff --git a/tests/HeatTransfer/2D/test4/CMakeLists.txt b/tests/HeatTransfer/2D/test4/CMakeLists.txt new file mode 100644 index 00000000..5283bf1b --- /dev/null +++ b/tests/HeatTransfer/2D/test4/CMakeLists.txt @@ -0,0 +1,6 @@ +create_test("HeatTransfer2Dtest4" "HeatTransfer2Dtest4" FALSE "2D EXE HEAT" 1) +configure_file(${CMAKE_SOURCE_DIR}/tests/tools/convergence_study.py ${CMAKE_CURRENT_BINARY_DIR}/convergence_study.py COPYONLY) +create_col_comparison_convergence("CompareHeatTransfer2Dtest4Convergence" "convergence_output_ref.csv" "convergence_output.csv" -1 absolute 1e-16 FALSE "HeatTransfer2Dtest4" "2D HEAT" 1 100000) + + + diff --git a/tests/HeatTransfer/2D/test4/main.cpp b/tests/HeatTransfer/2D/test4/main.cpp new file mode 100644 index 00000000..eba3bd9a --- /dev/null +++ b/tests/HeatTransfer/2D/test4/main.cpp @@ -0,0 +1,191 @@ +/** + * @file main.cpp + * @author ci230846 (clement.introini@cea.fr) + * @brief 2D axisymmetric heat transfer problem + * @version 0.1 + * @date 2024-09-3 + * + * @copyright Copyright (c) 2024 + * + */ +#include +#include +#include +#include +#include +#include + +#include "Sloth/sloth.hpp" +#include "Sloth/tests.hpp" +///--------------- +/// Main program +///--------------- +int main(int argc, char* argv[]) { + //--------------------------------------- + // Initialize MPI and HYPRE + //--------------------------------------- + setVerbosity(Verbosity::Debug); + mfem::Mpi::Init(argc, argv); + mfem::Hypre::Init(); + // + //--------------------------------------- + // Profiling + Profiling::getInstance().enable(); + //--------------------------------------- + ///////////////////////// + const int DIM = 2; + using FECollection = Test::FECollection; + using VARS = Test::VARS; + using VAR = Test::VAR; + using PST = Test::PST; + using SPA = Test::SPA; + using BCS = Test::BCS; + // + using OPE = TransientOperator; + using PB = Problem; + + // ########################################### + // ########################################### + // Spatial Discretization // + // ########################################### + // ########################################### + // ############################## + // Meshing // + // ############################## + const std::string& mesh_type = "InlineSquareWithQuadrangles"; // type of mesh + const int order_fe = 1; // finite element order + const int refinement_level = 1; // number of levels of uniform refinement + const double L = 1.0; + std::vector vect_NN{30, 60, 120, 240}; + for (const auto NN : vect_NN) { + const int order = 1; // finite element order + const int n_elements = NN; + + SPA spatial(mesh_type, order_fe, refinement_level, + std::make_tuple(n_elements, n_elements, L, L)); + + // ############################## + // Boundary conditions // + // ############################## + auto boundaries = {Boundary("lower", 0, "Neumann"), Boundary("right", 1, "Neumann"), + Boundary("upper", 2, "Neumann"), Boundary("left", 3, "Neumann")}; + auto Xboundaries = {Boundary("lower", 0, "Neumann"), Boundary("right", 1, "Neumann"), + Boundary("upper", 2, "Neumann"), Boundary("left", 3, "Neumann")}; + auto Tbcs = BCS(&spatial, boundaries); + // ########################################### + // ########################################### + // Physical models // + // ########################################### + // ########################################### + // #################### + // parameters // + // #################### + // Heat + const auto& rho(1.); + const auto& cp(1.); + const auto& cond(1.); + Coefficient density(Glossary::Concentration, rho); + Coefficient heat_capacity(Glossary::Cp, cp); + Coefficient conductivity(Glossary::Conductivity, cond); + + // ############################ + // variables IC + SRC // + // ########################### + // Heat + double qtop = 1.0; + auto T_analytical = std::function( + [qtop]([[maybe_unused]] const mfem::Vector& vcoord, [[maybe_unused]] double time) { + const double r = vcoord[0]; + const double z = vcoord[1]; + + return std::exp(-time) * std::cos(M_PI * r * r) * std::cos(M_PI * z) - qtop * z; + }); + + Coefficient robin_b(Glossary::Neumann, -1.0); + Coefficient robin_bb(Glossary::Neumann, 1.0); + robin_b.set_bdr_index_coef(std::vector{2}); + robin_bb.set_bdr_index_coef(std::vector{0}); + + auto heat_vars = + VARS(VAR(&spatial, Tbcs, "T", Glossary::Temperature, 2, + AnalyticalFunctions(T_analytical), AnalyticalFunctions(T_analytical))); + // ########################################### + // ########################################### + // Post-processing // + // ########################################### + // ########################################### + + const std::string& main_folder_path = + "Saves_order_" + std::to_string(order) + "_Nx" + std::to_string(NN); + const auto& level_of_detail = 1; + const auto& frequency = 1; + // Heat + const std::string& calculation_path = "Problem1"; + auto p_pst = Parameters( + Parameter("main_folder_path", main_folder_path), + Parameter("calculation_path", calculation_path), Parameter("frequency", frequency), + Parameter("level_of_detail", level_of_detail), Parameter("enable_compute_energies", false)); + auto pst = PST(&spatial, p_pst); + + // #################### + // Problems // + // #################### + + // Heat: + Coefficients coef_pb(density, heat_capacity, conductivity, robin_b, robin_bb); + std::vector spatials{&spatial}; + + auto src_func = std::function( + [rho, cp, cond]([[maybe_unused]] const mfem::Vector& vcoord, [[maybe_unused]] double time) { + const double r = vcoord[0]; + const double z = vcoord[1]; + + const double crr = std::cos(M_PI * r * r); + const double srr = std::sin(M_PI * r * r); + const double cz = std::cos(M_PI * z); + + return r * std::exp(-time) * + (-rho * cp * crr * cz + + cond * (4.0 * M_PI * srr * cz + 4.0 * M_PI * M_PI * r * r * crr * cz + + M_PI * M_PI * crr * cz)); + }); + + std::vector > src_term; + src_term.emplace_back(AnalyticalFunctions(src_func)); + + OPE oper(spatials, {"Fourier"}, TimeScheme::EulerImplicit, "HeatTimeDerivative", src_term); + + oper.overload_nl_solver( + NLSolverType::NEWTON, + Parameters(Parameter("description", "Newton solver "), Parameter("print_level", 1), + Parameter("rel_tol", 1.e-9), Parameter("abs_tol", 1.e-9))); + + PB Heat_pb("Heat", oper, heat_vars, {coef_pb}, pst); + Heat_pb.setGeometry(Geometry::Axisymmetric); + // Coupling 1 + auto cc = Coupling("Heat transfer", Heat_pb); + // ########################################### + // ########################################### + // Time-integration // + // ########################################### + // ########################################### + const auto& t_initial = 0.0; + const auto& t_final = 0.01; + const auto& dt = 0.001; + auto time_params = Parameters(Parameter("initial_time", t_initial), + Parameter("final_time", t_final), Parameter("time_step", dt)); + auto time = TimeDiscretization(time_params, cc); + + time.solve(); + //--------------------------------------- + // Profiling stop + //--------------------------------------- + Profiling::getInstance().print(); + } + //--------------------------------------- + // Finalize MPI + //--------------------------------------- + MPI_Finalize(); + //--------------------------------------- + return 0; +} diff --git a/tests/HeatTransfer/2D/test4/ref/convergence_output_ref.csv b/tests/HeatTransfer/2D/test4/ref/convergence_output_ref.csv new file mode 100644 index 00000000..f1cc90b6 --- /dev/null +++ b/tests/HeatTransfer/2D/test4/ref/convergence_output_ref.csv @@ -0,0 +1,4 @@ +order 1,0.004166666666666667,3.2037e-06 +order 1,0.008333333333333333,1.68338e-05 +order 1,0.016666666666666666,7.26958e-05 +order 1,0.03333333333333333,0.000296269