added changes in pde_solver.* from idris-etal-aviation-2024

src/physics/pde_solver.cpp

@@ -1,3 +1,4 @@
+#include <cstddef>
 #include "finite_element_dual.hpp"
 #include "mfem.hpp"
 
@@ -149,7 +150,7 @@ MISOMesh::~MISOMesh()
    /// If we started PCU and we're the last one using it, close it
    if (!PCU_previously_initialized && pumi_mesh_count == 0)
    {
-#ifdef HAVE_EGADS
+#ifdef MFEM_USE_EGADS
       gmi_egads_stop();
 #endif
 #ifdef HAVE_SIMMETRIX
@@ -204,12 +205,17 @@ MISOMesh constructMesh(MPI_Comm comm,
    {
       throw MISOException("Unrecognized mesh file extension!\n");
    }
+   // auto *nodes = mesh.mesh->GetNodes();
+   // if (nodes == nullptr)
+   // {
+   //    mesh.mesh->SetCurvature(1, false, 3, mfem::Ordering::byVDIM);
+   // }
    mesh.mesh->EnsureNodes();
 
-   if (!keep_boundaries)
-   {
-      mesh.mesh->RemoveInternalBoundaries();
-   }
+   // if (!keep_boundaries)
+   // {
+   //    mesh.mesh->RemoveInternalBoundaries();
+   // }
 
    // std::ofstream file("pde_solver_mesh.mesh");
    // mesh.mesh->Print(file);
@@ -339,37 +345,76 @@ PDESolver::PDESolver(MPI_Comm incomm,
    mesh_(constructMesh(comm, options["mesh"], std::move(smesh))),
    materials(material_library)
 {
+   /// loop over all components specified in options and add their specified
+   /// materials to the solver's known material library
+   if (solver_options.contains("components"))
+   {
+      for (const auto &component : solver_options["components"])
+      {
+         const auto &material = component["material"];
+         if (material.is_string())
+         {
+            continue;
+         }
+         else
+         {
+            const auto &material_name = material["name"].get<std::string>();
+            materials[material_name].merge_patch(material);
+         }
+      }
+   }
+
    fields.emplace(
        "state",
        FiniteElementState(mesh(), options["space-dis"], num_states, "state"));
+   fields.emplace(
+       "dirichlet_bc",
+       FiniteElementState(
+           mesh(), options["space-dis"], num_states, "dirichlet_bc"));
    fields.emplace(
        "adjoint",
        FiniteElementState(mesh(), options["space-dis"], num_states, "adjoint"));
    duals.emplace(
        "residual",
        FiniteElementDual(mesh(), options["space-dis"], num_states, "residual"));
 
+   // // The pm demag constraint field
+   // fields.emplace(
+   //     "pm_demag_field",
+   //     FiniteElementState(mesh(), options["space-dis"], num_states,
+   //     "pm_demag_field"));
+
    setUpExternalFields();
 }
 
-PDESolver::PDESolver(MPI_Comm incomm,
-                     const nlohmann::json &solver_options,
-                     std::function<int(const nlohmann::json &, int)> num_states,
-                     std::unique_ptr<mfem::Mesh> smesh)
+PDESolver::PDESolver(
+    MPI_Comm incomm,
+    const nlohmann::json &solver_options,
+    const std::function<int(const nlohmann::json &, int)> &num_states,
+    std::unique_ptr<mfem::Mesh> smesh)
  : AbstractSolver2(incomm, solver_options),
    mesh_(constructMesh(comm, options["mesh"], std::move(smesh))),
    materials(material_library)
 {
    int ns = num_states(solver_options, mesh().SpaceDimension());
    fields.emplace(
        "state", FiniteElementState(mesh(), options["space-dis"], ns, "state"));
+   fields.emplace(
+       "dirichlet_bc",
+       FiniteElementState(mesh(), options["space-dis"], ns, "dirichlet_bc"));
    fields.emplace(
        "adjoint",
        FiniteElementState(mesh(), options["space-dis"], ns, "adjoint"));
    duals.emplace(
        "residual",
        FiniteElementDual(mesh(), options["space-dis"], ns, "residual"));
 
+   // // The pm demag constraint field
+   // fields.emplace(
+   //     "pm_demag_field",
+   //     FiniteElementState(mesh(), options["space-dis"], ns,
+   //     "pm_demag_field"));
+
    setUpExternalFields();
 }
 
@@ -427,6 +472,11 @@ void PDESolver::setState_(std::any function,
        { fields.at(name).project(vec_fun, state); },
        [&](mfem::VectorCoefficient *vec_coeff)
        { fields.at(name).project(*vec_coeff, state); });
+
+   if (name == "state")
+   {
+      PDESolver::setState_(function, "dirichlet_bc", state);
+   }
 }
 
 double PDESolver::calcStateError_(std::any ex_sol,
@@ -473,6 +523,8 @@ void PDESolver::initialHook(const mfem::Vector &state)
    int inverted_elems = mesh().CheckElementOrientation(false);
    if (inverted_elems > 0)
    {
+      mesh().PrintVTU("inverted_mesh", mfem::VTKFormat::BINARY, true);
+      std::cout << "!!!!inverted mesh!!!!\n";
       throw MISOException("Mesh contains inverted elements!\n");
    }
    else
@@ -498,6 +550,7 @@ void PDESolver::terminalHook(int iter,
                              const mfem::Vector &state)
 {
    AbstractSolver2::terminalHook(iter, t_final, state);
+   derivedPDETerminalHook(iter, t_final, state);
 }
 
 }  // namespace miso

src/physics/pde_solver.hpp

@@ -47,8 +47,8 @@ struct MISOMesh
    MISOMesh(const MISOMesh &) = delete;
    MISOMesh &operator=(const MISOMesh &) = delete;
 
-   MISOMesh(MISOMesh &&) noexcept;
-   MISOMesh &operator=(MISOMesh &&) noexcept;
+   MISOMesh(MISOMesh && /*other*/) noexcept;
+   MISOMesh &operator=(MISOMesh && /*other*/) noexcept;
 
    ~MISOMesh();
 #endif
@@ -95,16 +95,16 @@ class PDESolver : public AbstractSolver2
    /// available to determine the number of states.
    PDESolver(MPI_Comm incomm,
              const nlohmann::json &solver_options,
-             std::function<int(const nlohmann::json &, int)> num_states,
+             const std::function<int(const nlohmann::json &, int)> &num_states,
              std::unique_ptr<mfem::Mesh> smesh = nullptr);
 
 protected:
    /// object defining the mfem computational mesh
    MISOMesh mesh_;
 
    /// Reference to solver state vector
-   mfem::ParMesh &mesh() { return *mesh_.mesh; }
-   const mfem::ParMesh &mesh() const { return *mesh_.mesh; }
+   mfem::ParMesh &mesh() const { return *mesh_.mesh; }
+   // const mfem::ParMesh &mesh() const { return *mesh_.mesh; }
 
    /// solver material properties
    nlohmann::json materials;
@@ -125,6 +125,11 @@ class PDESolver : public AbstractSolver2
    FiniteElementDual &res_vec() { return duals.at("residual"); }
    const FiniteElementDual &res_vec() const { return duals.at("residual"); }
 
+   /// Reference to the pm demag constraint field
+   // FiniteElementState &pm_demag_field() { return fields.at("pm_demag_field");
+   // } const FiniteElementState &pm_demag_field() const { return
+   // fields.at("pm_demag_field"); }
+
    /// Reference to the state vectors finite element space
    mfem::ParFiniteElementSpace &fes() { return state().space(); }
    const mfem::ParFiniteElementSpace &fes() const { return state().space(); }
@@ -146,7 +151,7 @@ class PDESolver : public AbstractSolver2
    /// client overwrites this definition; however, there is a call to the
    /// virtual function derivedPDEinitialHook(state) that the client can
    /// overwrite.
-   virtual void initialHook(const mfem::Vector &state) override final;
+   void initialHook(const mfem::Vector &state) final;
 
    /// Code in a derived class that should be executed before time-stepping
    /// \param[in] state - the current state
@@ -159,10 +164,10 @@ class PDESolver : public AbstractSolver2
    /// \param[in] state - the current state
    /// \note This is `final` because we want to ensure that
    /// AbstractSolver2::iterationHook() is called.
-   virtual void iterationHook(int iter,
-                              double t,
-                              double dt,
-                              const mfem::Vector &state) override final;
+   void iterationHook(int iter,
+                      double t,
+                      double dt,
+                      const mfem::Vector &state) final;
 
    /// Code in a derived class that should be executed each time step
    /// \param[in] iter - the current iteration
@@ -179,9 +184,7 @@ class PDESolver : public AbstractSolver2
    /// \param[in] iter - the terminal iteration
    /// \param[in] t_final - the final time
    /// \param[in] state - the current state
-   virtual void terminalHook(int iter,
-                             double t_final,
-                             const mfem::Vector &state) override final;
+   void terminalHook(int iter, double t_final, const mfem::Vector &state) final;
 
    /// Code in a derived class that should be executed after time stepping ends
    /// \param[in] iter - the terminal iteration

test/regression/CMakeLists.txt

@@ -14,16 +14,22 @@ endfunction(create_tests)
 
 # group together all regression tests
 set(REGRESSION_TEST_SRCS
-   test_steady_vortex
    test_navier_stokes_mms
    test_navier_stokes_3Dmms
-   test_inviscid_mms
-   test_FarField3D_exact
    test_thermal_cube
+   test_steady_thermal_cube
    test_miso_inputs
    test_magnetostatic_box
+   test_magnetostatic_box2d
+   #test_meshmovement_box
+   #test_meshmovement_annulus
    test_ac_loss
    test_magnetostatic_box_new
+   # test_2d_magnet_in_box
+   test_magnetostatic_residual
+   test_thermal_residual
+   test_weak_boundary
+   test_thermal_solver
 )
 
 create_tests("${REGRESSION_TEST_SRCS}" regression_data.cpp)