further progress on compliling

Author: jehicken

src/physics/CMakeLists.txt

@@ -11,7 +11,7 @@ set(MISO_PHYSICS_HEADERS
    mfem_common_integ.hpp
    pde_solver.hpp
    physics.hpp
-   solver.hpp
+   # solver.hpp
 )
 
 target_sources(miso
@@ -25,7 +25,7 @@ target_sources(miso
       miso_nonlinearform.cpp
       mfem_common_integ.cpp
       pde_solver.cpp
-      solver.cpp
+      # solver.cpp
       ${MISO_PHYSICS_HEADERS}
 )
 

src/physics/electromagnetics/current_source_functions.cpp

@@ -211,6 +211,24 @@ void box2_current(const xdouble *x, xdouble *J)
    J[2] = 6 * y;
 }
 
+template <typename xdouble = double>
+xdouble box1_current2D(const xdouble *x)
+{
+   auto y = x[1] - .5;
+
+   // J[2] = -current_density*6*y*(1/(M_PI*4e-7)); // for real scaled problem
+   return -6 * y;
+}
+
+template <typename xdouble = double>
+xdouble box2_current2D(const xdouble *x)
+{
+   auto y = x[1] - .5;
+
+   // J[2] = current_density*6*y*(1/(M_PI*4e-7)); // for real scaled problem
+   return 6 * y;
+}
+
 /// function to get the sign of a number
 template <typename T>
 int sgn(T val)
@@ -536,6 +554,66 @@ void box2CurrentSourceRevDiff(adept::Stack &diff_stack,
    source_jac.MultTranspose(V_bar, x_bar);
 }
 
