Merge pull request #68 from OptimalDesignLab/miso-merge

Miso merge

miso/mesh_warper.py

@@ -1,16 +1,16 @@
 import numpy as np
 import openmdao.api as om
 
-from .pyMach import MeshWarper
+from .pyMISO import MeshWarper
 
-class MachMeshWarper(om.ImplicitComponent):
+class MISOMeshWarper(om.ImplicitComponent):
     def initialize(self):
         self.options.declare("warper", types=MeshWarper, recordable=False)
 
     def setup(self):
         warper = self.options["warper"]
 
-        # hold map of vector-valued I/O names -> contiguous vectors to pass to Mach
+        # hold map of vector-valued I/O names -> contiguous vectors to pass to MISO
         self.vectors = dict()
 
         local_surf_mesh_size = warper.getSurfaceCoordsSize()

sandbox/CMakeLists.txt

@@ -14,7 +14,7 @@ function(create_sandbox source_list)
 endfunction(create_sandbox)
 
 set(SANDBOX_SRCS
-   air-gap
+   # air-gap
    steady_vortex
    steady_vortex_adjoint
    unsteady_vortex

src/physics/thermal/thermal_residual.cpp

@@ -1,21 +1,26 @@
+#include <iostream>
+#include <memory>
 #include <string>
 
 #include "mfem.hpp"
 #include "nlohmann/json.hpp"
 
+#include "electromag_integ.hpp"
 #include "coefficient.hpp"
 #include "miso_input.hpp"
-#include "mfem_common_integ.hpp"
 #include "thermal_integ.hpp"
-
 #include "thermal_residual.hpp"
+#include "utils.hpp"
 
 namespace miso
 {
 int getSize(const ThermalResidual &residual) { return getSize(residual.res); }
 
 void setInputs(ThermalResidual &residual, const miso::MISOInputs &inputs)
 {
+   residual.kappa->setInputs(inputs);
+   residual.rho->setInputs(inputs);
+
    setInputs(residual.res, inputs);
 
    setVectorFromInputs(inputs, "thermal_load", residual.load);
@@ -70,6 +75,15 @@ void setUpAdjointSystem(ThermalResidual &residual,
    setUpAdjointSystem(residual.res, adj_solver, inputs, state_bar, adjoint);
 }
 
+void finalizeAdjointSystem(ThermalResidual &residual,
+                           mfem::Solver &adj_solver,
+                           const MISOInputs &inputs,
+                           mfem::Vector &state_bar,
+                           mfem::Vector &adjoint)
+{
+   finalizeAdjointSystem(residual.res, adj_solver, inputs, state_bar, adjoint);
+}
+
 double jacobianVectorProduct(ThermalResidual &residual,
                              const mfem::Vector &wrt_dot,
                              const std::string &wrt)
@@ -83,6 +97,26 @@ void jacobianVectorProduct(ThermalResidual &residual,
                            const std::string &wrt,
                            mfem::Vector &res_dot)
 {
+   // if wrt starts with prefix "thermal_load"
+   if (wrt.rfind("thermal_load", 0) == 0)
+   {
+      /// TODO: Test this implementation
+      // Thermal load is subtracted from the (stiffness matrix * thermal) state
+      // term Therefore derivative of residual w/r/t thermal_load is negative
+      // identity matrix (represented as negative of input vector in matrix-free
+      // form)
+      res_dot = 0.0;
+      res_dot.Add(-1.0, wrt_dot);
+      return;
+   }
+   // if wrt starts with prefix "temperature"
+   /// TODO: Determine if this should be "state" instead of temperature
+   else if (wrt.rfind("temperature", 0) == 0)
+   {
+      /// NOTE: Derivative of the thermal residual wrt the thermal state is
+      /// Jacobian (already implemented above)
+      return;
+   }
    jacobianVectorProduct(residual.res, wrt_dot, wrt, res_dot);
 }
 
@@ -99,6 +133,26 @@ void vectorJacobianProduct(ThermalResidual &residual,
                            const std::string &wrt,
                            mfem::Vector &wrt_bar)
 {
+   // if wrt starts with prefix "thermal_load"
+   if (wrt.rfind("thermal_load", 0) == 0)
+   {
+      /// TODO: Test this implementation
+      // Thermal load is subtracted from the (stiffness matrix * thermal) state
+      // term Therefore derivative of residual w/r/t thermal_load is negative
+      // identity matrix (represented as negative of input vector in matrix-free
+      // form)
+      wrt_bar = 0.0;
+      wrt_bar.Add(-1.0, res_bar);
+      return;
+   }
+   // if wrt starts with prefix "temperature"
+   /// TODO: Determine if this should be "state" instead of temperature
+   else if (wrt.rfind("temperature", 0) == 0)
+   {
+      /// NOTE: Derivative of the thermal residual wrt the thermal state is
+      /// Jacobian (already implemented above)
+      return;
+   }
    vectorJacobianProduct(residual.res, res_bar, wrt, wrt_bar);
 }
 
@@ -113,35 +167,205 @@ ThermalResidual::ThermalResidual(
     const nlohmann::json &options,
     const nlohmann::json &materials)
  : res(fes, fields),
+   g(std::make_unique<MeshDependentCoefficient>(
+       std::make_unique<mfem::GridFunctionCoefficient>(
+           &fields.at("dirichlet_bc").gridFunc()))),
    kappa(
        constructMaterialCoefficient("kappa", options["components"], materials)),
    rho(constructMaterialCoefficient("rho", options["components"], materials)),
-   cv(constructMaterialCoefficient("cv", options["components"], materials)),
+   // cv(constructMaterialCoefficient("cv", options["components"], materials)),
    prec(constructPreconditioner(fes, options["lin-prec"]))
 
 {
+   res.addDomainIntegrator(new NonlinearDiffusionIntegrator(*kappa));
+
    const auto &basis_type =
        options["space-dis"]["basis-type"].get<std::string>();
-   if (basis_type != "H1" && basis_type != "h1" && basis_type != "CG" &&
-       basis_type != "cg")
+
+   if (basis_type == "L2" || basis_type == "l2" || basis_type == "DG" ||
+       basis_type == "dg")
    {
-      throw MISOException(
-          "Thermal residual currently only supports H1 state field!\n");
+      auto mu = options["space-dis"].value("sipg-penalty", -1.0);
+      if (mu < 0)
+      {
+         auto degree = options["space-dis"]["degree"].get<double>();
+         mu = pow(degree + 1, 2);
+      }
+      res.addInteriorFaceIntegrator(
+          new DGInteriorFaceDiffusionIntegrator(*kappa, mu));
+      std::cout << "adding sipg integ!\n";
    }
 
-   res.addDomainIntegrator(new mfem::DiffusionIntegrator(*kappa));
-
    if (options.contains("bcs"))
    {
-      auto &bcs = options["bcs"];
+      const auto &bcs = options["bcs"];
 
       // convection heat transfer boundary condition
       if (bcs.contains("convection"))
       {
-         const auto &bdr_attr_marker =
-             bcs["convection"].get<std::vector<int>>();
+         auto &bc = bcs["convection"];
+         std::vector<int> bdr_attr_marker;
+         double h = 1.0;
+         double theta_f = 1.0;
+         if (bc.is_array())
+         {
+            bdr_attr_marker = bc.get<std::vector<int>>();
+         }
+         else if (bc.is_object())
+         {
+            bdr_attr_marker = bc["attrs"].get<std::vector<int>>();
+            h = bc.value("h", h);
+            theta_f = bc.value("fluid_temp", theta_f);
+         }
+
+         res.addBdrFaceIntegrator(new ConvectionBCIntegrator(h, theta_f),
+                                  bdr_attr_marker);
+      }
+
+      if (bcs.contains("outflux"))
+      {
+         auto &bc = bcs["outflux"];
+         std::vector<int> bdr_attr_marker;
+         double flux = 1.0;
+         if (bc.is_array())
+         {
+            bdr_attr_marker = bc.get<std::vector<int>>();
+         }
+         else if (bc.is_object())
+         {
+            bdr_attr_marker = bc["attrs"].get<std::vector<int>>();
+            flux = bc.value("flux", flux);
+         }
+
+         res.addBdrFaceIntegrator(new OutfluxBCIntegrator(flux, 1.0),
+                                  bdr_attr_marker);
+      }
+
+      // dirichlet boundary condition
+      if (bcs.contains("essential"))
+      {
+         if (basis_type == "L2" || basis_type == "l2" || basis_type == "DG" ||
+             basis_type == "dg")
+         {
+            std::vector<int> bdr_attr_marker;
+            if (bcs["essential"].is_array())
+            {
+               bdr_attr_marker = bcs["essential"].get<std::vector<int>>();
+            }
+            else
+            {
+               throw MISOException(
+                   "Unrecognized JSON type for boundary attrs!");
+            }
+
+            auto mu = options["space-dis"].value("sipg-penalty", -1.0);
+            if (mu < 0)
+            {
+               auto degree = options["space-dis"]["degree"].get<double>();
+               mu = pow(degree + 1, 2);
+            }
+            std::cout << "mu: " << mu << "\n";
+            res.addBdrFaceIntegrator(
+                new NonlinearDGDiffusionIntegrator(*kappa, *g, mu),
+                bdr_attr_marker);
+         }
+      }
+
+      if (bcs.contains("weak-essential"))
+      {
+         std::vector<int> bdr_attr_marker;
+         if (bcs["weak-essential"].is_array())
+         {
+            bdr_attr_marker = bcs["weak-essential"].get<std::vector<int>>();
+         }
+         else
+         {
+            throw MISOException("Unrecognized JSON type for boundary attrs!");
+         }
+
+         auto mu = options["space-dis"].value("sipg-penalty", -1.0);
+         if (mu < 0)
+         {
+            auto degree = options["space-dis"]["degree"].get<double>();
+            mu = pow(degree + 1, 2);
+         }
+         std::cout << "mu: " << mu << "\n";
+         res.addBdrFaceIntegrator(
+             new NonlinearDGDiffusionIntegrator(*kappa, *g, mu),
+             bdr_attr_marker);
+      }
+   }
+
+   if (options.contains("interfaces"))
+   {
+      const auto &interfaces = options["interfaces"];
+
+      for (const auto &[kind, interface] : interfaces.items())
+      {
+         std::cout << "kind: " << kind << "\n";
+         std::cout << "interface: " << interface << "\n";
+         // thermal contact resistance interface
+         if (kind == "thermal_contact_resistance")
+         {
+            for (const auto &[name, intr] : interface.items())
+            {
+               const auto &attrs = intr["attrs"].get<std::vector<int>>();
+
+               auto mu = options["space-dis"].value("sipg-penalty", -1.0);
+               if (mu < 0)
+               {
+                  auto degree = options["space-dis"]["degree"].get<double>();
+                  mu = pow(degree + 1, 2);
+               }
+
+               res.addInternalBoundaryFaceIntegrator(
+                   new DGInteriorFaceDiffusionIntegrator(*kappa, mu, -1),
+                   attrs);
+
+               const auto h = intr["h"].get<double>();
+
+               auto *integ = new ThermalContactResistanceIntegrator(h, name);
+               // setInputs(*integ, {{"h", h}});
+               res.addInternalBoundaryFaceIntegrator(integ, attrs);
+
+               std::cout << "adding " << name << " TCR integrator!\n";
+               std::cout << "with attrs: " << intr["attrs"] << "\n";
+            }
+         }
+         else if (kind == "convection")
+         {
+            for (const auto &[name, intr] : interface.items())
+            {
+               const auto &attrs = intr["attrs"].get<std::vector<int>>();
+
+               auto mu = options["space-dis"].value("sipg-penalty", -1.0);
+               if (mu < 0)
+               {
+                  auto degree = options["space-dis"]["degree"].get<double>();
+                  mu = pow(degree + 1, 2);
+               }
+
+               res.addInternalBoundaryFaceIntegrator(
+                   new DGInteriorFaceDiffusionIntegrator(*kappa, mu, -1),
+                   attrs);
+
+               const auto h = intr["h"].get<double>();
+               const auto theta_f = intr["fluid_temp"].get<double>();
+
+               auto *integ =
+                   new InternalConvectionInterfaceIntegrator(h, theta_f, name);
+               res.addInternalBoundaryFaceIntegrator(integ, attrs);
+
+               // auto *integ = new ThermalContactResistanceIntegrator(h,
+               // name);
+               // // setInputs(*integ, {{"h", h}});
+               // res.addInternalBoundaryFaceIntegrator(integ, attrs);
 
-         res.addBdrFaceIntegrator(new ConvectionBCIntegrator, bdr_attr_marker);
+               std::cout << "adding " << name
+                         << " internal convection integrator!\n";
+               std::cout << "with attrs: " << intr["attrs"] << "\n";
+            }
+         }
       }
    }
 }