From 684954d906126ab55df5a7c760d3c55d1a314a4d Mon Sep 17 00:00:00 2001 From: Edward Palmer Date: Fri, 28 Jun 2024 10:42:10 +0000 Subject: [PATCH] Removes EM formulations. --- include/formulations/CustomFormulation.h | 9 +- include/formulations/MFEMFormulation.h | 4 +- include/formulations/a_formulation.h | 55 ---- include/formulations/av_formulation.h | 53 ---- include/formulations/complex_a_formulation.h | 88 ------- include/formulations/complex_e_formulation.h | 90 ------- .../complex_em_formulation_interface.h | 66 ----- .../complex_maxwell_formulation.h | 92 ------- include/formulations/dual_formulation.h | 77 ------ include/formulations/e_formulation.h | 38 --- include/formulations/eb_dual_formulation.h | 44 ---- .../formulations/em_formulation_interface.h | 63 ----- include/formulations/formulation.h | 4 - .../frequency_domain_em_formulation.h | 19 -- include/formulations/h_formulation.h | 54 ---- include/formulations/hcurl_formulation.h | 49 ---- .../formulations/magnetostatic_formulation.h | 37 --- include/formulations/statics_formulation.h | 51 ---- .../steady_state_em_formulation.h | 16 -- .../formulations/time_domain_em_formulation.h | 16 -- src/formulations/CustomFormulation.C | 6 +- src/formulations/MFEMFormulation.C | 2 - src/formulations/a_formulation.C | 145 ----------- src/formulations/av_formulation.C | 194 -------------- src/formulations/complex_a_formulation.C | 137 ---------- src/formulations/complex_e_formulation.C | 131 ---------- .../complex_maxwell_formulation.C | 224 ---------------- src/formulations/dual_formulation.C | 241 ------------------ src/formulations/e_formulation.C | 91 ------- src/formulations/eb_dual_formulation.C | 91 ------- .../frequency_domain_em_formulation.C | 8 - src/formulations/h_formulation.C | 131 ---------- src/formulations/hcurl_formulation.C | 167 ------------ src/formulations/magnetostatic_formulation.C | 84 ------ src/formulations/statics_formulation.C | 159 ------------ .../steady_state_em_formulation.C | 8 - src/formulations/time_domain_em_formulation.C | 58 ----- 37 files changed, 9 insertions(+), 2793 deletions(-) delete mode 100644 include/formulations/a_formulation.h delete mode 100644 include/formulations/av_formulation.h delete mode 100644 include/formulations/complex_a_formulation.h delete mode 100644 include/formulations/complex_e_formulation.h delete mode 100644 include/formulations/complex_em_formulation_interface.h delete mode 100644 include/formulations/complex_maxwell_formulation.h delete mode 100644 include/formulations/dual_formulation.h delete mode 100644 include/formulations/e_formulation.h delete mode 100644 include/formulations/eb_dual_formulation.h delete mode 100644 include/formulations/em_formulation_interface.h delete mode 100644 include/formulations/formulation.h delete mode 100644 include/formulations/frequency_domain_em_formulation.h delete mode 100644 include/formulations/h_formulation.h delete mode 100644 include/formulations/hcurl_formulation.h delete mode 100644 include/formulations/magnetostatic_formulation.h delete mode 100644 include/formulations/statics_formulation.h delete mode 100644 include/formulations/steady_state_em_formulation.h delete mode 100644 include/formulations/time_domain_em_formulation.h delete mode 100644 src/formulations/a_formulation.C delete mode 100644 src/formulations/av_formulation.C delete mode 100644 src/formulations/complex_a_formulation.C delete mode 100644 src/formulations/complex_e_formulation.C delete mode 100644 src/formulations/complex_maxwell_formulation.C delete mode 100644 src/formulations/dual_formulation.C delete mode 100644 src/formulations/e_formulation.C delete mode 100644 src/formulations/eb_dual_formulation.C delete mode 100644 src/formulations/frequency_domain_em_formulation.C delete mode 100644 src/formulations/h_formulation.C delete mode 100644 src/formulations/hcurl_formulation.C delete mode 100644 src/formulations/magnetostatic_formulation.C delete mode 100644 src/formulations/statics_formulation.C delete mode 100644 src/formulations/steady_state_em_formulation.C delete mode 100644 src/formulations/time_domain_em_formulation.C diff --git a/include/formulations/CustomFormulation.h b/include/formulations/CustomFormulation.h index 6a8ce609..6d061e9d 100644 --- a/include/formulations/CustomFormulation.h +++ b/include/formulations/CustomFormulation.h @@ -1,7 +1,7 @@ #pragma once #include "MFEMFormulation.h" -#include "time_domain_em_formulation.h" +#include "time_domain_equation_system_problem_builder.h" class CustomFormulation : public MFEMFormulation { @@ -9,14 +9,15 @@ class CustomFormulation : public MFEMFormulation static InputParameters validParams(); CustomFormulation(const InputParameters & parameters); - virtual ~CustomFormulation(); + ~CustomFormulation() override = default; virtual void execute() override {} virtual void initialize() override {} virtual void finalize() override {} - std::shared_ptr getProblemBuilder() override { return formulation; } + /// Returns a shared pointer to the time-domain equation system problem builder. + std::shared_ptr getProblemBuilder() const override { return _formulation; } private: - std::shared_ptr formulation{nullptr}; + const std::shared_ptr _formulation{nullptr}; }; diff --git a/include/formulations/MFEMFormulation.h b/include/formulations/MFEMFormulation.h index d7e37c17..b2f194f2 100644 --- a/include/formulations/MFEMFormulation.h +++ b/include/formulations/MFEMFormulation.h @@ -9,13 +9,13 @@ class MFEMFormulation : public GeneralUserObject static InputParameters validParams(); MFEMFormulation(const InputParameters & parameters); - virtual ~MFEMFormulation(); + virtual ~MFEMFormulation() override = default; virtual void execute() override {} virtual void initialize() override {} virtual void finalize() override {} - virtual std::shared_ptr getProblemBuilder() + virtual std::shared_ptr getProblemBuilder() const { mooseError( "Base class MFEMFormulation cannot return a valid ProblemBuilder. Use a child class."); diff --git a/include/formulations/a_formulation.h b/include/formulations/a_formulation.h deleted file mode 100644 index 4fbf7628..00000000 --- a/include/formulations/a_formulation.h +++ /dev/null @@ -1,55 +0,0 @@ -#pragma once -#include "../common/pfem_extras.hpp" -#include "hcurl_formulation.h" -#include "inputs.h" - -namespace platypus -{ - -class AFormulation : public platypus::HCurlFormulation -{ -public: - AFormulation(const std::string & magnetic_reluctivity_name, - std::string magnetic_permeability_name, - const std::string & electric_conductivity_name, - const std::string & magnetic_vector_potential_name); - - ~AFormulation() override = default; - - // Enable auxiliary calculation of J ∈ H(div) - void RegisterCurrentDensityAux(const std::string & j_field_name, - const std::string & external_j_field_name = "") override; - - // Enable auxiliary calculation of B ∈ H(div) - void RegisterMagneticFluxDensityAux(const std::string & b_field_name, - const std::string & external_b_field_name = "") override; - - // Enable auxiliary calculation of E ∈ H(curl) - void RegisterElectricFieldAux(const std::string & e_field_name, - const std::string & external_e_field_name = "") override; - - // Enable auxiliary calculation of H ∈ H(curl) - void RegisterMagneticFieldAux(const std::string & h_field_name, - const std::string & external_h_field_name = "") override; - - // Enable auxiliary calculation of F ∈ L2 - void RegisterLorentzForceDensityAux(const std::string & f_field_name, - const std::string & b_field_name, - const std::string & j_field_name) override; - - // Enable auxiliary calculation of P ∈ L2 - void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name, - const std::string & j_field_name) override; - - void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name) override; - - void RegisterCoefficients() override; - -protected: - const std::string _magnetic_permeability_name; - const std::string & _magnetic_reluctivity_name = platypus::HCurlFormulation::_alpha_coef_name; - const std::string & _electric_conductivity_name = platypus::HCurlFormulation::_beta_coef_name; -}; -} // namespace platypus diff --git a/include/formulations/av_formulation.h b/include/formulations/av_formulation.h deleted file mode 100644 index eb34f1e7..00000000 --- a/include/formulations/av_formulation.h +++ /dev/null @@ -1,53 +0,0 @@ -#pragma once -#include "../common/pfem_extras.hpp" -#include "formulation.h" -#include "inputs.h" -#include "sources.h" - -namespace platypus -{ - -class AVFormulation : public TimeDomainEMFormulation -{ -public: - AVFormulation(std::string alpha_coef_name, - std::string inv_alpha_coef_name, - std::string beta_coef_name, - std::string vector_potential_name, - std::string scalar_potential_name); - - ~AVFormulation() override = default; - - void ConstructOperator() override; - - void RegisterGridFunctions() override; - - void RegisterCoefficients() override; - -protected: - const std::string _alpha_coef_name; - const std::string _inv_alpha_coef_name; - const std::string _beta_coef_name; - const std::string _vector_potential_name; - const std::string _scalar_potential_name; -}; - -class AVEquationSystem : public TimeDependentEquationSystem -{ -public: - AVEquationSystem(const platypus::InputParameters & params); - - ~AVEquationSystem() override = default; - - void Init(platypus::GridFunctions & gridfunctions, - const platypus::FESpaces & fespaces, - platypus::BCMap & bc_map, - platypus::Coefficients & coefficients) override; - void AddKernels() override; - - std::string _a_name, _v_name, _coupled_variable_name, _alpha_coef_name, _beta_coef_name, - _dtalpha_coef_name, _neg_beta_coef_name; - mfem::ConstantCoefficient _neg_coef; -}; - -} // namespace platypus diff --git a/include/formulations/complex_a_formulation.h b/include/formulations/complex_a_formulation.h deleted file mode 100644 index 4fbc855f..00000000 --- a/include/formulations/complex_a_formulation.h +++ /dev/null @@ -1,88 +0,0 @@ -#pragma once -#include "complex_maxwell_formulation.h" - -namespace platypus -{ -/* -Formulation for solving: -∇×(ν∇×A) + iωσA - ω²εA = J - -via the weak form: -(ν∇×A, ∇×A') + (iωσA, A') - (ω²εA, A') - <(ν∇×A)×n, A'> = (J, A') - -where -A ∈ H(curl) is the magnetic vector potential -J ∈ H(div) is the divergence-free source current density -ω is the angular frequency -ν is the reluctivity (reciprocal of magnetic permeability; 1/µ) -σ is the electric conductivity -ε is the electric permittivity - -Dirichlet boundaries strongly constrain n×n×A -Integrated boundaries weakly constrain (ν∇×A)×n = H×n -Robin boundaries weakly constrain (ν∇×A)×n + γ(n×n×A) = F - -Divergence cleaning (such as via Helmholtz projection) -should be performed on J before use in this operator. -*/ -class ComplexAFormulation : public platypus::ComplexMaxwellFormulation -{ - -public: - ComplexAFormulation(const std::string & magnetic_reluctivity_name, - const std::string & electric_conductivity_name, - const std::string & dielectric_permittivity_name, - const std::string & frequency_coef_name, - const std::string & magnetic_vector_potential_complex_name, - const std::string & magnetic_vector_potential_real_name, - const std::string & magnetic_vector_potential_imag_name); - - ~ComplexAFormulation() override = default; - - // Enable auxiliary calculation of J ∈ H(div) - //* Induced electric current, Jind = σE = -iωσA - void RegisterCurrentDensityAux(const std::string & j_field_real_name, - const std::string & j_field_imag_name, - const std::string & external_j_field_real_name = "", - const std::string & external_j_field_imag_name = "") override; - - //* Magnetic flux density B = curl A - void RegisterMagneticFluxDensityAux(const std::string & b_field_real_name, - const std::string & b_field_imag_name, - const std::string & external_b_field_real_name = "", - const std::string & external_b_field_imag_name = "") override; - - //* Electric field E =-dA/dt=-iωA - void RegisterElectricFieldAux(const std::string & e_field_real_name, - const std::string & e_field_imag_name, - const std::string & external_e_field_real_name = "", - const std::string & external_e_field_imag_name = "") override; - - // Enable auxiliary calculation of P ∈ L2 - // Time averaged Joule heating density E.J - void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_real_name, - const std::string & e_field_imag_name, - const std::string & j_field_real_name, - const std::string & j_field_imag_name) override; - - // Time averaged Joule heating density σ|E|^2 - void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_real_name, - const std::string & e_field_imag_name) override; - -protected: - const std::string & _magnetic_reluctivity_name = - platypus::ComplexMaxwellFormulation::_alpha_coef_name; - const std::string & _electric_conductivity_name = - platypus::ComplexMaxwellFormulation::_beta_coef_name; - const std::string & _electric_permittivity_name = - platypus::ComplexMaxwellFormulation::_zeta_coef_name; - const std::string & _magnetic_vector_potential_complex_name = - platypus::ComplexMaxwellFormulation::_h_curl_var_complex_name; - const std::string & _magnetic_vector_potential_real_name = - platypus::ComplexMaxwellFormulation::_h_curl_var_real_name; - const std::string & _magnetic_vector_potential_imag_name = - platypus::ComplexMaxwellFormulation::_h_curl_var_imag_name; -}; -} // namespace platypus diff --git a/include/formulations/complex_e_formulation.h b/include/formulations/complex_e_formulation.h deleted file mode 100644 index 3d21426a..00000000 --- a/include/formulations/complex_e_formulation.h +++ /dev/null @@ -1,90 +0,0 @@ -#pragma once -#include "complex_maxwell_formulation.h" - -namespace platypus -{ -/* -Formulation for solving: -∇×(ν∇×E꜀) + iωσE꜀ - ω²εE꜀ = -iωJ꜀ᵉ - -via the weak form: -(ν∇×E꜀, ∇×E꜀') + (iωσE꜀, E꜀') - (ω²εE꜀, E꜀') - <(ν∇×E꜀)×n, E꜀'> = -(iωJ꜀ᵉ, E꜀') - -where -E꜀ ∈ H(curl) is the complex electric field -J꜀ᵉ ∈ H(div) is the divergence-free source current density -ω is the angular frequency -ν is the reluctivity (reciprocal of magnetic permeability; 1/µ) -σ is the electric conductivity -ε is the electric permittivity - -Dirichlet boundaries strongly constrain n×n×E꜀ -Integrated boundaries weakly constrain (ν∇×E꜀)×n = n×dB꜀/dt -Robin boundaries weakly constrain (ν∇×E꜀)×n + γ(n×n×E꜀) = F - -Divergence cleaning (such as via Helmholtz projection) -should be performed on J꜀ᵉ before use in this operator. -*/ -class ComplexEFormulation : public platypus::ComplexMaxwellFormulation -{ -public: - ComplexEFormulation(const std::string & magnetic_reluctivity_name, - const std::string & electric_conductivity_name, - const std::string & dielectric_permittivity_name, - const std::string & frequency_coef_name, - const std::string & e_field_complex_name, - const std::string & e_field_real_name, - const std::string & e_field_imag_name); - - ~ComplexEFormulation() override = default; - - // Enable auxiliary calculation of J ∈ H(div) - //* Induced electric current, Jind = σE - void RegisterCurrentDensityAux(const std::string & j_field_real_name, - const std::string & j_field_imag_name, - const std::string & external_j_field_real_name = "", - const std::string & external_j_field_imag_name = "") override; - - //* Magnetic flux density B = (i/ω) curl E - //* (∇×E = -dB/dt = -iωB) - void RegisterMagneticFluxDensityAux(const std::string & b_field_real_name, - const std::string & b_field_imag_name, - const std::string & external_b_field_real_name = "", - const std::string & external_b_field_imag_name = "") override; - - //* Electric field is already a state variable solved for - void RegisterElectricFieldAux(const std::string & e_field_real_name, - const std::string & e_field_imag_name, - const std::string & external_e_field_real_name = "", - const std::string & external_e_field_imag_name = "") override - { - } - - // Enable auxiliary calculation of P ∈ L2 - // Time averaged Joule heating density E.J - void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_real_name, - const std::string & e_field_imag_name, - const std::string & j_field_real_name, - const std::string & j_field_imag_name) override; - - // Time averaged Joule heating density σ|E|^2 - void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_real_name, - const std::string & e_field_imag_name) override; - -protected: - const std::string & _magnetic_reluctivity_name = - platypus::ComplexMaxwellFormulation::_alpha_coef_name; - const std::string & _electric_conductivity_name = - platypus::ComplexMaxwellFormulation::_beta_coef_name; - const std::string & _electric_permittivity_name = - platypus::ComplexMaxwellFormulation::_zeta_coef_name; - const std::string & _electric_field_complex_name = - platypus::ComplexMaxwellFormulation::_h_curl_var_complex_name; - const std::string & _electric_field_real_name = - platypus::ComplexMaxwellFormulation::_h_curl_var_real_name; - const std::string & _electric_field_imag_name = - platypus::ComplexMaxwellFormulation::_h_curl_var_imag_name; -}; -} // namespace platypus diff --git a/include/formulations/complex_em_formulation_interface.h b/include/formulations/complex_em_formulation_interface.h deleted file mode 100644 index f998342b..00000000 --- a/include/formulations/complex_em_formulation_interface.h +++ /dev/null @@ -1,66 +0,0 @@ -#pragma once -#include "mfem.hpp" - -namespace platypus -{ - -// Specifies interface specific to EM formulations. -class ComplexEMFormulationInterface -{ -public: - ComplexEMFormulationInterface() = default; - - // Enable auxiliary calculation of J ∈ H(div) - virtual void RegisterCurrentDensityAux(const std::string & j_field_real_name, - const std::string & j_field_imag_name, - const std::string & external_j_field_real_name = "", - const std::string & external_j_field_imag_name = "") - { - MFEM_ABORT("Current density auxsolver not available for this formulation"); - } - - // Enable auxiliary calculation of B ∈ H(div) - virtual void RegisterMagneticFluxDensityAux(const std::string & b_field_real_name, - const std::string & b_field_imag_name, - const std::string & external_b_field_real_name = "", - const std::string & external_b_field_imag_name = "") - { - MFEM_ABORT("Magnetic flux density auxsolver not available for this formulation"); - } - - // Enable auxiliary calculation of E ∈ H(curl) - virtual void RegisterElectricFieldAux(const std::string & e_field_real_name, - const std::string & e_field_imag_name, - const std::string & external_e_field_real_name = "", - const std::string & external_e_field_imag_name = "") - { - MFEM_ABORT("Electric field auxsolver not available for this formulation"); - } - - // Enable auxiliary calculation of H ∈ H(curl) - virtual void RegisterMagneticFieldAux(const std::string & h_field_real_name, - const std::string & h_field_imag_name, - const std::string & external_h_field_real_name = "", - const std::string & external_h_field_imag_name = "") - { - MFEM_ABORT("Magnetic field auxsolver not available for this formulation"); - } - - // Enable auxiliary calculation of P ∈ L2 - virtual void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_real_name, - const std::string & e_field_imag_name, - const std::string & j_field_real_name, - const std::string & j_field_imag_name) - { - MFEM_ABORT("Joule heating auxsolver not available for this formulation"); - } - - virtual void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_real_name, - const std::string & e_field_imag_name) - { - MFEM_ABORT("Joule heating auxsolver not available for this formulation"); - } -}; -} // namespace platypus diff --git a/include/formulations/complex_maxwell_formulation.h b/include/formulations/complex_maxwell_formulation.h deleted file mode 100644 index 85a1790f..00000000 --- a/include/formulations/complex_maxwell_formulation.h +++ /dev/null @@ -1,92 +0,0 @@ -#pragma once -#include "frequency_domain_em_formulation.h" - -namespace platypus -{ -/* -Formulation for solving: -∇×(α∇×u) + iωβu - ω²γu = g - -via the weak form: -(α∇×u, ∇×u') + (iωβu, u') - (ω²ζu, u') - <(α∇×u)×n, u'> = (g, u') - -where -u ∈ H(curl) -g ∈ H(div) is a divergence-free source field -ω is the angular frequency -α is the stiffness coefficient -ωβ is the loss coefficient -ω²γ is the mass coefficient - -Dirichlet boundaries strongly constrain n×n×u -Integrated boundaries weakly constrain (α∇×u)×n -Robin boundaries weakly constrain (α∇×u)×n + γ(n×n×u) = F - -Divergence cleaning (such as via Helmholtz projection) -should be performed on g before use in this operator. -*/ -class ComplexMaxwellFormulation : public platypus::FrequencyDomainEMFormulation -{ -public: - ComplexMaxwellFormulation(std::string frequency_coef_name, - std::string alpha_coef_name, - std::string beta_coef_name, - std::string zeta_coef_name, - std::string h_curl_var_complex_name, - std::string h_curl_var_real_name, - std::string h_curl_var_imag_name); - - ~ComplexMaxwellFormulation() override = default; - - void ConstructJacobianSolver() override; - - void ConstructOperator() override; - - void RegisterGridFunctions() override; - - void RegisterCoefficients() override; - - // std::vector local_trial_vars, local_test_vars; -protected: - const std::string _alpha_coef_name; - const std::string _beta_coef_name; - const std::string _zeta_coef_name; - const std::string _frequency_coef_name; - const std::string _h_curl_var_complex_name; - const std::string _h_curl_var_real_name; - const std::string _h_curl_var_imag_name; - const std::string _mass_coef_name; - const std::string _loss_coef_name; -}; - -class ComplexMaxwellOperator : public ProblemOperator -{ -public: - ComplexMaxwellOperator(platypus::Problem & problem, - std::string h_curl_var_complex_name, - std::string h_curl_var_real_name, - std::string h_curl_var_imag_name, - std::string stiffness_coef_name, - std::string mass_coef_name, - std::string loss_coef_name); - - ~ComplexMaxwellOperator() override = default; - - void SetGridFunctions() override; - void Init(mfem::Vector & X) override; - void Solve(mfem::Vector & X) override; - - std::string _h_curl_var_complex_name, _h_curl_var_real_name, _h_curl_var_imag_name, - _stiffness_coef_name, _mass_coef_name, _loss_coef_name; - - mfem::ComplexOperator::Convention _conv{mfem::ComplexOperator::HERMITIAN}; - - std::unique_ptr _u{nullptr}; - mfem::Coefficient * _stiff_coef{nullptr}; // Dia/Paramagnetic Material Coefficient - mfem::Coefficient * _mass_coef{nullptr}; // -omega^2 epsilon - mfem::Coefficient * _loss_coef{nullptr}; // omega sigma - - mfem::Array _ess_bdr_tdofs; -}; - -} // namespace platypus diff --git a/include/formulations/dual_formulation.h b/include/formulations/dual_formulation.h deleted file mode 100644 index 185f8d85..00000000 --- a/include/formulations/dual_formulation.h +++ /dev/null @@ -1,77 +0,0 @@ -#pragma once -#include "../common/pfem_extras.hpp" -#include "formulation.h" -#include "inputs.h" - -namespace platypus -{ - -class DualFormulation : public TimeDomainEMFormulation -{ -public: - DualFormulation(std::string alpha_coef_name, - std::string beta_coef_name, - std::string h_curl_var_name, - std::string h_div_var_name); - - ~DualFormulation() override = default; - - void ConstructJacobianPreconditioner() override; - - void ConstructJacobianSolver() override; - - void ConstructOperator() override; - - void RegisterGridFunctions() override; - - void RegisterCoefficients() override; - -protected: - const std::string _alpha_coef_name; - const std::string _beta_coef_name; - const std::string _h_curl_var_name; - const std::string _h_div_var_name; -}; - -class WeakCurlEquationSystem : public TimeDependentEquationSystem -{ -public: - WeakCurlEquationSystem(const platypus::InputParameters & params); - ~WeakCurlEquationSystem() override = default; - - void Init(platypus::GridFunctions & gridfunctions, - const platypus::FESpaces & fespaces, - platypus::BCMap & bc_map, - platypus::Coefficients & coefficients) override; - void AddKernels() override; - - std::string _h_curl_var_name, _h_div_var_name, _alpha_coef_name, _beta_coef_name, - _dtalpha_coef_name; -}; - -class DualOperator : public TimeDomainEquationSystemProblemOperator -{ -public: - DualOperator(platypus::Problem & problem, - std::unique_ptr equation_system) - : TimeDomainEquationSystemProblemOperator(problem, std::move(equation_system)) - { - } - - void Init(mfem::Vector & X) override; - - void ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vector & dX_dt) override; - void SetGridFunctions() override; - - mfem::ParFiniteElementSpace * _h_curl_fe_space{nullptr}; - mfem::ParFiniteElementSpace * _h_div_fe_space{nullptr}; - - std::string _h_curl_var_name, _h_div_var_name; - - mfem::ParGridFunction * _u{nullptr}; // HCurl vector field - mfem::ParGridFunction * _dv{nullptr}; // HDiv vector field - -protected: - std::unique_ptr _curl; -}; -} // namespace platypus diff --git a/include/formulations/e_formulation.h b/include/formulations/e_formulation.h deleted file mode 100644 index 5942f98d..00000000 --- a/include/formulations/e_formulation.h +++ /dev/null @@ -1,38 +0,0 @@ -#pragma once -#include "../common/pfem_extras.hpp" -#include "hcurl_formulation.h" -#include "inputs.h" - -namespace platypus -{ - -class EFormulation : public platypus::HCurlFormulation -{ -public: - EFormulation(const std::string & magnetic_reluctivity_name, - std::string magnetic_permeability_name, - const std::string & electric_conductivity_name, - const std::string & e_field_name); - - ~EFormulation() override = default; - - void RegisterCoefficients() override; - - // Enable auxiliary calculation of J ∈ H(div) - void RegisterCurrentDensityAux(const std::string & j_field_name, - const std::string & external_j_field_name = "") override; - - // Enable auxiliary calculation of P ∈ L2 - void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name, - const std::string & j_field_name) override; - void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name) override; - -protected: - const std::string _magnetic_permeability_name; - const std::string & _magnetic_reluctivity_name = platypus::HCurlFormulation::_alpha_coef_name; - const std::string & _electric_conductivity_name = platypus::HCurlFormulation::_beta_coef_name; -}; - -} // namespace platypus diff --git a/include/formulations/eb_dual_formulation.h b/include/formulations/eb_dual_formulation.h deleted file mode 100644 index 8cd21003..00000000 --- a/include/formulations/eb_dual_formulation.h +++ /dev/null @@ -1,44 +0,0 @@ -#pragma once -#include "../common/pfem_extras.hpp" -#include "dual_formulation.h" -#include "inputs.h" - -namespace platypus -{ - -class EBDualFormulation : public platypus::DualFormulation -{ -public: - EBDualFormulation(const std::string & magnetic_reluctivity_name, - std::string magnetic_permeability_name, - const std::string & electric_conductivity_name, - const std::string & e_field_name, - const std::string & b_field_name); - - ~EBDualFormulation() override = default; - - // Enable auxiliary calculation of J ∈ H(div) - void RegisterCurrentDensityAux(const std::string & j_field_name, - const std::string & external_j_field_name = "") override; - - // Enable auxiliary calculation of F ∈ L2 - void RegisterLorentzForceDensityAux(const std::string & f_field_name, - const std::string & b_field_name, - const std::string & j_field_name) override; - - // Enable auxiliary calculation of P ∈ L2 - void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name, - const std::string & j_field_name) override; - void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name) override; - - void RegisterCoefficients() override; - -protected: - const std::string _magnetic_permeability_name; - const std::string & _magnetic_reluctivity_name = platypus::DualFormulation::_alpha_coef_name; - const std::string & _electric_conductivity_name = platypus::DualFormulation::_beta_coef_name; -}; - -} // namespace platypus diff --git a/include/formulations/em_formulation_interface.h b/include/formulations/em_formulation_interface.h deleted file mode 100644 index 03c29754..00000000 --- a/include/formulations/em_formulation_interface.h +++ /dev/null @@ -1,63 +0,0 @@ -#pragma once -#include "mfem.hpp" - -namespace platypus -{ - -// Specifies interface specific to EM formulations. -class EMFormulationInterface -{ -public: - EMFormulationInterface() = default; - - // Enable auxiliary calculation of J ∈ H(div) - virtual void RegisterCurrentDensityAux(const std::string & j_field_name, - const std::string & external_j_field_name = "") - { - MFEM_ABORT("Current density auxsolver not available for this formulation"); - } - - // Enable auxiliary calculation of B ∈ H(div) - virtual void RegisterMagneticFluxDensityAux(const std::string & b_field_name, - const std::string & external_b_field_name = "") - { - MFEM_ABORT("Magnetic flux density auxsolver not available for this formulation"); - } - - // Enable auxiliary calculation of E ∈ H(curl) - virtual void RegisterElectricFieldAux(const std::string & e_field_name, - const std::string & external_e_field_name = "") - { - MFEM_ABORT("Electric field auxsolver not available for this formulation"); - } - - // Enable auxiliary calculation of H ∈ H(curl) - virtual void RegisterMagneticFieldAux(const std::string & h_field_name, - const std::string & external_h_field_name = "") - { - MFEM_ABORT("Magnetic field auxsolver not available for this formulation"); - } - - // Enable auxiliary calculation of F ∈ L2 - virtual void RegisterLorentzForceDensityAux(const std::string & f_field_name, - const std::string & b_field_name, - const std::string & j_field_name) - { - MFEM_ABORT("Lorentz force auxsolver not available for this formulation"); - } - - // Enable auxiliary calculation of P ∈ L2 - virtual void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name, - const std::string & j_field_name) - { - MFEM_ABORT("Joule heating auxsolver not available for this formulation"); - } - - virtual void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name) - { - MFEM_ABORT("Joule heating auxsolver not available for this formulation"); - } -}; -} // namespace platypus diff --git a/include/formulations/formulation.h b/include/formulations/formulation.h deleted file mode 100644 index 940f72da..00000000 --- a/include/formulations/formulation.h +++ /dev/null @@ -1,4 +0,0 @@ -#pragma once -#include "frequency_domain_em_formulation.h" -#include "steady_state_em_formulation.h" -#include "time_domain_em_formulation.h" diff --git a/include/formulations/frequency_domain_em_formulation.h b/include/formulations/frequency_domain_em_formulation.h deleted file mode 100644 index 686cb1a9..00000000 --- a/include/formulations/frequency_domain_em_formulation.h +++ /dev/null @@ -1,19 +0,0 @@ -#pragma once -#include "complex_em_formulation_interface.h" -#include "steady_state_problem_builder.h" - -namespace platypus -{ - -// Specifies output interfaces of a frequency-domain EM formulation. -class FrequencyDomainEMFormulation : public platypus::SteadyStateProblemBuilder, - public platypus::ComplexEMFormulationInterface -{ -public: - FrequencyDomainEMFormulation(); - ~FrequencyDomainEMFormulation() override = default; - -protected: - mfem::ConstantCoefficient * _freq_coef{nullptr}; -}; -} // namespace platypus diff --git a/include/formulations/h_formulation.h b/include/formulations/h_formulation.h deleted file mode 100644 index 89bbc186..00000000 --- a/include/formulations/h_formulation.h +++ /dev/null @@ -1,54 +0,0 @@ -#pragma once -#include "../common/pfem_extras.hpp" -#include "hcurl_formulation.h" -#include "inputs.h" - -namespace platypus -{ - -class HFormulation : public platypus::HCurlFormulation -{ -public: - HFormulation(const std::string & electric_resistivity_name, - std::string electric_conductivity_name, - const std::string & magnetic_permeability_name, - const std::string & h_field_name); - - ~HFormulation() override = default; - - // Enable auxiliary calculation of J ∈ H(div) - void RegisterCurrentDensityAux(const std::string & j_field_name, - const std::string & external_j_field_name = "") override; - - // Enable auxiliary calculation of B ∈ H(div) - void RegisterMagneticFluxDensityAux(const std::string & b_field_name, - const std::string & external_b_field_name = "") override; - - // Enable auxiliary calculation of E ∈ H(curl) - void RegisterElectricFieldAux(const std::string & e_field_name, - const std::string & external_e_field_name = "") override; - - // Enable auxiliary calculation of H ∈ H(curl) - void RegisterMagneticFieldAux(const std::string & h_field_name, - const std::string & external_h_field_name = "") override; - - // Enable auxiliary calculation of F ∈ L2 - void RegisterLorentzForceDensityAux(const std::string & f_field_name, - const std::string & b_field_name, - const std::string & j_field_name) override; - - // Enable auxiliary calculation of P ∈ L2 - void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name, - const std::string & j_field_name) override; - void RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name) override; - - void RegisterCoefficients() override; - -protected: - const std::string _electric_conductivity_name; - const std::string & _electric_resistivity_name = platypus::HCurlFormulation::_alpha_coef_name; - const std::string & _magnetic_permeability_name = platypus::HCurlFormulation::_beta_coef_name; -}; -} // namespace platypus diff --git a/include/formulations/hcurl_formulation.h b/include/formulations/hcurl_formulation.h deleted file mode 100644 index ca3a7214..00000000 --- a/include/formulations/hcurl_formulation.h +++ /dev/null @@ -1,49 +0,0 @@ -#pragma once -#include "../common/pfem_extras.hpp" -#include "formulation.h" -#include "inputs.h" -#include "sources.h" - -namespace platypus -{ - -class HCurlFormulation : public TimeDomainEMFormulation -{ -public: - HCurlFormulation(std::string alpha_coef_name, - std::string beta_coef_name, - std::string h_curl_var_name); - - ~HCurlFormulation() override = default; - - void ConstructOperator() override; - - void ConstructJacobianPreconditioner() override; - - void ConstructJacobianSolver() override; - - void RegisterGridFunctions() override; - - void RegisterCoefficients() override; - -protected: - const std::string _alpha_coef_name; - const std::string _beta_coef_name; - const std::string _h_curl_var_name; -}; - -class CurlCurlEquationSystem : public TimeDependentEquationSystem -{ -public: - CurlCurlEquationSystem(const platypus::InputParameters & params); - - void Init(platypus::GridFunctions & gridfunctions, - const platypus::FESpaces & fespaces, - platypus::BCMap & bc_map, - platypus::Coefficients & coefficients) override; - void AddKernels() override; - - std::string _h_curl_var_name, _alpha_coef_name, _beta_coef_name, _dtalpha_coef_name; -}; - -} // namespace platypus diff --git a/include/formulations/magnetostatic_formulation.h b/include/formulations/magnetostatic_formulation.h deleted file mode 100644 index 1fe46ffe..00000000 --- a/include/formulations/magnetostatic_formulation.h +++ /dev/null @@ -1,37 +0,0 @@ -#pragma once -#include "../common/pfem_extras.hpp" -#include "inputs.h" -#include "statics_formulation.h" - -namespace platypus -{ - -class MagnetostaticFormulation : public platypus::StaticsFormulation -{ -public: - MagnetostaticFormulation(const std::string & magnetic_reluctivity_name, - std::string magnetic_permeability_name, - const std::string & magnetic_vector_potential_name); - - ~MagnetostaticFormulation() override = default; - - // Enable auxiliary calculation of B ∈ H(div) - void RegisterMagneticFluxDensityAux(const std::string & b_field_name, - const std::string & external_b_field_name = "") override; - - // Enable auxiliary calculation of H ∈ H(curl) - void RegisterMagneticFieldAux(const std::string & h_field_name, - const std::string & external_h_field_name = "") override; - - // Enable auxiliary calculation of F ∈ L2 - void RegisterLorentzForceDensityAux(const std::string & f_field_name, - const std::string & b_field_name, - const std::string & j_field_name) override; - - void RegisterCoefficients() override; - -protected: - const std::string _magnetic_permeability_name; - const std::string & _magnetic_reluctivity_name = platypus::StaticsFormulation::_alpha_coef_name; -}; -} // namespace platypus diff --git a/include/formulations/statics_formulation.h b/include/formulations/statics_formulation.h deleted file mode 100644 index 9b2da9b0..00000000 --- a/include/formulations/statics_formulation.h +++ /dev/null @@ -1,51 +0,0 @@ -#pragma once -#include "../common/pfem_extras.hpp" -#include "formulation.h" -#include "inputs.h" -#include "sources.h" - -namespace platypus -{ - -class StaticsFormulation : public SteadyStateEMFormulation -{ -public: - StaticsFormulation(std::string alpha_coef_name, std::string h_curl_var_name); - - ~StaticsFormulation() override = default; - - void ConstructJacobianPreconditioner() override; - - void ConstructJacobianSolver() override; - - void ConstructOperator() override; - - void RegisterGridFunctions() override; - - void RegisterCoefficients() override; - -protected: - const std::string _alpha_coef_name; - const std::string _h_curl_var_name; -}; - -class StaticsOperator : public ProblemOperator -{ -public: - StaticsOperator(platypus::Problem & problem, - std::string h_curl_var_name, - std::string stiffness_coef_name); - - ~StaticsOperator() override = default; - - void SetGridFunctions() override; - void Init(mfem::Vector & X) override; - void Solve(mfem::Vector & X) override; - -private: - std::string _h_curl_var_name, _stiffness_coef_name; - - mfem::Coefficient * _stiff_coef{nullptr}; // Stiffness Material Coefficient -}; - -} // namespace platypus diff --git a/include/formulations/steady_state_em_formulation.h b/include/formulations/steady_state_em_formulation.h deleted file mode 100644 index 60f6959c..00000000 --- a/include/formulations/steady_state_em_formulation.h +++ /dev/null @@ -1,16 +0,0 @@ -#pragma once -#include "em_formulation_interface.h" -#include "steady_state_problem_builder.h" - -namespace platypus -{ - -// Specifies output interfaces of a time-independent problem formulation. -class SteadyStateEMFormulation : public platypus::SteadyStateProblemBuilder, - public platypus::EMFormulationInterface -{ -public: - SteadyStateEMFormulation(); - ~SteadyStateEMFormulation() override = default; -}; -} // namespace platypus diff --git a/include/formulations/time_domain_em_formulation.h b/include/formulations/time_domain_em_formulation.h deleted file mode 100644 index 284b9e18..00000000 --- a/include/formulations/time_domain_em_formulation.h +++ /dev/null @@ -1,16 +0,0 @@ -#pragma once -#include "em_formulation_interface.h" -#include "time_domain_equation_system_problem_builder.h" - -namespace platypus -{ - -// Abstract Factory class of a time-domain EM formulation. -class TimeDomainEMFormulation : public platypus::TimeDomainEquationSystemProblemBuilder, - public platypus::EMFormulationInterface -{ -public: - TimeDomainEMFormulation(); - ~TimeDomainEMFormulation() override = default; -}; -} // namespace platypus diff --git a/src/formulations/CustomFormulation.C b/src/formulations/CustomFormulation.C index 76027e4d..c11130f1 100644 --- a/src/formulations/CustomFormulation.C +++ b/src/formulations/CustomFormulation.C @@ -10,9 +10,7 @@ CustomFormulation::validParams() } CustomFormulation::CustomFormulation(const InputParameters & parameters) - : MFEMFormulation(parameters) + : MFEMFormulation(parameters), + _formulation(std::make_shared()) { - formulation = std::make_shared(); } - -CustomFormulation::~CustomFormulation() {} diff --git a/src/formulations/MFEMFormulation.C b/src/formulations/MFEMFormulation.C index f1ef4897..3cf71952 100644 --- a/src/formulations/MFEMFormulation.C +++ b/src/formulations/MFEMFormulation.C @@ -13,5 +13,3 @@ MFEMFormulation::validParams() MFEMFormulation::MFEMFormulation(const InputParameters & parameters) : GeneralUserObject(parameters) { } - -MFEMFormulation::~MFEMFormulation() {} diff --git a/src/formulations/a_formulation.C b/src/formulations/a_formulation.C deleted file mode 100644 index 305c21bd..00000000 --- a/src/formulations/a_formulation.C +++ /dev/null @@ -1,145 +0,0 @@ -//* Solves: -//* ∇×(ν∇×A) + σdA/dt = Jᵉ -//* -//* in weak form -//* (ν∇×A, ∇×A') + (σdA/dt, A') - (Jᵉ, A') - <(ν∇×A)×n, A'> = 0 - -//* where: -//* reluctivity ν = 1/μ -//* electrical_conductivity σ=1/ρ -//* Magnetic vector potential A -//* Electric field, E = ρJᵉ -dA/dt -//* Magnetic flux density, B = ∇×A -//* Magnetic field H = ν∇×A -//* Current density J = Jᵉ -σdA/dt - -#include "a_formulation.h" - -#include - -namespace platypus -{ - -AFormulation::AFormulation(const std::string & magnetic_reluctivity_name, - std::string magnetic_permeability_name, - const std::string & electric_conductivity_name, - const std::string & magnetic_vector_potential_name) - : HCurlFormulation( - magnetic_reluctivity_name, electric_conductivity_name, magnetic_vector_potential_name), - _magnetic_permeability_name(std::move(magnetic_permeability_name)) -{ -} - -void -AFormulation::RegisterCurrentDensityAux(const std::string & j_field_name, - const std::string & external_j_field_name) -{ - //* Current density J = Jᵉ -σdA/dt - //* Induced electric field, Jind = -σdA/dt - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(j_field_name, - std::make_shared( - GetTimeDerivativeName(_h_curl_var_name), - j_field_name, - _electric_conductivity_name, - -1.0, - 1.0, - external_j_field_name)); -} - -void -AFormulation::RegisterMagneticFluxDensityAux(const std::string & b_field_name, - const std::string & external_b_field_name) -{ - //* Magnetic flux density, B = ∇×A - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(b_field_name, - std::make_shared(_h_curl_var_name, b_field_name)); -} - -void -AFormulation::RegisterElectricFieldAux(const std::string & e_field_name, - const std::string & external_e_field_name) -{ - //* Total electric field, E = ρJᵉ -dA/dt - //* Induced electric field, Eind = -dA/dt - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(e_field_name, - std::make_shared( - GetTimeDerivativeName(_h_curl_var_name), - e_field_name, - "_one", - -1.0, - 1.0, - external_e_field_name)); -} - -void -AFormulation::RegisterMagneticFieldAux(const std::string & h_field_name, - const std::string & external_h_field_name) -{ - //* Magnetic field H = ν∇×A - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(h_field_name, - std::make_shared( - _h_curl_var_name, h_field_name, _magnetic_reluctivity_name)); -} - -void -AFormulation::RegisterLorentzForceDensityAux(const std::string & f_field_name, - const std::string & b_field_name, - const std::string & j_field_name) -{ - //* Lorentz force density = J x B - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(f_field_name, - std::make_shared( - f_field_name, f_field_name, j_field_name, b_field_name)); - - auxsolvers.Get(f_field_name)->SetPriority(2); -} - -void -AFormulation::RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name, - const std::string & j_field_name) -{ - //* Joule heating density = E.J - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(p_field_name, - std::make_shared( - p_field_name, p_field_name, "", e_field_name, j_field_name)); - auxsolvers.Get(p_field_name)->SetPriority(2); -} - -void -AFormulation::RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name) -{ - //* Joule heating density = E.J - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - p_field_name, - std::make_shared( - p_field_name, p_field_name, _electric_conductivity_name, e_field_name, e_field_name)); - auxsolvers.Get(p_field_name)->SetPriority(2); -} - -void -AFormulation::RegisterCoefficients() -{ - platypus::Coefficients & coefficients = GetProblem()->_coefficients; - if (!coefficients._scalars.Has(_magnetic_permeability_name)) - { - MFEM_ABORT(_magnetic_permeability_name + " coefficient not found."); - } - if (!coefficients._scalars.Has(_electric_conductivity_name)) - { - MFEM_ABORT(_electric_conductivity_name + " coefficient not found."); - } - coefficients._scalars.Register( - _magnetic_reluctivity_name, - std::make_shared( - &_one_coef, coefficients._scalars.Get(_magnetic_permeability_name), fracFunc)); -} -} // namespace platypus diff --git a/src/formulations/av_formulation.C b/src/formulations/av_formulation.C deleted file mode 100644 index 10f9958c..00000000 --- a/src/formulations/av_formulation.C +++ /dev/null @@ -1,194 +0,0 @@ -// Solves: -// ∇×(ν∇×A) + σ(dA/dt + ∇ V) = Jᵉ -// ∇·(σ(dA/dt + ∇ V))= 0 - -// where -// Jᵉ ∈ H(div) source field -// A ∈ H(curl) -// V ∈ H1 - -//* in weak form -//* (ν∇×A, ∇×A') + (σ(dA/dt + ∇ V), A') - (Jᵉ, A') - <(ν∇×A) × n, A'> = 0 -//* (σ(dA/dt + ∇ V), ∇V') - <σ(dA/dt + ∇ V)·n, V'> =0 -//* -//* where: -//* reluctivity ν = 1/μ -//* electrical_conductivity σ=1/ρ -//* Magnetic vector potential A -//* Scalar electric potential V -//* Electric field, E = -dA/dt -∇V -//* Magnetic flux density, B = ∇×A -//* Magnetic field H = ν∇×A -//* Current density J = -σ(dA/dt + ∇ V) - -//* Either: -//* B.n (or E×n) at boundary: A×n (Dirichlet) -//* H×n at boundary: ν∇×A (Integrated) -//* -σ(dA/dt + ∇ V)·n (J·n, Neumann), V (potential, Dirichlet) - -#include "av_formulation.h" - -#include - -namespace platypus -{ - -AVFormulation::AVFormulation(std::string alpha_coef_name, - std::string inv_alpha_coef_name, - std::string beta_coef_name, - std::string vector_potential_name, - std::string scalar_potential_name) - : _alpha_coef_name(std::move(alpha_coef_name)), - _inv_alpha_coef_name(std::move(inv_alpha_coef_name)), - _beta_coef_name(std::move(beta_coef_name)), - _vector_potential_name(std::move(vector_potential_name)), - _scalar_potential_name(std::move(scalar_potential_name)) -{ -} - -void -AVFormulation::ConstructOperator() -{ - platypus::InputParameters av_system_params; - av_system_params.SetParam("VectorPotentialName", _vector_potential_name); - av_system_params.SetParam("ScalarPotentialName", _scalar_potential_name); - av_system_params.SetParam("AlphaCoefName", _alpha_coef_name); - av_system_params.SetParam("BetaCoefName", _beta_coef_name); - - auto equation_system = std::make_unique(av_system_params); - - GetProblem()->SetOperator(std::make_unique( - *GetProblem(), std::move(equation_system))); -} - -void -AVFormulation::RegisterGridFunctions() -{ - int & myid = GetProblem()->_myid; - platypus::GridFunctions & gridfunctions = GetProblem()->_gridfunctions; - - // Register default ParGridFunctions of state gridfunctions if not provided - if (!gridfunctions.Has(_vector_potential_name)) - { - if (myid == 0) - { - MFEM_WARNING(_vector_potential_name - << " not found in gridfunctions: building gridfunction from " - "defaults"); - } - AddFESpace(std::string("_HCurlFESpace"), std::string("ND_3D_P2")); - AddGridFunction(_vector_potential_name, std::string("_HCurlFESpace")); - } - - // Register default ParGridFunctions of state gridfunctions if not provided - if (!gridfunctions.Has(_scalar_potential_name)) - { - if (myid == 0) - { - MFEM_WARNING(_scalar_potential_name - << " not found in gridfunctions: building gridfunction from " - "defaults"); - } - AddFESpace(std::string("_H1FESpace"), std::string("H1_3D_P2")); - AddGridFunction(_scalar_potential_name, std::string("_H1FESpace")); - } - - // Register time derivatives - TimeDomainEquationSystemProblemBuilder::RegisterGridFunctions(); -}; - -void -AVFormulation::RegisterCoefficients() -{ - platypus::Coefficients & coefficients = GetProblem()->_coefficients; - if (!coefficients._scalars.Has(_inv_alpha_coef_name)) - { - MFEM_ABORT(_inv_alpha_coef_name + " coefficient not found."); - } - if (!coefficients._scalars.Has(_beta_coef_name)) - { - MFEM_ABORT(_beta_coef_name + " coefficient not found."); - } - - coefficients._scalars.Register( - _alpha_coef_name, - std::make_shared( - &_one_coef, coefficients._scalars.Get(_inv_alpha_coef_name), fracFunc)); -} - -AVEquationSystem::AVEquationSystem(const platypus::InputParameters & params) - : _a_name(params.GetParam("VectorPotentialName")), - _v_name(params.GetParam("ScalarPotentialName")), - _alpha_coef_name(params.GetParam("AlphaCoefName")), - _beta_coef_name(params.GetParam("BetaCoefName")), - _dtalpha_coef_name(std::string("dt_") + _alpha_coef_name), - _neg_beta_coef_name(std::string("negative_") + _beta_coef_name), - _neg_coef(-1.0) -{ -} - -void -AVEquationSystem::Init(platypus::GridFunctions & gridfunctions, - const platypus::FESpaces & fespaces, - platypus::BCMap & bc_map, - platypus::Coefficients & coefficients) -{ - coefficients._scalars.Register( - _dtalpha_coef_name, - std::make_shared( - &_dt_coef, coefficients._scalars.Get(_alpha_coef_name), prodFunc)); - - coefficients._scalars.Register( - _neg_beta_coef_name, - std::make_shared( - &_neg_coef, coefficients._scalars.Get(_beta_coef_name), prodFunc)); - - TimeDependentEquationSystem::Init(gridfunctions, fespaces, bc_map, coefficients); -} - -void -AVEquationSystem::AddKernels() -{ - AddTrialVariableNameIfMissing(_a_name); - std::string da_dt_name = GetTimeDerivativeName(_a_name); - AddTrialVariableNameIfMissing(_v_name); - std::string dv_dt_name = GetTimeDerivativeName(_v_name); - - // (α∇×A_{n}, ∇×dA'/dt) - platypus::InputParameters weak_curl_curl_params; - weak_curl_curl_params.SetParam("CoupledVariableName", _a_name); - weak_curl_curl_params.SetParam("CoefficientName", _alpha_coef_name); - AddKernel(da_dt_name, std::make_shared(weak_curl_curl_params)); - - // (αdt∇×dA/dt_{n+1}, ∇×dA'/dt) - platypus::InputParameters curl_curl_params; - curl_curl_params.SetParam("CoefficientName", _dtalpha_coef_name); - AddKernel(da_dt_name, std::make_shared(curl_curl_params)); - - // (βdA/dt_{n+1}, dA'/dt) - platypus::InputParameters vector_fe_mass_params; - vector_fe_mass_params.SetParam("CoefficientName", _beta_coef_name); - AddKernel(da_dt_name, std::make_shared(vector_fe_mass_params)); - - // (σ ∇ V, dA'/dt) - platypus::InputParameters mixed_vector_gradient_params; - mixed_vector_gradient_params.SetParam("CoefficientName", _beta_coef_name); - AddKernel(_v_name, - da_dt_name, - std::make_shared(mixed_vector_gradient_params)); - - // (σ ∇ V, ∇ V') - platypus::InputParameters diffusion_params; - diffusion_params.SetParam("CoefficientName", _beta_coef_name); - AddKernel(_v_name, std::make_shared(diffusion_params)); - - // (σdA/dt, ∇ V') - platypus::InputParameters vector_fe_weak_divergence_params; - vector_fe_weak_divergence_params.SetParam("CoefficientName", _beta_coef_name); - AddKernel( - da_dt_name, - _v_name, - std::make_shared(vector_fe_weak_divergence_params)); -} - -} // namespace platypus diff --git a/src/formulations/complex_a_formulation.C b/src/formulations/complex_a_formulation.C deleted file mode 100644 index 7359d4c2..00000000 --- a/src/formulations/complex_a_formulation.C +++ /dev/null @@ -1,137 +0,0 @@ -#include "complex_a_formulation.h" - -namespace platypus -{ - -ComplexAFormulation::ComplexAFormulation(const std::string & magnetic_reluctivity_name, - const std::string & electric_conductivity_name, - const std::string & dielectric_permittivity_name, - const std::string & frequency_coef_name, - const std::string & magnetic_vector_potential_complex_name, - const std::string & magnetic_vector_potential_real_name, - const std::string & magnetic_vector_potential_imag_name) - : ComplexMaxwellFormulation(magnetic_reluctivity_name, - electric_conductivity_name, - dielectric_permittivity_name, - frequency_coef_name, - magnetic_vector_potential_complex_name, - magnetic_vector_potential_real_name, - magnetic_vector_potential_imag_name) -{ -} - -// Enable auxiliary calculation of J ∈ H(div) -void -ComplexAFormulation::RegisterCurrentDensityAux(const std::string & j_field_real_name, - const std::string & j_field_imag_name, - const std::string & external_j_field_real_name, - const std::string & external_j_field_imag_name) -{ - //* Current density J = Jᵉ + σE - //* Induced current density Jind = σE = -iωσA - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - j_field_imag_name, - std::make_shared(_magnetic_vector_potential_real_name, - j_field_imag_name, - _loss_coef_name, - -1.0, - 1.0, - external_j_field_imag_name)); - auxsolvers.Register( - j_field_real_name, - std::make_shared(_magnetic_vector_potential_imag_name, - j_field_real_name, - _loss_coef_name, - 1.0, - 1.0, - external_j_field_real_name)); -}; - -void -ComplexAFormulation::RegisterMagneticFluxDensityAux(const std::string & b_field_real_name, - const std::string & b_field_imag_name, - const std::string & external_b_field_real_name, - const std::string & external_b_field_imag_name) -{ - //* Magnetic flux density B = curl A - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(b_field_real_name, - std::make_shared( - _magnetic_vector_potential_real_name, b_field_real_name)); - auxsolvers.Register(b_field_imag_name, - std::make_shared( - _magnetic_vector_potential_imag_name, b_field_imag_name)); -} - -void -ComplexAFormulation::RegisterElectricFieldAux(const std::string & e_field_real_name, - const std::string & e_field_imag_name, - const std::string & external_e_field_real_name, - const std::string & external_e_field_imag_name) -{ - //* Electric field E =-dA/dt=-iωA - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - e_field_imag_name, - std::make_shared(_magnetic_vector_potential_real_name, - e_field_imag_name, - "_angular_frequency", - -1.0, - 1.0, - external_e_field_imag_name)); - auxsolvers.Register( - e_field_real_name, - std::make_shared(_magnetic_vector_potential_imag_name, - e_field_real_name, - "_angular_frequency", - 1.0, - 1.0, - external_e_field_real_name)); -} - -// Enable auxiliary calculation of P ∈ L2 -// Time averaged Joule heating density E.J -void -ComplexAFormulation::RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_real_name, - const std::string & e_field_imag_name, - const std::string & j_field_real_name, - const std::string & j_field_imag_name) -{ - //* Time averaged Joule heating density = E.J - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(p_field_name, - std::make_shared( - p_field_name, - p_field_name, - _loss_coef_name, - _magnetic_vector_potential_real_name, - j_field_real_name, - _magnetic_vector_potential_imag_name, - j_field_imag_name, - true)); - auxsolvers.Get(p_field_name)->SetPriority(2); -} - -//* Time averaged Joule heating density = σ|E|^2 -void -ComplexAFormulation::RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_real_name, - const std::string & e_field_imag_name) -{ - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - p_field_name, - std::make_shared(p_field_name, - p_field_name, - _electric_conductivity_name, - e_field_real_name, - e_field_real_name, - e_field_imag_name, - e_field_imag_name, - true)); - auxsolvers.Get(p_field_name)->SetPriority(2); -} - -} // namespace platypus diff --git a/src/formulations/complex_e_formulation.C b/src/formulations/complex_e_formulation.C deleted file mode 100644 index d795fb0c..00000000 --- a/src/formulations/complex_e_formulation.C +++ /dev/null @@ -1,131 +0,0 @@ -#include "complex_e_formulation.h" - -namespace platypus -{ - -ComplexEFormulation::ComplexEFormulation(const std::string & magnetic_reluctivity_name, - const std::string & electric_conductivity_name, - const std::string & dielectric_permittivity_name, - const std::string & frequency_coef_name, - const std::string & e_field_complex_name, - const std::string & e_field_real_name, - const std::string & e_field_imag_name) - : platypus::ComplexMaxwellFormulation(magnetic_reluctivity_name, - electric_conductivity_name, - dielectric_permittivity_name, - frequency_coef_name, - e_field_complex_name, - e_field_real_name, - e_field_imag_name) -{ -} - -// Enable auxiliary calculation of J ∈ H(div) -void -ComplexEFormulation::RegisterCurrentDensityAux(const std::string & j_field_real_name, - const std::string & j_field_imag_name, - const std::string & external_j_field_real_name, - const std::string & external_j_field_imag_name) -{ - //* Current density J = Jᵉ + σE - //* Induced electric current, Jind = σE - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - j_field_real_name, - std::make_shared( - _electric_field_real_name, - j_field_real_name, - _electric_conductivity_name, - 1.0, - 1.0, - external_j_field_real_name, - platypus::InputParameters({{"Tolerance", float(1.0e-12)}, - {"MaxIter", (unsigned int)200}, - {"PrintLevel", 2}}))); // GetGlobalPrintLevel()}}))); - auxsolvers.Register( - j_field_imag_name, - std::make_shared( - _electric_field_imag_name, - j_field_imag_name, - _electric_conductivity_name, - 1.0, - 1.0, - external_j_field_imag_name, - platypus::InputParameters({{"Tolerance", float(1.0e-12)}, - {"MaxIter", (unsigned int)200}, - {"PrintLevel", 2}}))); // GetGlobalPrintLevel()}}))); -}; - -void -ComplexEFormulation::RegisterMagneticFluxDensityAux(const std::string & b_field_real_name, - const std::string & b_field_imag_name, - const std::string & external_b_field_real_name, - const std::string & external_b_field_imag_name) -{ - //* Magnetic flux density B = (i/ω) curl E - //* (∇×E = -dB/dt = -iωB) - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - b_field_imag_name, - std::make_shared(_electric_field_real_name, - b_field_imag_name, - "_inv_angular_frequency", - 1.0, - 1.0, - external_b_field_imag_name)); - auxsolvers.Register( - b_field_real_name, - std::make_shared(_electric_field_imag_name, - b_field_real_name, - "_inv_angular_frequency", - -1.0, - 1.0, - external_b_field_real_name)); -} - -//* Enable auxiliary calculation of P ∈ L2 -//* Time averaged Joule heating density = E.J -void -ComplexEFormulation::RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_real_name, - const std::string & e_field_imag_name, - const std::string & j_field_real_name, - const std::string & j_field_imag_name) -{ - - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - p_field_name, - std::make_shared(p_field_name, - p_field_name, - "", - _electric_field_real_name, - j_field_real_name, - _electric_field_imag_name, - j_field_imag_name, - true)); - auxsolvers.Get(p_field_name)->SetPriority(2); -} - -//* Enable auxiliary calculation of P ∈ L2 -//* Time averaged Joule heating density σ|E|^2 -void -ComplexEFormulation::RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_real_name, - const std::string & e_field_imag_name) -{ - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - p_field_name, - std::make_shared(p_field_name, - p_field_name, - _electric_conductivity_name, - _electric_field_real_name, - _electric_field_real_name, - _electric_field_imag_name, - _electric_field_imag_name, - true)); - auxsolvers.Get(p_field_name)->SetPriority(2); -} - -} // namespace platypus diff --git a/src/formulations/complex_maxwell_formulation.C b/src/formulations/complex_maxwell_formulation.C deleted file mode 100644 index b293ac9f..00000000 --- a/src/formulations/complex_maxwell_formulation.C +++ /dev/null @@ -1,224 +0,0 @@ -#include "complex_maxwell_formulation.h" - -#include - -namespace platypus -{ - -ComplexMaxwellFormulation::ComplexMaxwellFormulation(std::string alpha_coef_name, - std::string beta_coef_name, - std::string zeta_coef_name, - std::string frequency_coef_name, - std::string h_curl_var_complex_name, - std::string h_curl_var_real_name, - std::string h_curl_var_imag_name) - : _alpha_coef_name(std::move(alpha_coef_name)), - _beta_coef_name(std::move(beta_coef_name)), - _zeta_coef_name(std::move(zeta_coef_name)), - _frequency_coef_name(std::move(frequency_coef_name)), - _h_curl_var_complex_name(std::move(h_curl_var_complex_name)), - _h_curl_var_real_name(std::move(h_curl_var_real_name)), - _h_curl_var_imag_name(std::move(h_curl_var_imag_name)), - _mass_coef_name(std::string("maxwell_mass")), - _loss_coef_name(std::string("maxwell_loss")) -{ -} - -void -ComplexMaxwellFormulation::ConstructJacobianSolver() -{ - ConstructJacobianSolverWithOptions(SolverType::SUPER_LU); -} - -void -ComplexMaxwellFormulation::ConstructOperator() -{ - auto new_operator = std::make_unique(*GetProblem(), - _h_curl_var_complex_name, - _h_curl_var_real_name, - _h_curl_var_imag_name, - _alpha_coef_name, - _mass_coef_name, - _loss_coef_name); - - GetProblem()->SetOperator(std::move(new_operator)); -} - -void -ComplexMaxwellFormulation::RegisterGridFunctions() -{ - int & myid = GetProblem()->_myid; - platypus::GridFunctions & gridfunctions = GetProblem()->_gridfunctions; - platypus::FESpaces & fespaces = GetProblem()->_fespaces; - - // Register default ParGridFunctions of state gridfunctions if not provided - if (!gridfunctions.Has(_h_curl_var_real_name)) - { - if (myid == 0) - { - MFEM_WARNING(_h_curl_var_real_name << " not found in gridfunctions: building " - "gridfunction from defaults"); - } - AddFESpace(std::string("_HCurlFESpace"), std::string("ND_3D_P1")); - AddGridFunction(_h_curl_var_real_name, std::string("_HCurlFESpace")); - AddGridFunction(_h_curl_var_imag_name, std::string("_HCurlFESpace")); - } -} - -void -ComplexMaxwellFormulation::RegisterCoefficients() -{ - platypus::Coefficients & coefficients = GetProblem()->_coefficients; - - if (!coefficients._scalars.Has(_frequency_coef_name)) - { - MFEM_ABORT(_frequency_coef_name + " coefficient not found."); - } - if (!coefficients._scalars.Has("magnetic_permeability")) - { - MFEM_ABORT("Magnetic permeability coefficient not found."); - } - if (!coefficients._scalars.Has(_beta_coef_name)) - { - MFEM_ABORT(_beta_coef_name + " coefficient not found."); - } - if (!coefficients._scalars.Has(_zeta_coef_name)) - { - MFEM_ABORT(_zeta_coef_name + " coefficient not found."); - } - - _freq_coef = coefficients._scalars.Get(_frequency_coef_name); - - // define transformed - coefficients._scalars.Register( - "_angular_frequency", - std::make_shared(2.0 * M_PI * _freq_coef->constant)); - coefficients._scalars.Register( - "_neg_angular_frequency", - std::make_shared(-2.0 * M_PI * _freq_coef->constant)); - coefficients._scalars.Register( - "_angular_frequency_sq", - std::make_shared(pow(2.0 * M_PI * _freq_coef->constant, 2))); - coefficients._scalars.Register( - "_neg_angular_frequency_sq", - std::make_shared(-pow(2.0 * M_PI * _freq_coef->constant, 2))); - - coefficients._scalars.Register("_inv_angular_frequency", - std::make_shared( - 1.0, coefficients._scalars.GetRef("_angular_frequency"))); - - coefficients._scalars.Register(_mass_coef_name, - std::make_shared( - coefficients._scalars.Get("_neg_angular_frequency_sq"), - coefficients._scalars.Get(_zeta_coef_name), - prodFunc)); - - coefficients._scalars.Register(_loss_coef_name, - std::make_shared( - coefficients._scalars.Get("_angular_frequency"), - coefficients._scalars.Get(_beta_coef_name), - prodFunc)); - - coefficients._scalars.Register( - _alpha_coef_name, - std::make_shared( - &_one_coef, coefficients._scalars.Get("magnetic_permeability"), fracFunc)); -} - -ComplexMaxwellOperator::ComplexMaxwellOperator(platypus::Problem & problem, - std::string h_curl_var_complex_name, - std::string h_curl_var_real_name, - std::string h_curl_var_imag_name, - std::string stiffness_coef_name, - std::string mass_coef_name, - std::string loss_coef_name) - : ProblemOperator(problem), - _h_curl_var_complex_name(std::move(h_curl_var_complex_name)), - _h_curl_var_real_name(std::move(h_curl_var_real_name)), - _h_curl_var_imag_name(std::move(h_curl_var_imag_name)), - _stiffness_coef_name(std::move(stiffness_coef_name)), - _mass_coef_name(std::move(mass_coef_name)), - _loss_coef_name(std::move(loss_coef_name)) -{ -} - -void -ComplexMaxwellOperator::SetGridFunctions() -{ - _trial_var_names.push_back(_h_curl_var_real_name); - _trial_var_names.push_back(_h_curl_var_imag_name); - - ProblemOperator::SetGridFunctions(); - - _u = std::make_unique(_trial_variables.at(0)->ParFESpace()); - *_u = std::complex(0.0, 0.0); -}; - -void -ComplexMaxwellOperator::Init(mfem::Vector & X) -{ - ProblemOperator::Init(X); - - _stiff_coef = _problem._coefficients._scalars.Get(_stiffness_coef_name); - - if (_problem._coefficients._scalars.Has(_mass_coef_name)) - _mass_coef = _problem._coefficients._scalars.Get(_mass_coef_name); - if (_problem._coefficients._scalars.Has(_loss_coef_name)) - _loss_coef = _problem._coefficients._scalars.Get(_loss_coef_name); -} - -void -ComplexMaxwellOperator::Solve(mfem::Vector & X) -{ - mfem::OperatorHandle jac; - mfem::Vector u, rhs; - mfem::OperatorHandle pc_op; - - mfem::Vector zero_vec(3); - zero_vec = 0.0; - mfem::VectorConstantCoefficient zero_coef(zero_vec); - - mfem::ParSesquilinearForm sqlf(_u->ParFESpace(), _conv); - sqlf.AddDomainIntegrator(new mfem::CurlCurlIntegrator(*_stiff_coef), nullptr); - if (_mass_coef) - { - sqlf.AddDomainIntegrator(new mfem::VectorFEMassIntegrator(*_mass_coef), nullptr); - } - if (_loss_coef) - { - sqlf.AddDomainIntegrator(nullptr, new mfem::VectorFEMassIntegrator(*_loss_coef)); - } - - mfem::ParLinearForm lf_real(_u->ParFESpace()); - mfem::ParLinearForm lf_imag(_u->ParFESpace()); - lf_real = 0.0; - lf_imag = 0.0; - - _problem._sources.Apply(&lf_real); - - mfem::ParComplexLinearForm lf(_u->ParFESpace(), _conv); - _problem._bc_map.ApplyEssentialBCs( - _h_curl_var_complex_name, _ess_bdr_tdofs, *_u, _problem._pmesh.get()); - _problem._bc_map.ApplyIntegratedBCs(_h_curl_var_complex_name, lf, _problem._pmesh.get()); - _problem._bc_map.ApplyIntegratedBCs(_h_curl_var_complex_name, sqlf, _problem._pmesh.get()); - - sqlf.Assemble(); - sqlf.Finalize(); - - lf.Assemble(); - lf.real() += lf_real; - lf.imag() += lf_imag; - - sqlf.FormLinearSystem(_ess_bdr_tdofs, *_u, lf, jac, u, rhs); - - auto * jac_z = jac.As()->GetSystemMatrix(); - - _problem._jacobian_solver->SetOperator(*jac_z); - _problem._jacobian_solver->Mult(rhs, u); - sqlf.RecoverFEMSolution(u, lf, *_u); - - _problem._gridfunctions.GetRef(_trial_var_names.at(0)) = _u->real(); - _problem._gridfunctions.GetRef(_trial_var_names.at(1)) = _u->imag(); -} - -} // namespace platypus diff --git a/src/formulations/dual_formulation.C b/src/formulations/dual_formulation.C deleted file mode 100644 index 02632703..00000000 --- a/src/formulations/dual_formulation.C +++ /dev/null @@ -1,241 +0,0 @@ -// Solves the equations -// ∇⋅s0 = 0 -// ∇×(αv) - βu = s0 -// dv/dt = -∇×u - -// where -// s0 ∈ H(div) source field -// v ∈ H(div) -// u ∈ H(curl) -// p ∈ H1 - -// Weak form (Space discretisation) -// -(s0, ∇ p') + = 0 -// (αv, ∇×u') - (βu, u') - (s0, u') - <(αv) × n, u'> = 0 -// (dv/dt, v') + (∇×u, v') = 0 - -// Time discretisation using implicit scheme: -// Unknowns -// s0_{n+1} ∈ H(div) source field, where s0 = -β∇p -// dv/dt_{n+1} ∈ H(div) -// u_{n+1} ∈ H(curl) -// p_{n+1} ∈ H1 - -// Fully discretised equations -// -(s0_{n+1}, ∇ p') + = 0 -// (αv_{n}, ∇×u') - (αdt∇×u_{n+1}, ∇×u') - (βu_{n+1}, u') - (s0_{n+1}, u') - -// <(αv) × n, u'> = 0 -// (dv/dt_{n+1}, v') + (∇×u_{n+1}, v') = 0 -// using -// v_{n+1} = v_{n} + dt dv/dt_{n+1} = v_{n} - dt ∇×u_{n+1} - -// Rewritten as -// a0(p_{n+1}, p') = b0(p') -// a1(u_{n+1}, u') = b1(u') -// dv/dt_{n+1} = -∇×u - -// where -// a0(p, p') = (β ∇ p, ∇ p') -// b0(p') = -// a1(u, u') = (βu, u') + (αdt∇×u, ∇×u') -// b1(u') = (s0_{n+1}, u') + (αv_{n}, ∇×u') + <(αdt∇×u_{n+1}) × n, u'> - -#include "dual_formulation.h" - -#include - -namespace platypus -{ - -DualFormulation::DualFormulation(std::string alpha_coef_name, - std::string beta_coef_name, - std::string h_curl_var_name, - std::string h_div_var_name) - : _alpha_coef_name(std::move(alpha_coef_name)), - _beta_coef_name(std::move(beta_coef_name)), - _h_curl_var_name(std::move(h_curl_var_name)), - _h_div_var_name(std::move(h_div_var_name)) -{ -} - -void -DualFormulation::ConstructJacobianPreconditioner() -{ - auto precond = - std::make_shared(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0)); - - precond->SetSingularProblem(); - precond->SetPrintLevel(-1); - - GetProblem()->_jacobian_preconditioner = precond; -} - -void -DualFormulation::ConstructJacobianSolver() -{ - ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG, {._max_iteration = 1000}); -} - -void -DualFormulation::ConstructOperator() -{ - platypus::InputParameters weak_form_params; - weak_form_params.SetParam("HCurlVarName", _h_curl_var_name); - weak_form_params.SetParam("HDivVarName", _h_div_var_name); - weak_form_params.SetParam("AlphaCoefName", _alpha_coef_name); - weak_form_params.SetParam("BetaCoefName", _beta_coef_name); - - auto equation_system = std::make_unique(weak_form_params); - - GetProblem()->SetOperator( - std::make_unique(*GetProblem(), std::move(equation_system))); -} - -void -DualFormulation::RegisterGridFunctions() -{ - int & myid = GetProblem()->_myid; - platypus::GridFunctions & gridfunctions = GetProblem()->_gridfunctions; - - // Register default ParGridFunctions of state gridfunctions if not provided - if (!gridfunctions.Has(_h_curl_var_name)) - { - if (myid == 0) - { - MFEM_WARNING(_h_curl_var_name << " not found in gridfunctions: building " - "gridfunction from defaults"); - } - AddFESpace(std::string("_HCurlFESpace"), std::string("ND_3D_P2")); - AddGridFunction(_h_curl_var_name, std::string("_HCurlFESpace")); - } - - // Register default ParGridFunctions of state gridfunctions if not provided - if (!gridfunctions.Has(_h_div_var_name)) - { - if (myid == 0) - { - MFEM_WARNING(_h_div_var_name << " not found in gridfunctions: building " - "gridfunction from defaults"); - } - AddFESpace(std::string("_HDivFESpace"), std::string("RT_3D_P3")); - AddGridFunction(_h_div_var_name, std::string("_HDivFESpace")); - } - - // Register time derivatives - TimeDomainEquationSystemProblemBuilder::RegisterGridFunctions(); -}; - -void -DualFormulation::RegisterCoefficients() -{ - platypus::Coefficients & coefficients = GetProblem()->_coefficients; - - if (!coefficients._scalars.Has(_alpha_coef_name)) - { - MFEM_ABORT(_alpha_coef_name + " coefficient not found."); - } - if (!coefficients._scalars.Has(_beta_coef_name)) - { - MFEM_ABORT(_beta_coef_name + " coefficient not found."); - } -} - -WeakCurlEquationSystem::WeakCurlEquationSystem(const platypus::InputParameters & params) - : _h_curl_var_name(params.GetParam("HCurlVarName")), - _h_div_var_name(params.GetParam("HDivVarName")), - _alpha_coef_name(params.GetParam("AlphaCoefName")), - _beta_coef_name(params.GetParam("BetaCoefName")), - _dtalpha_coef_name(std::string("dt_") + _alpha_coef_name) -{ -} - -void -WeakCurlEquationSystem::Init(platypus::GridFunctions & gridfunctions, - const platypus::FESpaces & fespaces, - platypus::BCMap & bc_map, - platypus::Coefficients & coefficients) -{ - coefficients._scalars.Register( - _dtalpha_coef_name, - std::make_shared( - &_dt_coef, coefficients._scalars.Get(_alpha_coef_name), prodFunc)); - TimeDependentEquationSystem::Init(gridfunctions, fespaces, bc_map, coefficients); -} - -void -WeakCurlEquationSystem::AddKernels() -{ - AddTrialVariableNameIfMissing(_h_curl_var_name); - AddTrialVariableNameIfMissing(_h_div_var_name); - std::string dh_curl_var_dt = GetTimeDerivativeName(_h_curl_var_name); - - // (αv_{n}, ∇×u') - platypus::InputParameters weak_curl_params; - weak_curl_params.SetParam("HCurlVarName", _h_curl_var_name); - weak_curl_params.SetParam("HDivVarName", _h_div_var_name); - weak_curl_params.SetParam("CoefficientName", _alpha_coef_name); - AddKernel(_h_curl_var_name, std::make_shared(weak_curl_params)); - - // (αdt∇×u_{n+1}, ∇×u') - platypus::InputParameters curl_curl_params; - curl_curl_params.SetParam("CoefficientName", _dtalpha_coef_name); - AddKernel(_h_curl_var_name, std::make_shared(curl_curl_params)); - - // (βu_{n+1}, u') - platypus::InputParameters vector_fe_mass_params; - vector_fe_mass_params.SetParam("CoefficientName", _beta_coef_name); - AddKernel(_h_curl_var_name, - std::make_shared(vector_fe_mass_params)); -} - -void -DualOperator::Init(mfem::Vector & X) -{ - TimeDomainEquationSystemProblemOperator::Init(X); - auto * eqs = dynamic_cast(GetEquationSystem()); - - _h_curl_var_name = eqs->_h_curl_var_name; - _h_div_var_name = eqs->_h_div_var_name; - - _u = _problem._gridfunctions.Get(_h_curl_var_name); - _dv = _problem._gridfunctions.Get(GetTimeDerivativeName(_h_div_var_name)); - - _h_curl_fe_space = _u->ParFESpace(); - _h_div_fe_space = _dv->ParFESpace(); - - _curl = std::make_unique(_h_curl_fe_space, _h_div_fe_space); - _curl->AddDomainInterpolator(new mfem::CurlInterpolator); - _curl->Assemble(); -} - -void -DualOperator::ImplicitSolve(const double dt, const mfem::Vector & X, mfem::Vector & dX_dt) -{ - TimeDomainEquationSystemProblemOperator::ImplicitSolve(dt, X, dX_dt); - // Subtract off contribution from source - _problem._sources.SubtractSources(_u); - // dv/dt_{n+1} = -∇×u - _curl->Mult(*_u, *_dv); - *_dv *= -1.0; -} - -void -DualOperator::SetGridFunctions() -{ - TimeDomainEquationSystemProblemOperator::SetGridFunctions(); - // Blocks for solution vector are smaller than the operator size - // for DualOperator, as curl is stored separately. - // Block operator only has the HCurl TrueVSize; - _block_true_offsets.SetSize(_trial_variables.size()); - _block_true_offsets[0] = 0; - for (unsigned int ind = 0; ind < _trial_variables.size() - 1; ++ind) - { - _block_true_offsets[ind + 1] = _trial_variables.at(ind)->ParFESpace()->TrueVSize(); - } - _block_true_offsets.PartialSum(); - - _true_x.Update(_block_true_offsets); - _true_rhs.Update(_block_true_offsets); -} - -} // namespace platypus diff --git a/src/formulations/e_formulation.C b/src/formulations/e_formulation.C deleted file mode 100644 index 98667651..00000000 --- a/src/formulations/e_formulation.C +++ /dev/null @@ -1,91 +0,0 @@ -//* Solves: -//* ∇×(ν∇×E) + σdE/dt = -dJᵉ/dt -//* -//* in weak form -//* (ν∇×E, ∇×E') + (σdE/dt, E') + (dJᵉ/dt, E') - <(ν∇×E)×n, E'> = 0 - -//* where: -//* reluctivity ν = 1/μ -//* electrical_conductivity σ=1/ρ -//* Electric Field E -//* Current density J = σE -//* Magnetic flux density, dB/dt = -∇×E -//* Magnetic field dH/dt = -ν∇×E - -#include "e_formulation.h" - -#include - -namespace platypus -{ - -EFormulation::EFormulation(const std::string & magnetic_reluctivity_name, - std::string magnetic_permeability_name, - const std::string & electric_conductivity_name, - const std::string & e_field_name) - : HCurlFormulation(magnetic_reluctivity_name, electric_conductivity_name, e_field_name), - _magnetic_permeability_name(std::move(magnetic_permeability_name)) -{ -} - -void -EFormulation::RegisterCoefficients() -{ - platypus::Coefficients & coefficients = GetProblem()->_coefficients; - if (!coefficients._scalars.Has(_magnetic_permeability_name)) - { - MFEM_ABORT(_magnetic_permeability_name + " coefficient not found."); - } - if (!coefficients._scalars.Has(_electric_conductivity_name)) - { - MFEM_ABORT(_electric_conductivity_name + " coefficient not found."); - } - coefficients._scalars.Register( - _magnetic_reluctivity_name, - std::make_shared( - &_one_coef, coefficients._scalars.Get(_magnetic_permeability_name), fracFunc)); -} - -void -EFormulation::RegisterCurrentDensityAux(const std::string & j_field_name, - const std::string & external_j_field_name) -{ - //* Current density J = Jᵉ + σE - //* Induced electric field, Jind = σE - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - j_field_name, - std::make_shared(_h_curl_var_name, - j_field_name, - _electric_conductivity_name, - 1.0, - 1.0, - external_j_field_name)); -} - -void -EFormulation::RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name, - const std::string & j_field_name) -{ - //* Joule heating density = E.J - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(p_field_name, - std::make_shared( - p_field_name, p_field_name, "", e_field_name, j_field_name)); - auxsolvers.Get(p_field_name)->SetPriority(2); -} -void -EFormulation::RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name) -{ - //* Joule heating density = E.J - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - p_field_name, - std::make_shared( - p_field_name, p_field_name, _electric_conductivity_name, e_field_name, e_field_name)); - auxsolvers.Get(p_field_name)->SetPriority(2); -} - -} // namespace platypus diff --git a/src/formulations/eb_dual_formulation.C b/src/formulations/eb_dual_formulation.C deleted file mode 100644 index 03a0d6d4..00000000 --- a/src/formulations/eb_dual_formulation.C +++ /dev/null @@ -1,91 +0,0 @@ -#include "eb_dual_formulation.h" - -#include - -namespace platypus -{ - -EBDualFormulation::EBDualFormulation(const std::string & magnetic_reluctivity_name, - std::string magnetic_permeability_name, - const std::string & electric_conductivity_name, - const std::string & e_field_name, - const std::string & b_field_name) - : DualFormulation( - magnetic_reluctivity_name, electric_conductivity_name, e_field_name, b_field_name), - _magnetic_permeability_name(std::move(magnetic_permeability_name)) -{ -} - -void -EBDualFormulation::RegisterCurrentDensityAux(const std::string & j_field_name, - const std::string & external_j_field_name) -{ - //* Current density J = Jᵉ + σE - //* Induced electric field, Jind = σE - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - j_field_name, - std::make_shared(_h_curl_var_name, - j_field_name, - _electric_conductivity_name, - 1.0, - 1.0, - external_j_field_name)); -} - -void -EBDualFormulation::RegisterLorentzForceDensityAux(const std::string & f_field_name, - const std::string & b_field_name, - const std::string & j_field_name) -{ - //* Lorentz force density = J x B - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(f_field_name, - std::make_shared( - f_field_name, f_field_name, j_field_name, b_field_name)); - - auxsolvers.Get(f_field_name)->SetPriority(2); -} - -void -EBDualFormulation::RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name, - const std::string & j_field_name) -{ - //* Joule heating density = E.J - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(p_field_name, - std::make_shared( - p_field_name, p_field_name, "", e_field_name, j_field_name)); - auxsolvers.Get(p_field_name)->SetPriority(2); -} - -void -EBDualFormulation::RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name) -{ - //* Joule heating density = E.J - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - p_field_name, - std::make_shared( - p_field_name, p_field_name, _electric_conductivity_name, e_field_name, e_field_name)); - auxsolvers.Get(p_field_name)->SetPriority(2); -} - -void -EBDualFormulation::RegisterCoefficients() -{ - platypus::Coefficients & coefficients = GetProblem()->_coefficients; - if (!coefficients._scalars.Has(_magnetic_permeability_name)) - { - MFEM_ABORT(_magnetic_permeability_name + " coefficient not found."); - } - coefficients._scalars.Register( - _magnetic_reluctivity_name, - std::make_shared( - &_one_coef, coefficients._scalars.Get(_magnetic_permeability_name), fracFunc)); - DualFormulation::RegisterCoefficients(); -} - -} // namespace platypus diff --git a/src/formulations/frequency_domain_em_formulation.C b/src/formulations/frequency_domain_em_formulation.C deleted file mode 100644 index d6d17290..00000000 --- a/src/formulations/frequency_domain_em_formulation.C +++ /dev/null @@ -1,8 +0,0 @@ -#include "frequency_domain_em_formulation.h" - -namespace platypus -{ - -FrequencyDomainEMFormulation::FrequencyDomainEMFormulation() = default; - -} // namespace platypus diff --git a/src/formulations/h_formulation.C b/src/formulations/h_formulation.C deleted file mode 100644 index 75d048c0..00000000 --- a/src/formulations/h_formulation.C +++ /dev/null @@ -1,131 +0,0 @@ -//* Solves: -//* ∇×(ρ∇×H) + μdH/dt = -dBᵉ/dt -//* -//* in weak form -//* (ρ∇×H, ∇×v) + (μdH/dt, v) + (dBᵉ/dt, v) - <(ρ∇×H)×n, v> = 0 - -//* where: -//* magnetic permeability μ = 1/ν -//* electrical resistivity ρ=1/σ -//* Electric field, E = ρ∇×H -//* Magnetic flux density, B = Bᵉ + μH -//* Current density J = ∇×H - -#include "h_formulation.h" - -#include - -namespace platypus -{ - -HFormulation::HFormulation(const std::string & electric_resistivity_name, - std::string electric_conductivity_name, - const std::string & magnetic_permeability_name, - const std::string & h_field_name) - : HCurlFormulation(electric_resistivity_name, magnetic_permeability_name, h_field_name), - _electric_conductivity_name(std::move(electric_conductivity_name)) -{ -} - -void -HFormulation::RegisterCurrentDensityAux(const std::string & j_field_name, - const std::string & external_j_field_name) -{ - //* Current density J = ∇×H - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(j_field_name, - std::make_shared(_h_curl_var_name, j_field_name)); -} - -void -HFormulation::RegisterMagneticFluxDensityAux(const std::string & b_field_name, - const std::string & external_b_field_name) -{ - //* Magnetic flux density, B = Bᵉ + μH - //* Induced flux density, B = μH - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - b_field_name, - std::make_shared(_h_curl_var_name, - b_field_name, - _magnetic_permeability_name, - 1.0, - 1.0, - external_b_field_name)); -} - -void -HFormulation::RegisterElectricFieldAux(const std::string & e_field_name, - const std::string & external_e_field_name) -{ - //* Electric field, E = ρ∇×H - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(e_field_name, - std::make_shared( - _h_curl_var_name, e_field_name, _electric_resistivity_name)); -} - -void -HFormulation::RegisterMagneticFieldAux(const std::string & h_field_name, - const std::string & external_h_field_name) - -{ - //* Magnetic field H is a state variable; no additional calculation needed -} - -void -HFormulation::RegisterLorentzForceDensityAux(const std::string & f_field_name, - const std::string & b_field_name, - const std::string & j_field_name) -{ - //* Lorentz force density = J x B - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(f_field_name, - std::make_shared( - f_field_name, f_field_name, j_field_name, b_field_name)); - - auxsolvers.Get(f_field_name)->SetPriority(2); -} - -void -HFormulation::RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name, - const std::string & j_field_name) -{ - //* Joule heating density = E.J - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(p_field_name, - std::make_shared( - p_field_name, p_field_name, "", e_field_name, j_field_name)); - auxsolvers.Get(p_field_name)->SetPriority(2); -} - -void -HFormulation::RegisterJouleHeatingDensityAux(const std::string & p_field_name, - const std::string & e_field_name) -{ - //* Joule heating density = E.J - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - p_field_name, - std::make_shared( - p_field_name, p_field_name, _electric_conductivity_name, e_field_name, e_field_name)); - auxsolvers.Get(p_field_name)->SetPriority(2); -} - -void -HFormulation::RegisterCoefficients() -{ - platypus::Coefficients & coefficients = GetProblem()->_coefficients; - if (!coefficients._scalars.Has(_electric_conductivity_name)) - { - MFEM_ABORT(_electric_conductivity_name + " coefficient not found."); - } - coefficients._scalars.Register( - _electric_resistivity_name, - std::make_shared( - &_one_coef, coefficients._scalars.Get(_electric_conductivity_name), fracFunc)); - HCurlFormulation::RegisterCoefficients(); -} - -} // namespace platypus diff --git a/src/formulations/hcurl_formulation.C b/src/formulations/hcurl_formulation.C deleted file mode 100644 index fbdfe5f3..00000000 --- a/src/formulations/hcurl_formulation.C +++ /dev/null @@ -1,167 +0,0 @@ -// Solves the equations -// ∇⋅s0 = 0 -// ∇×(α∇×u) + βdu/dt = s0 - -// where -// s0 ∈ H(div) source field -// u ∈ H(curl) -// p ∈ H1 - -// Dirichlet boundaries constrain du/dt -// Integrated boundaries constrain (α∇×u) × n - -// Weak form (Space discretisation) -// -(s0, ∇ p') + = 0 -// (α∇×u, ∇×u') + (βdu/dt, u') - (s0, u') - <(α∇×u) × n, u'> = 0 - -// Time discretisation using implicit scheme: -// Unknowns -// s0_{n+1} ∈ H(div) source field, where s0 = -β∇p -// du/dt_{n+1} ∈ H(curl) -// p_{n+1} ∈ H1 - -// Fully discretised equations -// -(s0_{n+1}, ∇ p') + = 0 -// (α∇×u_{n}, ∇×u') + (αdt∇×du/dt_{n+1}, ∇×u') + (βdu/dt_{n+1}, u') -// - (s0_{n+1}, u') - <(α∇×u_{n+1}) × n, u'> = 0 -// using -// u_{n+1} = u_{n} + dt du/dt_{n+1} - -// Rewritten as -// a0(p_{n+1}, p') = b0(p') -// a1(du/dt_{n+1}, u') = b1(u') - -// where -// a0(p, p') = (β ∇ p, ∇ p') -// b0(p') = -// a1(u, u') = (βu, u') + (αdt∇×u, ∇×u') -// b1(u') = (s0_{n+1}, u') - (α∇×u_{n}, ∇×u') + <(α∇×u_{n+1}) × n, u'> -#include "hcurl_formulation.h" - -#include - -namespace platypus -{ - -HCurlFormulation::HCurlFormulation(std::string alpha_coef_name, - std::string beta_coef_name, - std::string h_curl_var_name) - : _alpha_coef_name(std::move(alpha_coef_name)), - _beta_coef_name(std::move(beta_coef_name)), - _h_curl_var_name(std::move(h_curl_var_name)) -{ -} - -void -HCurlFormulation::ConstructOperator() -{ - platypus::InputParameters weak_form_params; - weak_form_params.SetParam("HCurlVarName", _h_curl_var_name); - weak_form_params.SetParam("AlphaCoefName", _alpha_coef_name); - weak_form_params.SetParam("BetaCoefName", _beta_coef_name); - - auto equation_system = std::make_unique(weak_form_params); - - GetProblem()->SetOperator(std::make_unique( - *GetProblem(), std::move(equation_system))); -} - -void -HCurlFormulation::ConstructJacobianPreconditioner() -{ - auto precond = - std::make_shared(GetProblem()->GetEquationSystem()->_test_pfespaces.at(0)); - - precond->SetSingularProblem(); - precond->SetPrintLevel(-1); - - GetProblem()->_jacobian_preconditioner = precond; -} - -void -HCurlFormulation::ConstructJacobianSolver() -{ - ConstructJacobianSolverWithOptions(SolverType::HYPRE_PCG); -} - -void -HCurlFormulation::RegisterGridFunctions() -{ - int & myid = GetProblem()->_myid; - platypus::GridFunctions & gridfunctions = GetProblem()->_gridfunctions; - platypus::FESpaces & fespaces = GetProblem()->_fespaces; - - // Register default ParGridFunctions of state gridfunctions if not provided - if (!gridfunctions.Has(_h_curl_var_name)) - { - if (myid == 0) - { - MFEM_WARNING(_h_curl_var_name << " not found in gridfunctions: building " - "gridfunction from defaults"); - } - AddFESpace(std::string("_HCurlFESpace"), std::string("ND_3D_P2")); - AddGridFunction(_h_curl_var_name, std::string("_HCurlFESpace")); - }; - // Register time derivatives - TimeDomainEquationSystemProblemBuilder::RegisterGridFunctions(); -}; - -CurlCurlEquationSystem::CurlCurlEquationSystem(const platypus::InputParameters & params) - : _h_curl_var_name(params.GetParam("HCurlVarName")), - _alpha_coef_name(params.GetParam("AlphaCoefName")), - _beta_coef_name(params.GetParam("BetaCoefName")), - _dtalpha_coef_name(std::string("dt_") + _alpha_coef_name) -{ -} - -void -CurlCurlEquationSystem::Init(platypus::GridFunctions & gridfunctions, - const platypus::FESpaces & fespaces, - platypus::BCMap & bc_map, - platypus::Coefficients & coefficients) -{ - coefficients._scalars.Register( - _dtalpha_coef_name, - std::make_shared( - &_dt_coef, coefficients._scalars.Get(_alpha_coef_name), prodFunc)); - TimeDependentEquationSystem::Init(gridfunctions, fespaces, bc_map, coefficients); -} - -void -CurlCurlEquationSystem::AddKernels() -{ - AddTrialVariableNameIfMissing(_h_curl_var_name); - std::string dh_curl_var_dt = GetTimeDerivativeName(_h_curl_var_name); - - // (α∇×u_{n}, ∇×u') - platypus::InputParameters weak_curl_curl_params; - weak_curl_curl_params.SetParam("CoupledVariableName", _h_curl_var_name); - weak_curl_curl_params.SetParam("CoefficientName", _alpha_coef_name); - AddKernel(dh_curl_var_dt, std::make_shared(weak_curl_curl_params)); - - // (αdt∇×du/dt_{n+1}, ∇×u') - platypus::InputParameters curl_curl_params; - curl_curl_params.SetParam("CoefficientName", _dtalpha_coef_name); - AddKernel(dh_curl_var_dt, std::make_shared(curl_curl_params)); - - // (βdu/dt_{n+1}, u') - platypus::InputParameters vector_fe_mass_params; - vector_fe_mass_params.SetParam("CoefficientName", _beta_coef_name); - AddKernel(dh_curl_var_dt, std::make_shared(vector_fe_mass_params)); -} - -void -HCurlFormulation::RegisterCoefficients() -{ - platypus::Coefficients & coefficients = GetProblem()->_coefficients; - if (!coefficients._scalars.Has(_alpha_coef_name)) - { - MFEM_ABORT(_alpha_coef_name + " coefficient not found."); - } - if (!coefficients._scalars.Has(_beta_coef_name)) - { - MFEM_ABORT(_beta_coef_name + " coefficient not found."); - } -} - -} // namespace platypus diff --git a/src/formulations/magnetostatic_formulation.C b/src/formulations/magnetostatic_formulation.C deleted file mode 100644 index a921161b..00000000 --- a/src/formulations/magnetostatic_formulation.C +++ /dev/null @@ -1,84 +0,0 @@ -//* Solves: -//* ∇×(ν∇×A) = Jᵉ -//* -//* in weak form -//* (ν∇×A, ∇×A') - (Jᵉ, A') - <(ν∇×A)×n, A'> = 0 - -//* where: -//* reluctivity ν = 1/μ -//* Magnetic vector potential A -//* Electric field, E = ρJᵉ -//* Magnetic flux density, B = ∇×A -//* Magnetic field H = ν∇×A -//* Current density J = Jᵉ - -#include "magnetostatic_formulation.h" - -#include - -namespace platypus -{ - -MagnetostaticFormulation::MagnetostaticFormulation( - const std::string & magnetic_reluctivity_name, - std::string magnetic_permeability_name, - const std::string & magnetic_vector_potential_name) - : StaticsFormulation(magnetic_reluctivity_name, magnetic_vector_potential_name), - _magnetic_permeability_name(std::move(magnetic_permeability_name)) -{ -} - -void -MagnetostaticFormulation::RegisterMagneticFluxDensityAux(const std::string & b_field_name, - const std::string & external_b_field_name) -{ - //* Magnetic flux density, B = ∇×A - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(b_field_name, - std::make_shared(_h_curl_var_name, b_field_name)); -} - -void -MagnetostaticFormulation::RegisterMagneticFieldAux(const std::string & h_field_name, - const std::string & external_h_field_name) -{ - //* Magnetic field H = ν∇×A - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register( - h_field_name, - std::make_shared(_h_curl_var_name, - h_field_name, - _magnetic_reluctivity_name, - 1.0, - 1.0, - external_h_field_name)); -} - -void -MagnetostaticFormulation::RegisterLorentzForceDensityAux(const std::string & f_field_name, - const std::string & b_field_name, - const std::string & j_field_name) -{ - //* Lorentz force density = J x B - platypus::AuxSolvers & auxsolvers = GetProblem()->_postprocessors; - auxsolvers.Register(f_field_name, - std::make_shared( - f_field_name, f_field_name, j_field_name, b_field_name)); - - auxsolvers.Get(f_field_name)->SetPriority(2); -} - -void -MagnetostaticFormulation::RegisterCoefficients() -{ - platypus::Coefficients & coefficients = GetProblem()->_coefficients; - if (!coefficients._scalars.Has(_magnetic_permeability_name)) - { - MFEM_ABORT(_magnetic_permeability_name + " coefficient not found."); - } - coefficients._scalars.Register( - _magnetic_reluctivity_name, - std::make_shared( - &_one_coef, coefficients._scalars.Get(_magnetic_permeability_name), fracFunc)); -} -} // namespace platypus diff --git a/src/formulations/statics_formulation.C b/src/formulations/statics_formulation.C deleted file mode 100644 index 44f8ed8b..00000000 --- a/src/formulations/statics_formulation.C +++ /dev/null @@ -1,159 +0,0 @@ -// Solves the equations -// ∇⋅s0 = 0 -// ∇×(α∇×u) = s0 - -// where -// s0 ∈ H(div) source field -// u ∈ H(curl) - -// Dirichlet boundaries constrain u -// Integrated boundaries constrain (α∇×u) × n - -// Weak form (Space discretisation) -// (α∇×u, ∇×u') - (s0, u') - <(α∇×u) × n, u'> = 0 - -// Time discretisation using implicit scheme: -// Unknowns -// u ∈ H(curl) - -// Fully discretised equations -// (α∇×u, ∇×u') - (s0, u') - <(α∇×u) × n, u'> = 0 - -// Rewritten as -// a1(u, u') = b1(u') - -// where -// a1(u, u') = (α∇×u, ∇×u') -// b1(u') = (s0, u') + <(α∇×u) × n, u'> -#include "statics_formulation.h" - -#include - -namespace platypus -{ - -StaticsFormulation::StaticsFormulation(std::string alpha_coef_name, std::string h_curl_var_name) - : _alpha_coef_name(std::move(alpha_coef_name)), _h_curl_var_name(std::move(h_curl_var_name)) -{ -} - -void -StaticsFormulation::ConstructJacobianPreconditioner() -{ - std::shared_ptr precond{std::make_shared( - GetProblem()->_gridfunctions.Get(_h_curl_var_name)->ParFESpace())}; - - precond->SetSingularProblem(); - precond->SetPrintLevel(-1); - - GetProblem()->_jacobian_preconditioner = precond; -} - -void -StaticsFormulation::ConstructJacobianSolver() -{ - ConstructJacobianSolverWithOptions(SolverType::HYPRE_FGMRES, - {._max_iteration = 100, ._k_dim = 10}); -} - -void -StaticsFormulation::ConstructOperator() -{ - auto new_operator = std::make_unique( - *GetProblem(), _h_curl_var_name, _alpha_coef_name); - - GetProblem()->SetOperator(std::move(new_operator)); -} - -void -StaticsFormulation::RegisterGridFunctions() -{ - int & myid = GetProblem()->_myid; - platypus::GridFunctions & gridfunctions = GetProblem()->_gridfunctions; - platypus::FESpaces & fespaces = GetProblem()->_fespaces; - - // Register default ParGridFunctions of state gridfunctions if not provided - if (!gridfunctions.Has(_h_curl_var_name)) - { - if (myid == 0) - { - MFEM_WARNING(_h_curl_var_name << " not found in gridfunctions: building " - "gridfunction from defaults"); - } - AddFESpace(std::string("_HCurlFESpace"), std::string("ND_3D_P2")); - AddGridFunction(_h_curl_var_name, std::string("_HCurlFESpace")); - }; -}; - -void -StaticsFormulation::RegisterCoefficients() -{ - platypus::Coefficients & coefficients = GetProblem()->_coefficients; - if (!coefficients._scalars.Has(_alpha_coef_name)) - { - MFEM_ABORT(_alpha_coef_name + " coefficient not found."); - } -} - -StaticsOperator::StaticsOperator(platypus::Problem & problem, - std::string h_curl_var_name, - std::string stiffness_coef_name) - : ProblemOperator(problem), - _h_curl_var_name(std::move(h_curl_var_name)), - _stiffness_coef_name(std::move(stiffness_coef_name)) -{ -} - -void -StaticsOperator::SetGridFunctions() -{ - _trial_var_names.push_back(_h_curl_var_name); - ProblemOperator::SetGridFunctions(); -}; - -void -StaticsOperator::Init(mfem::Vector & X) -{ - ProblemOperator::Init(X); - _stiff_coef = _problem._coefficients._scalars.Get(_stiffness_coef_name); -} - -/* -This is the main method that solves for u. - -Unknowns -s0 ∈ H(div) divergence-free source field -u ∈ H(curl) - -Fully discretised equations -(α∇×u, ∇×u') - (s0, u') - <(α∇×u) × n, u'> = 0 -*/ -void -StaticsOperator::Solve(mfem::Vector & X) -{ - mfem::ParGridFunction & gf(*_trial_variables.at(0)); - gf = 0.0; - mfem::ParLinearForm lf(gf.ParFESpace()); - lf = 0.0; - mfem::Array ess_bdr_tdofs; - _problem._bc_map.ApplyEssentialBCs(_h_curl_var_name, ess_bdr_tdofs, gf, _problem._pmesh.get()); - _problem._bc_map.ApplyIntegratedBCs(_h_curl_var_name, lf, _problem._pmesh.get()); - lf.Assemble(); - _problem._sources.Apply(&lf); - mfem::ParBilinearForm blf(gf.ParFESpace()); - blf.AddDomainIntegrator(new mfem::CurlCurlIntegrator(*_stiff_coef)); - blf.Assemble(); - blf.Finalize(); - mfem::HypreParMatrix curl_mu_inv_curl; - mfem::HypreParVector sol_tdofs(gf.ParFESpace()); - mfem::HypreParVector rhs_tdofs(gf.ParFESpace()); - blf.FormLinearSystem(ess_bdr_tdofs, gf, lf, curl_mu_inv_curl, sol_tdofs, rhs_tdofs); - - // Define and apply a parallel FGMRES solver for AX=B with the AMS - // preconditioner from hypre. - _problem._jacobian_solver->SetOperator(curl_mu_inv_curl); - _problem._jacobian_solver->Mult(rhs_tdofs, sol_tdofs); - blf.RecoverFEMSolution(sol_tdofs, lf, gf); -} - -} // namespace platypus diff --git a/src/formulations/steady_state_em_formulation.C b/src/formulations/steady_state_em_formulation.C deleted file mode 100644 index 80c241a4..00000000 --- a/src/formulations/steady_state_em_formulation.C +++ /dev/null @@ -1,8 +0,0 @@ -#include "steady_state_em_formulation.h" - -namespace platypus -{ - -SteadyStateEMFormulation::SteadyStateEMFormulation() = default; - -} // namespace platypus diff --git a/src/formulations/time_domain_em_formulation.C b/src/formulations/time_domain_em_formulation.C deleted file mode 100644 index 726ec2cb..00000000 --- a/src/formulations/time_domain_em_formulation.C +++ /dev/null @@ -1,58 +0,0 @@ -#include "time_domain_em_formulation.h" - -namespace platypus -{ - -// enum EMField { -// ElectricField, -// MagneticField, -// MagneticFluxDensity, -// CurrentDensity, -// LorentzForceDensity, -// JouleHeatingPower -// }; - -// enum EMPotential { -// MagneticVectorPotential, -// ElectricScalarPotential, -// ElectricVectorPotential, -// MagneticScalarPotential -// }; - -// enum EMMaterialProperty { -// ElectricConductivity, -// ElectricResistivity, -// MagneticPermeability, -// MagneticReluctivity, -// ElectricPermittivity -// }; - -// // req_coef = blah -// // check - -// Struct: motive: to automate AddAuxSolver -// and more easily validate fields - -// Formulation(EMFieldSet, EMMaterialProperties) -// struct EMFieldSet: -// active_fields = std::map - -// std::string enum_to_longname(EMField field) { -// switch (field) { -// case ElectricField: -// return "electric_field"; -// case MagneticField: -// return "magnetic_field"; -// case MagneticFluxDensity: -// return "magnetic_flux_density"; -// case CurrentDensity: -// return "current_density"; -// default: -// return ""; -// } -// } - -TimeDomainEMFormulation::TimeDomainEMFormulation() = default; -; - -} // namespace platypus