From bf4f6bdba933088ffcb96942d8a1e403b2b4b0cd Mon Sep 17 00:00:00 2001 From: Habeeb Idris Date: Sat, 7 Sep 2024 13:16:45 -0400 Subject: [PATCH 1/3] Fixed error in sandbox --- sandbox/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sandbox/CMakeLists.txt b/sandbox/CMakeLists.txt index f3db91a0..aad8272f 100644 --- a/sandbox/CMakeLists.txt +++ b/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 From 41e5a37ce85b8944cf77fc6b3e1838ff60a1e376 Mon Sep 17 00:00:00 2001 From: Habeeb Idris Date: Sat, 7 Sep 2024 13:24:44 -0400 Subject: [PATCH 2/3] Added changes to mesh_warper and thermal_residual.* from idris-etal-aviation-2024 --- miso/mesh_warper.py | 6 +- src/physics/thermal/thermal_residual.cpp | 253 +++++++++++++++++++++-- src/physics/thermal/thermal_residual.hpp | 41 +++- 3 files changed, 275 insertions(+), 25 deletions(-) diff --git a/miso/mesh_warper.py b/miso/mesh_warper.py index 127c47db..a64ac2bd 100644 --- a/miso/mesh_warper.py +++ b/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() diff --git a/src/physics/thermal/thermal_residual.cpp b/src/physics/thermal/thermal_residual.cpp index b76ff669..4a3e175b 100644 --- a/src/physics/thermal/thermal_residual.cpp +++ b/src/physics/thermal/thermal_residual.cpp @@ -1,14 +1,17 @@ +#include +#include #include #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 { @@ -16,11 +19,18 @@ 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); if (residual.load.Size() != 0) { + // std::cout << "residual.load.Size() != 0 and is:\n"; + // residual.load.Print(mfem::out, 25); + // std::cout << "That has been the thermal load vector\n"; + residual.load.SetSubVector(residual.res.getEssentialDofs(), 0.0); } } @@ -70,6 +80,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 +102,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 +138,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 +172,205 @@ ThermalResidual::ThermalResidual( const nlohmann::json &options, const nlohmann::json &materials) : res(fes, fields), + g(std::make_unique( + std::make_unique( + &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(); - 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(); + 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>(); + auto &bc = bcs["convection"]; + std::vector bdr_attr_marker; + double h = 1.0; + double theta_f = 1.0; + if (bc.is_array()) + { + bdr_attr_marker = bc.get>(); + } + else if (bc.is_object()) + { + bdr_attr_marker = bc["attrs"].get>(); + 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 bdr_attr_marker; + double flux = 1.0; + if (bc.is_array()) + { + bdr_attr_marker = bc.get>(); + } + else if (bc.is_object()) + { + bdr_attr_marker = bc["attrs"].get>(); + 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 bdr_attr_marker; + if (bcs["essential"].is_array()) + { + bdr_attr_marker = bcs["essential"].get>(); + } + 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(); + 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 bdr_attr_marker; + if (bcs["weak-essential"].is_array()) + { + bdr_attr_marker = bcs["weak-essential"].get>(); + } + 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(); + 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>(); + + auto mu = options["space-dis"].value("sipg-penalty", -1.0); + if (mu < 0) + { + auto degree = options["space-dis"]["degree"].get(); + mu = pow(degree + 1, 2); + } + + res.addInternalBoundaryFaceIntegrator( + new DGInteriorFaceDiffusionIntegrator(*kappa, mu, -1), + attrs); + + const auto h = intr["h"].get(); + + 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>(); + + auto mu = options["space-dis"].value("sipg-penalty", -1.0); + if (mu < 0) + { + auto degree = options["space-dis"]["degree"].get(); + mu = pow(degree + 1, 2); + } + + res.addInternalBoundaryFaceIntegrator( + new DGInteriorFaceDiffusionIntegrator(*kappa, mu, -1), + attrs); + + const auto h = intr["h"].get(); + const auto theta_f = intr["fluid_temp"].get(); + + 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"; + } + } } } } diff --git a/src/physics/thermal/thermal_residual.hpp b/src/physics/thermal/thermal_residual.hpp index 458ba3d0..c248df39 100644 --- a/src/physics/thermal/thermal_residual.hpp +++ b/src/physics/thermal/thermal_residual.hpp @@ -1,5 +1,5 @@ -#ifndef MISO_THERMAL_RESIDUAL -#define MISO_THERMAL_RESIDUAL +#ifndef MISO_MAGNETOSTATIC_RESIDUAL +#define MISO_MAGNETOSTATIC_RESIDUAL #include #include @@ -47,6 +47,12 @@ class ThermalResidual final mfem::Vector &state_bar, mfem::Vector &adjoint); + friend void finalizeAdjointSystem(ThermalResidual &residual, + mfem::Solver &adj_solver, + const MISOInputs &inputs, + mfem::Vector &state_bar, + mfem::Vector &adjoint); + friend double jacobianVectorProduct(ThermalResidual &residual, const mfem::Vector &wrt_dot, const std::string &wrt); @@ -75,13 +81,15 @@ class ThermalResidual final private: /// Nonlinear form that handles the weak form MISONonlinearForm res; - /// Material dependent coefficient representing thermal conductivity + /// coefficient for weakly imposed boundary conditions + std::unique_ptr g; + // Material dependent coefficient representing thermal conductivity std::unique_ptr kappa; /// Material dependent coefficient representing density std::unique_ptr rho; - /// Material dependent coefficient representing specific heat - /// (at constant volume) - std::unique_ptr cv; + // /// Material dependent coefficient representing specific heat + // /// (at constant volume) + // std::unique_ptr cv; /// Right-hand-side load vector to apply to residual mfem::Vector load; @@ -89,13 +97,26 @@ class ThermalResidual final /// preconditioner for inverting residual's state Jacobian std::unique_ptr prec; - std::unique_ptr constructPreconditioner( + static std::unique_ptr constructPreconditioner( mfem::ParFiniteElementSpace &fes, const nlohmann::json &prec_options) { - auto amg = std::make_unique(); - amg->SetPrintLevel(prec_options["printlevel"].get()); - return amg; + auto prec_type = prec_options["type"].get(); + if (prec_type == "hypreboomeramg") + { + auto amg = std::make_unique(); + amg->SetPrintLevel(prec_options["printlevel"].get()); + return amg; + } + else if (prec_type == "hypreilu") + { + auto ilu = std::make_unique(); + HYPRE_ILUSetType(*ilu, prec_options["ilu-type"]); + HYPRE_ILUSetLevelOfFill(*ilu, prec_options["lev-fill"]); + HYPRE_ILUSetLocalReordering(*ilu, prec_options["ilu-reorder"]); + HYPRE_ILUSetPrintLevel(*ilu, prec_options["printlevel"]); + } + return nullptr; } }; From 7d7550f2a371ef669256060630be5ce83b5b6571 Mon Sep 17 00:00:00 2001 From: Habeeb Idris Date: Mon, 9 Sep 2024 20:00:29 -0400 Subject: [PATCH 3/3] Regrouped headers in thermal_residual.cpp --- src/physics/thermal/thermal_residual.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/physics/thermal/thermal_residual.cpp b/src/physics/thermal/thermal_residual.cpp index 4a3e175b..c2a1ebc6 100644 --- a/src/physics/thermal/thermal_residual.cpp +++ b/src/physics/thermal/thermal_residual.cpp @@ -9,7 +9,6 @@ #include "coefficient.hpp" #include "miso_input.hpp" #include "thermal_integ.hpp" - #include "thermal_residual.hpp" #include "utils.hpp" @@ -27,10 +26,6 @@ void setInputs(ThermalResidual &residual, const miso::MISOInputs &inputs) setVectorFromInputs(inputs, "thermal_load", residual.load); if (residual.load.Size() != 0) { - // std::cout << "residual.load.Size() != 0 and is:\n"; - // residual.load.Print(mfem::out, 25); - // std::cout << "That has been the thermal load vector\n"; - residual.load.SetSubVector(residual.res.getEssentialDofs(), 0.0); } }