Got the optimization running, with some other additions

src/common/data_logging.hpp

@@ -120,10 +120,13 @@ class ParaViewLogger
                   double time,
                   int rank)
    {
-      fields.at(fieldname)->SetFromTrueDofs(state);
-      pv.SetCycle(timestep);
-      pv.SetTime(time);
-      pv.Save();
+      if (fields.count(fieldname) > 0)
+      {
+         fields.at(fieldname)->SetFromTrueDofs(state);
+         pv.SetCycle(timestep);
+         pv.SetTime(time);
+         pv.Save();
+      }
    }
 
    void registerField(const std::string &name, mfem::ParGridFunction &field)

src/physics/CMakeLists.txt

@@ -11,8 +11,10 @@ set(MISO_PHYSICS_HEADERS
    mfem_common_integ.hpp
    pde_solver.hpp
    physics.hpp
+   solver.hpp
 )
 
+
 target_sources(miso
    PRIVATE
       common_outputs.cpp
@@ -24,13 +26,14 @@ target_sources(miso
       miso_nonlinearform.cpp
       mfem_common_integ.cpp
       pde_solver.cpp
+      solver.cpp
       ${MISO_PHYSICS_HEADERS}
 )
 
-# get_target_property(MACH_PUBLIC_HEADERS mach PUBLIC_HEADER)
-# set(MACH_PUBLIC_HEADERS ${MACH_PUBLIC_HEADERS} ${MACH_PHYSICS_HEADERS})
-# set_target_properties(mach PROPERTIES
-#    PUBLIC_HEADER "${MACH_PUBLIC_HEADERS}"
+# get_target_property(MISO_PUBLIC_HEADERS miso PUBLIC_HEADER)
+# set(MISO_PUBLIC_HEADERS ${MISO_PUBLIC_HEADERS} ${MISO_PHYSICS_HEADERS})
+# set_target_properties(miso PROPERTIES
+#    PUBLIC_HEADER "${MISO_PUBLIC_HEADERS}"
 # )
 
 add_public_headers_to_target(miso "${MISO_PHYSICS_HEADERS}")
@@ -43,4 +46,6 @@ target_include_directories(miso
 
 add_subdirectory(electromagnetics)
 add_subdirectory(fluidflow)
+# add_subdirectory(joule)
+add_subdirectory(meshmove)
 add_subdirectory(thermal)

src/physics/electromagnetics/CMakeLists.txt

@@ -38,6 +38,12 @@ target_sources(miso
       ${MISO_PHYSICS_ELECTROMAGNETICS_HEADERS}
 )
 
+# get_target_property(MISO_PUBLIC_HEADERS miso PUBLIC_HEADER)
+# set(MISO_PUBLIC_HEADERS ${MISO_PUBLIC_HEADERS} ${MISO_PHYSICS_ELECTROMAGNETICS_HEADERS})
+# set_target_properties(miso PROPERTIES
+#    PUBLIC_HEADER "${MISO_PUBLIC_HEADERS}"
+# )
+
 add_public_headers_to_target(miso "${MISO_PHYSICS_ELECTROMAGNETICS_HEADERS}")
 
 target_include_directories(miso

src/physics/electromagnetics/electromag_integ.cpp

@@ -404,7 +404,7 @@ void NonlinearDiffusionIntegratorMeshRevSens::AssembleRHSElementVect(
 
       /// const double pointflux_norm = pointflux.Norml2();
       pointflux_bar = 0.0;
-      if (fabs(pointflux_norm) > 1e-14)
+      if (fabs(pointflux_norm) > 1e-16)
       {
          add(pointflux_bar,
              pointflux_norm_bar / pointflux_norm,
@@ -1651,14 +1651,14 @@ void DGInteriorFaceDiffusionIntegratorMeshRevSens::AssembleRHSElementVect(
 
       /// const double ip_flux2_norm = ip_flux2.Norml2();
       ip_flux2_bar = 0.0;
-      if (fabs(ip_flux2_norm) > 1e-14)
+      if (fabs(ip_flux2_norm) > 1e-16)
       {
          ip_flux2_bar.Add(ip_flux2_norm_bar / ip_flux2_norm, ip_flux2);
       }
 
       /// const double ip_flux1_norm = ip_flux1.Norml2();
       ip_flux1_bar = 0.0;
-      if (fabs(ip_flux1_norm) > 1e-14)
+      if (fabs(ip_flux1_norm) > 1e-16)
       {
          ip_flux1_bar.Add(ip_flux1_norm_bar / ip_flux1_norm, ip_flux1);
       }

src/physics/joule/CMakeLists.txt

@@ -0,0 +1,10 @@
+target_sources(miso
+   PRIVATE
+      joule.cpp
+      joule.hpp
+)
+
+target_include_directories(miso
+   PUBLIC
+      "${CMAKE_CURRENT_SOURCE_DIR}"
+)

src/physics/joule/joule.cpp

@@ -0,0 +1,151 @@
+#include <memory>
+
+#include "joule.hpp"
+#include "magnetostatic.hpp"
+#include "thermal.hpp"
+
+using namespace std;
+using namespace mfem;
+
+namespace miso
+{
+/// TODO: read options file and construct options files for EM and
+JouleSolver::JouleSolver(const std::string &opt_file_name,
+                         std::unique_ptr<mfem::Mesh> smesh,
+                         MPI_Comm comm)
+ // : AbstractSolver(opt_file_name, move(smesh))
+ : AbstractSolver(opt_file_name, comm)
+{
+   nlohmann::json em_opts = options["em-opts"];
+   nlohmann::json thermal_opts = options["thermal-opts"];
+
+   auto mesh_file = options["mesh"]["file"].get<std::string>();
+   auto model_file = options["mesh"]["model-file"].get<std::string>();
+   auto mesh_out_file = options["mesh"]["out-file"].get<std::string>();
+
+   /// get mesh file name and extension
+   std::string mesh_name;
+   std::string mesh_ext;
+   {
+      size_t i = mesh_file.rfind('.', mesh_file.length());
+      if (i != string::npos)
+      {
+         mesh_name = (mesh_file.substr(0, i));
+         mesh_ext = (mesh_file.substr(i + 1, mesh_file.length() - i));
+      }
+      else
+      {
+         throw MISOException(
+             "JouleSolver::JouleSolver()\n"
+             "\tMesh file has no extension!\n");
+      }
+   }
+
+   em_opts["mesh"]["file"] = mesh_name + "_em." + mesh_ext;
+   thermal_opts["mesh"]["file"] = mesh_name + "_thermal." + mesh_ext;
+
+   /// get model file name and extension
+   std::string model_name;
+   std::string model_ext;
+   {
+      size_t i = model_file.rfind('.', model_file.length());
+      if (i != string::npos)
+      {
+         model_name = (model_file.substr(0, i));
+         model_ext = (model_file.substr(i + 1, model_file.length() - i));
+      }
+      else
+      {
+         throw MISOException(
+             "JouleSolver::JouleSolver()\n"
+             "\tModel file has no extension!\n");
+      }
+   }
+
+   em_opts["mesh"]["model-file"] = model_name + "_em." + model_ext;
+   thermal_opts["mesh"]["model-file"] = model_name + "_thermal." + model_ext;
+
+   em_opts["mesh"]["out-file"] = mesh_out_file + "_em";
+   thermal_opts["mesh"]["out-file"] = mesh_out_file + "_thermal";
+
+   em_opts["components"] = options["components"];
+   thermal_opts["components"] = options["components"];
+
+   em_opts["problem-opts"] = options["problem-opts"];
+   thermal_opts["problem-opts"] = options["problem-opts"];
+
+   /// TODO: need to do this until Magnetostatic solver is updated to support
+   /// newer abstract solver construction model
+   std::string em_opt_filename = "em_opt_file.json";
+   {
+      std::ofstream em_opt_file(em_opt_filename, std::ios::trunc);
+      em_opt_file << em_opts;
+      em_opt_file.close();
+   }
+
+   *out << "EM options:\n";
+   *out << setw(3) << em_opts << endl;
+
+   em_solver =
+       std::make_unique<MagnetostaticSolver>(em_opt_filename, nullptr, comm);
+   /// TODO: this should be moved to an init derived when a factory is made
+   // em_solver->initDerived();
+
+   thermal_solver =
+       std::make_unique<ThermalSolver>(thermal_opts, nullptr, comm);
+   // thermal_solver->initDerived();
+}
+
+void JouleSolver::initDerived()
+{
+   em_solver->initDerived();
+   em_fields = em_solver->getFields();
+   thermal_solver->setAField(em_fields[0]);
+}
+
+/// TODO: Change this in AbstractSolver to mark a flag so that unsteady
+/// solutions can be saved
+void JouleSolver::printSolution(const std::string &filename, int refine)
+{
+   std::string em_filename = filename + "_em";
+   std::string thermal_filename = filename + "_thermal";
+   em_solver->printSolution(thermal_filename, refine);
+   thermal_solver->printSolution(thermal_filename, refine);
+}
+
+void JouleSolver::setInitialCondition(
+    const std::function<double(const mfem::Vector &)> &u_init)
+{
+   thermal_init = u_init;
+}
+
+std::vector<GridFunType *> JouleSolver::getFields()
+{
+   return {thermal_fields[0], em_fields[0], em_fields[1]};
+}
+
+void JouleSolver::solveForState()
+{
+   em_solver->solveForState();
+   thermal_solver->initDerived();
+   thermal_fields = thermal_solver->getFields();
+   thermal_solver->setInitialCondition(thermal_init);
+   thermal_solver->solveForState();
+}
+
+void JouleSolver::solveForAdjoint(const std::string &fun)
+{
+   if (fun == "temp-agg")
+   {
+      thermal_solver->solveForAdjoint(fun);
+      Vector psiTdRtdA;
+      thermal_solver->getASensitivity(psiTdRtdA);
+      em_solver->solveForAdjoint(psiTdRtdA);
+   }
+   else if (fun == "co-energy")
+   {
+      em_solver->solveForAdjoint(fun);
+   }
+}
+
+}  // namespace miso

src/physics/joule/joule.hpp

@@ -0,0 +1,68 @@
+#ifndef MISO_JOULE
+#define MISO_JOULE
+
+#include "solver.hpp"
+
+namespace miso
+{
+class MagnetostaticSolver;
+class ThermalSolver;
+
+class JouleSolver : public AbstractSolver
+{
+public:
+   /// Class constructor.
+   /// \param[in] opt_file_name - file where options are stored
+   /// \param[in] smesh - if provided, defines the mesh for the problem
+   /// \param[in] comm - MPI communicator for parallel operations
+   JouleSolver(const std::string &opt_file_name,
+               std::unique_ptr<mfem::Mesh> smesh,
+               MPI_Comm comm);
+
+   /// Fully initialize the Joule Solver and its sub-solvers
+   void initDerived() override;
+
+   /// Write the solutions of both the EM and thermal problems to a vtk file
+   /// \param[in] file_name - prefix file name **without** .vtk extension
+   ///                        EM vtk file will be filename_em,
+   ///                        Thermal vtk file will be filename_thermal
+   /// \param[in] refine - if >=0, indicates the number of refinements to make
+   /// \todo make this work for parallel!
+   /// \note the `refine` argument is useful for high-order meshes and
+   /// solutions; it divides the elements up so it is possible to visualize.
+   void printSolution(const std::string &filename, int refine = -1) override;
+
+   /// Initializes the state variable to a given function.
+   /// \param[in] u_init - function that defines the initial condition
+   void setInitialCondition(
+       const std::function<double(const mfem::Vector &)> &u_init) override;
+
+   /// \brief Returns a vector of pointers to grid functions that define fields
+   /// returns {T, A, B}
+   std::vector<GridFunType *> getFields() override;
+
+   void solveForState() override;
+
+   void solveForAdjoint(const std::string &fun) override;
+
+   // mfem::Vector* getMeshSensitivities() override;
+
+   int getNumState() override { return 0; }
+
+private:
+   std::unique_ptr<MagnetostaticSolver> em_solver;
+   std::unique_ptr<ThermalSolver> thermal_solver;
+
+   std::vector<GridFunType *> em_fields;
+   std::vector<GridFunType *> thermal_fields;
+
+   std::function<double(const mfem::Vector &)> thermal_init;
+
+   friend SolverPtr createSolver<JouleSolver>(const std::string &opt_file_name,
+                                              std::unique_ptr<mfem::Mesh> smesh,
+                                              MPI_Comm comm);
+};
+
+}  // namespace miso
+
+#endif

src/physics/meshmove/CMakeLists.txt

@@ -0,0 +1,26 @@
+set(MISO_PHYSICS_MESHMOVE_HEADERS
+   mesh_move_integ.hpp
+   mesh_movement.hpp
+)
+
+target_sources(miso
+   PRIVATE
+      mesh_move_integ.cpp
+      mesh_movement.cpp
+      ${MISO_PHYSICS_MESHMOVE_HEADERS}
+)
+
+# get_target_property(MISO_PUBLIC_HEADERS miso PUBLIC_HEADER)
+# set(MISO_PUBLIC_HEADERS ${MISO_PUBLIC_HEADERS} ${MISO_PHYSICS_MESHMOVE_HEADERS})
+# set_target_properties(miso PROPERTIES
+#    PUBLIC_HEADER "${MISO_PUBLIC_HEADERS}"
+# )
+
+add_public_headers_to_target(miso "${MISO_PHYSICS_MESHMOVE_HEADERS}")
+
+target_include_directories(miso
+   PUBLIC
+     "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>"
+     "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDE_DIR}>"
+)
+