+double box1CurrentSource2D(const mfem::Vector &x)
+{
+   return box1_current2D(x.GetData());
+}
+
+/// TODO: this is how Adept should be used here, eventually
+/// should update the rest of the use cases to match this
+/// (it's more efficient)
+void box1CurrentSource2DRevDiff(adept::Stack &diff_stack,
+                                const mfem::Vector &x,
+                                double J_bar,
+                                mfem::Vector &x_bar)
+{
+   std::array<adouble, 2> x_a;
+   // copy data from mfem::Vector
+   adept::set_values(x_a.data(), x.Size(), x.GetData());
+
+   // start recording
+   diff_stack.new_recording();
+
+   // the depedent variable must be declared after the recording
+   auto J_a = box1_current2D<adouble>(x_a.data());
+
+   J_a.set_gradient(J_bar);
+   diff_stack.compute_adjoint();
+
+   // calculate the vector jacobian product w.r.t position
+   adept::get_gradients(x_a.data(), x.Size(), x_bar.GetData());
+}
+
+double box2CurrentSource2D(const mfem::Vector &x)
+{
+   return box2_current2D(x.GetData());
+}
+
+/// TODO: this is how Adept should be used here, eventually
+/// should update the rest of the use cases to match this
+/// (it's more efficient)
+void box2CurrentSource2DRevDiff(adept::Stack &diff_stack,
+                                const mfem::Vector &x,
+                                double J_bar,
+                                mfem::Vector &x_bar)
+{
+   std::array<adouble, 2> x_a;
+   // copy data from mfem::Vector
+   adept::set_values(x_a.data(), x.Size(), x.GetData());
+
+   // start recording
+   diff_stack.new_recording();
+
+   // the depedent variable must be declared after the recording
+   auto J_a = box2_current2D<adouble>(x_a.data());
+
+   J_a.set_gradient(J_bar);
+   diff_stack.compute_adjoint();
+
+   // calculate the vector jacobian product w.r.t position
+   adept::get_gradients(x_a.data(), x.Size(), x_bar.GetData());
+}
+
 void team13CurrentSource(const mfem::Vector &x, mfem::Vector &J)
 {
    team13_current(x.GetData(), J.GetData());
@@ -566,7 +644,7 @@ void team13CurrentSourceRevDiff(adept::Stack &diff_stack,
 
 }  // anonymous namespace
 
-namespace mach
+namespace miso
 {
 void CurrentDensityCoefficient::cacheCurrentDensity()
 {
@@ -592,7 +670,7 @@ void CurrentDensityCoefficient::resetCurrentDensityFromCache()
    }
 }
 
-bool setInputs(CurrentDensityCoefficient &current, const MachInputs &inputs)
+bool setInputs(CurrentDensityCoefficient &current, const MISOInputs &inputs)
 {
    bool updated = false;
    for (auto &[group, coeff] : current.group_map)
@@ -654,7 +732,7 @@ CurrentDensityCoefficient::CurrentDensityCoefficient(
             auto n_slots = cached_inputs.at("n_slots");
             cached_inputs.emplace("stack_length", 0.345);
             auto stack_length = cached_inputs.at("stack_length");
-            for (auto &attr : attrs)
+            for (const auto &attr : attrs)
             {
                source_coeffs.emplace(
                    attr,
@@ -689,7 +767,7 @@ CurrentDensityCoefficient::CurrentDensityCoefficient(
             auto n_slots = cached_inputs.at("n_slots");
             cached_inputs.emplace("stack_length", 0.345);
             auto stack_length = cached_inputs.at("stack_length");
-            for (auto &attr : attrs)
+            for (const auto &attr : attrs)
             {
                source_coeffs.emplace(
                    attr,
@@ -720,7 +798,7 @@ CurrentDensityCoefficient::CurrentDensityCoefficient(
          }
          else if (source == "x")
          {
-            for (auto &attr : attrs)
+            for (const auto &attr : attrs)
             {
                source_coeffs.emplace(
                    attr,
@@ -743,7 +821,7 @@ CurrentDensityCoefficient::CurrentDensityCoefficient(
          }
          else if (source == "y")
          {
-            for (auto &attr : attrs)
+            for (const auto &attr : attrs)
             {
                source_coeffs.emplace(
                    attr,
@@ -766,7 +844,7 @@ CurrentDensityCoefficient::CurrentDensityCoefficient(
          }
          else if (source == "z")
          {
-            for (auto &attr : attrs)
+            for (const auto &attr : attrs)
             {
                source_coeffs.emplace(
                    attr,
@@ -789,7 +867,7 @@ CurrentDensityCoefficient::CurrentDensityCoefficient(
          }
          else if (source == "-z")
          {
-            for (auto &attr : attrs)
+            for (const auto &attr : attrs)
             {
                source_coeffs.emplace(
                    attr,
@@ -812,7 +890,7 @@ CurrentDensityCoefficient::CurrentDensityCoefficient(
          }
          else if (source == "ring")
          {
-            for (auto &attr : attrs)
+            for (const auto &attr : attrs)
             {
                source_coeffs.emplace(
                    attr,
@@ -834,7 +912,7 @@ CurrentDensityCoefficient::CurrentDensityCoefficient(
          }
          else if (source == "box1")
          {
-            for (auto &attr : attrs)
+            for (const auto &attr : attrs)
             {
                source_coeffs.emplace(
                    attr,
@@ -856,7 +934,7 @@ CurrentDensityCoefficient::CurrentDensityCoefficient(
          }
          else if (source == "box2")
          {
-            for (auto &attr : attrs)
+            for (const auto &attr : attrs)
             {
                source_coeffs.emplace(
                    attr,
@@ -878,7 +956,7 @@ CurrentDensityCoefficient::CurrentDensityCoefficient(
          }
          else if (source == "team13")
          {
-            for (auto &attr : attrs)
+            for (const auto &attr : attrs)
             {
                source_coeffs.emplace(
                    attr,
@@ -903,4 +981,168 @@ CurrentDensityCoefficient::CurrentDensityCoefficient(
    }
 }
 
-}  // namespace mach
+void CurrentDensityCoefficient2D::cacheCurrentDensity()
+{
+   for (auto &[group, coeff] : group_map)
+   {
+      cached_inputs.at(group) = coeff.constant;
+   }
+}
+
+void CurrentDensityCoefficient2D::zeroCurrentDensity()
+{
+   for (auto &[group, coeff] : group_map)
+   {
+      coeff.constant = 0.0;
+   }
+}
+
+void CurrentDensityCoefficient2D::resetCurrentDensityFromCache()
+{
+   for (auto &[group, value] : cached_inputs)
+   {
+      group_map.at(group).constant = value;
+   }
+}
+
+bool setInputs(CurrentDensityCoefficient2D &current, const MISOInputs &inputs)
+{
+   bool updated = false;
+   for (auto &[group, coeff] : current.group_map)
+   {
+      auto old_const = coeff.constant;
+      std::string cd_group_id = "current_density:" + group;
+      setValueFromInputs(inputs, cd_group_id, coeff.constant);
+      if (coeff.constant != old_const)
+      {
+         updated = true;
+      }
+   }
+
+   for (auto &[input, value] : current.cached_inputs)
+   {
+      auto old_value = value;
+      setValueFromInputs(inputs, input, value);
+      if (value != old_value)
+      {
+         updated = true;
+      }
+   }
+   return updated;
+}
+
+double CurrentDensityCoefficient2D::Eval(mfem::ElementTransformation &trans,
+                                         const mfem::IntegrationPoint &ip)
+{
+   return -1.0 * current_coeff.Eval(trans, ip);
+}
+
+void CurrentDensityCoefficient2D::EvalRevDiff(
+    double Q_bar,
+    mfem::ElementTransformation &trans,
+    const mfem::IntegrationPoint &ip,
+    mfem::DenseMatrix &PointMat_bar)
+{
+   Q_bar *= -1.0;
+   current_coeff.EvalRevDiff(Q_bar, trans, ip, PointMat_bar);
+}
+
+CurrentDensityCoefficient2D::CurrentDensityCoefficient2D(
+    adept::Stack &diff_stack,
+    const nlohmann::json &current_options)
+{
+   for (auto &[group, group_details] : current_options.items())
+   {
+      group_map.emplace(group, mfem::ConstantCoefficient{1.0});
+      auto &group_coeff = group_map.at(group);
+
+      cached_inputs.emplace(group, 0.0);
+
+      for (auto &[source, attrs] : group_details.items())
+      {
+         if (source == "z")
+         {
+            for (const auto &attr : attrs)
+            {
+               source_coeffs.emplace(
+                   attr,
+                   mfem::FunctionCoefficient([](const mfem::Vector &x)
+                                             { return 1.0; },
+                                             [](const mfem::Vector &x,
+                                                const double Q_bar,
+                                                mfem::Vector &x_bar) {}));
+               auto &source_coeff = source_coeffs.at(attr);
+
+               current_coeff.addCoefficient(
+                   attr,
+                   std::make_unique<mfem::ProductCoefficient>(group_coeff,
+                                                              source_coeff));
+            }
+         }
+         else if (source == "-z")
+         {
+            for (const auto &attr : attrs)
+            {
+               // source_coeffs.emplace(attr, mfem::FunctionCoefficient(-1.0));
+               source_coeffs.emplace(
+                   attr,
+                   mfem::FunctionCoefficient([](const mfem::Vector &)
+                                             { return -1.0; },
+                                             [](const mfem::Vector &x,
+                                                const double Q_bar,
+                                                mfem::Vector &x_bar) {}));
+               auto &source_coeff = source_coeffs.at(attr);
+
+               current_coeff.addCoefficient(
+                   attr,
+                   std::make_unique<mfem::ProductCoefficient>(group_coeff,
+                                                              source_coeff));
+            }
+         }
+         else if (source == "box1")
+         {
+            for (const auto &attr : attrs)
+            {
+               source_coeffs.emplace(attr,
+                                     mfem::FunctionCoefficient(
+                                         box1CurrentSource2D,
+                                         [&diff_stack](const mfem::Vector &x,
+                                                       const double &J_bar,
+                                                       mfem::Vector &x_bar) {
+                                            box1CurrentSource2DRevDiff(
+                                                diff_stack, x, J_bar, x_bar);
+                                         }));
+               auto &source_coeff = source_coeffs.at(attr);
+
+               current_coeff.addCoefficient(
+                   attr,
+                   std::make_unique<mfem::ProductCoefficient>(group_coeff,
+                                                              source_coeff));
+            }
+         }
+         else if (source == "box2")
+         {
+            for (const auto &attr : attrs)
+            {
+               source_coeffs.emplace(attr,
+                                     mfem::FunctionCoefficient(
+                                         box2CurrentSource2D,
+                                         [&diff_stack](const mfem::Vector &x,
+                                                       const double &J_bar,
+                                                       mfem::Vector &x_bar) {
+                                            box2CurrentSource2DRevDiff(
+                                                diff_stack, x, J_bar, x_bar);
+                                         }));
+               auto &source_coeff = source_coeffs.at(attr);
+
+               current_coeff.addCoefficient(
+                   attr,
+                   std::make_unique<mfem::ProductCoefficient>(group_coeff,
+                                                              source_coeff));
+            }
+         }
+      }
+   }
+}
+
+}  // namespace miso