Merge branch 'main' into flux_profile_error

Author: ktbolt

Code/Source/solver/Parameters.cpp

@@ -63,6 +63,7 @@
 #include "Parameters.h"
 #include "consts.h"
 #include "LinearAlgebra.h"
+#include "ustruct.h"
 
 #include <iostream>
 #include <regex>
@@ -165,6 +166,9 @@ void Parameters::read_xml(std::string file_name)
   // Set Add_mesh values.
   set_mesh_values(root_element);
 
+  // Set Precomputed_solution values.
+  set_precomputed_solution_values(root_element);
+
   // Set mesh projection parameters.
   set_projection_values(root_element);
 
@@ -217,6 +221,16 @@ void Parameters::set_mesh_values(tinyxml2::XMLElement* root_element)
   }
 }
 
+void Parameters::set_precomputed_solution_values(tinyxml2::XMLElement* root_element)
+{
+  auto add_pre_sol_item = root_element->FirstChildElement(PrecomputedSolutionParameters::xml_element_name_.c_str());
+  if (add_pre_sol_item == nullptr) { 
+    return;
+  }
+
+  precomputed_solution_parameters.set_values(add_pre_sol_item);
+}
+
 void Parameters::set_projection_values(tinyxml2::XMLElement* root_element)
 {
   auto add_proj_item = root_element->FirstChildElement(ProjectionParameters::xml_element_name_.c_str());
@@ -791,6 +805,20 @@ void ConstitutiveModelParameters::set_values(tinyxml2::XMLElement* xml_elem)
   value_set = true;
 }
 
+/// @brief Check if a constitutive model is valid for the given equation.
+//
+void ConstitutiveModelParameters::check_constitutive_model(const Parameter<std::string>& eq_type_str)
+{
+  auto eq_type = consts::equation_name_to_type.at(eq_type_str.value());
+  auto model = consts::constitutive_model_name_to_type.at(type.value());
+
+  if (eq_type == consts::EquationType::phys_ustruct) {
+    if (! ustruct::constitutive_model_is_valid(model)) {
+      throw std::runtime_error("The " + type.value() + " constitutive model is not valid for ustruct equations.");
+    }
+  }
+}
+
 //////////////////////////////////////////////////////////
 //                  CoupleCplBCParameters               //
 //////////////////////////////////////////////////////////
@@ -1864,7 +1892,6 @@ void EquationParameters::set_values(tinyxml2::XMLElement* eq_elem)
   //
   while (item != nullptr) {
     auto name = std::string(item->Value());
-    //std::cout << "[EquationParameters::set_values] name: " << name << std::endl;
 
     if (name == BodyForceParameters::xml_element_name_) {
       auto bf_params = new BodyForceParameters();
@@ -1878,6 +1905,7 @@ void EquationParameters::set_values(tinyxml2::XMLElement* eq_elem)
 
     } else if (name == ConstitutiveModelParameters::xml_element_name_) {
       default_domain->constitutive_model.set_values(item);
+      default_domain->constitutive_model.check_constitutive_model(type);
 
     } else if (name == CoupleCplBCParameters::xml_element_name_) {
       couple_to_cplBC.set_values(item);
@@ -2004,12 +2032,8 @@ GeneralSimulationParameters::GeneralSimulationParameters()
   set_parameter("Starting time step", 0, !required, starting_time_step);
 
   set_parameter("Time_step_size", 0.0, required, time_step_size);
-  set_parameter("Precomputed_time_step_size", 0.0, !required, precomputed_time_step_size);
   set_parameter("Verbose", false, !required, verbose);
   set_parameter("Warning", false, !required, warning);
-  set_parameter("Use_precomputed_solution", false, !required, use_precomputed_solution);
-  set_parameter("Precomputed_solution_file_path", "", !required, precomputed_solution_file_path);
-  set_parameter("Precomputed_solution_field_name", "", !required, precomputed_solution_field_name);
 }
 
 void GeneralSimulationParameters::print_parameters()
@@ -2301,6 +2325,40 @@ void MeshParameters::set_values(tinyxml2::XMLElement* mesh_elem)
   }
 }
 
+/////////////////////////////////////////////////////////////////////////////
+//        P r e c o m p u t e d S o l u t i o n P a r a m e t e r s       //
+/////////////////////////////////////////////////////////////////////////////
+
+// The PrecomputedSolutionParameters class stores parameters for the
+// 'Precomputed_solution' XML element used to read in the data from a 
+// precomputed solution for the simulation state.
+
+const std::string PrecomputedSolutionParameters::xml_element_name_ = "Precomputed_solution";
+
+PrecomputedSolutionParameters::PrecomputedSolutionParameters()
+{
+  // A parameter that must be defined.
+  bool required = true;
+
+  set_parameter("Field_name", "", required, field_name);
+  set_parameter("File_path", "", required, file_path);
+  set_parameter("Time_step", 0.0, !required, time_step);
+  set_parameter("Use_precomputed_solution", false, !required, use_precomputed_solution);
+}
+
+void PrecomputedSolutionParameters::set_values(tinyxml2::XMLElement* xml_elem)
+{
+  using namespace tinyxml2;
+  std::string error_msg = "Unknown " + xml_element_name_ + " XML element '";
+  using std::placeholders::_1;
+  using std::placeholders::_2;
+
+  std::function<void(const std::string&, const std::string&)> ftpr =
+      std::bind( &PrecomputedSolutionParameters::set_parameter_value, *this, _1, _2);
+
+  xml_util_set_parameters(ftpr, xml_elem, error_msg);
+}
+
 //////////////////////////////////////////////////////////
 //        P r o j e c t i o n P a r a m e t e r s       //
 //////////////////////////////////////////////////////////

Code/Source/solver/Parameters.h

@@ -510,6 +510,7 @@ class ConstitutiveModelParameters : public ParameterLists
   public:
     ConstitutiveModelParameters();
     void print_parameters();
+    void check_constitutive_model(const Parameter<std::string>& eq_type);
     bool defined() const { return value_set; };
     void set_values(tinyxml2::XMLElement* modl_params);
     static const std::string xml_element_name_;
@@ -752,6 +753,30 @@ class OutputParameters : public ParameterLists
     std::vector<Parameter<std::string>> alias_list;
 };
 
+/// @brief The PrecomputedSolutionParameters class stores parameters for the
+/// 'Precomputed_solution' XML element used to read in the data from a precomputed solution 
+/// for the simulation state
+/// \code {.xml}
+/// <Precomputed_solution>
+///   <Project_from_face> lumen_wall </Project_from_face>
+/// </Precomputed_solution>
+/// \endcode
+
+class PrecomputedSolutionParameters: public ParameterLists
+{
+  public:
+    PrecomputedSolutionParameters();
+
+    void set_values(tinyxml2::XMLElement* xml_elem);
+
+    static const std::string xml_element_name_;
+
+    Parameter<std::string> field_name;
+    Parameter<std::string> file_path;
+    Parameter<double> time_step;
+    Parameter<bool> use_precomputed_solution;
+};
+
 /// @brief The ProjectionParameters class stores parameters for the
 /// 'Add_projection' XML element used for fluid-structure interaction 
 /// simulations.
@@ -1303,11 +1328,9 @@ class GeneralSimulationParameters : public ParameterLists
     Parameter<bool> start_averaging_from_zero;
     Parameter<bool> verbose;
     Parameter<bool> warning;
-    Parameter<bool> use_precomputed_solution;
 
     Parameter<double> spectral_radius_of_infinite_time_step;
     Parameter<double> time_step_size;
-    Parameter<double> precomputed_time_step_size;
 
     Parameter<int> increment_in_saving_restart_files;
     Parameter<int> increment_in_saving_vtk_files;
@@ -1322,8 +1345,6 @@ class GeneralSimulationParameters : public ParameterLists
     Parameter<std::string> searched_file_name_to_trigger_stop; 
     Parameter<std::string> save_results_in_folder; 
     Parameter<std::string> simulation_initialization_file_path;
-    Parameter<std::string> precomputed_solution_file_path;
-    Parameter<std::string> precomputed_solution_field_name;
 };
 
 /// @brief The FaceParameters class is used to store parameters for the
@@ -1427,6 +1448,7 @@ class Parameters {
     void set_contact_values(tinyxml2::XMLElement* root_element);
     void set_equation_values(tinyxml2::XMLElement* root_element);
     void set_mesh_values(tinyxml2::XMLElement* root_element);
+    void set_precomputed_solution_values(tinyxml2::XMLElement* root_element);
     void set_projection_values(tinyxml2::XMLElement* root_element);
 
     // Objects representing each parameter section of XML file.
@@ -1435,6 +1457,7 @@ class Parameters {
     std::vector<MeshParameters*> mesh_parameters;
     std::vector<EquationParameters*> equation_parameters;
     std::vector<ProjectionParameters*> projection_parameters;
+    PrecomputedSolutionParameters precomputed_solution_parameters;
 };
 
 #endif

Code/Source/solver/Simulation.cpp

@@ -100,10 +100,12 @@ void Simulation::set_module_parameters()
   com_mod.stFileIncr = general.increment_in_saving_restart_files.value();
   com_mod.rmsh.isReqd = general.simulation_requires_remeshing.value();
 
-  com_mod.usePrecomp = general.use_precomputed_solution.value();
-  com_mod.precompFileName = general.precomputed_solution_file_path.value();
-  com_mod.precompFieldName = general.precomputed_solution_field_name.value();
-  com_mod.precompDt = general.precomputed_time_step_size.value();
+  auto& precomp_sol = parameters.precomputed_solution_parameters;
+  com_mod.usePrecomp = precomp_sol.use_precomputed_solution.value();
+  com_mod.precompFileName = precomp_sol.file_path.value();
+  com_mod.precompFieldName = precomp_sol.field_name.value();
+  com_mod.precompDt = precomp_sol.time_step.value();
+
   if ((com_mod.precompDt == 0.0) && (com_mod.usePrecomp)) {
     std::cout << "Precomputed time step size is zero. Setting to simulation time step size." << std::endl;
     com_mod.precompDt = com_mod.dt;

Code/Source/solver/initialize.cpp

@@ -328,19 +328,21 @@ void initialize(Simulation* simulation, Vector<double>& timeP)
   dmsg.banner();
   #endif
 
-  std::string hist_file_name;
+  // Setup logging to history file.
+  if (!simulation->com_mod.cm.slv(simulation->cm_mod)) {
+    std::string hist_file_name;
 
-  if (chnl_mod.appPath != "") { 
-    auto mkdir_arg = std::string("mkdir -p ") + chnl_mod.appPath;
-    std::system(mkdir_arg.c_str());
-    hist_file_name = chnl_mod.appPath + "/" + simulation->history_file_name;
-  } else {
-    hist_file_name = simulation->history_file_name;
-  }
+    if (chnl_mod.appPath != "") { 
+      auto mkdir_arg = std::string("mkdir -p ") + chnl_mod.appPath;
+      std::system(mkdir_arg.c_str());
+      hist_file_name = chnl_mod.appPath + "/" + simulation->history_file_name;
+    } else {
+      hist_file_name = simulation->history_file_name;
+    }
 
-  // Setup logging to history file.
-  bool output_to_cout = true;
-  simulation->logger.initialize(hist_file_name, output_to_cout);
+    bool output_to_cout = true;
+    simulation->logger.initialize(hist_file_name, output_to_cout);
+  }
 
   #ifdef debug_initialize
   dmsg << "Set time " << std::endl;

Code/Source/solver/mat_models.cpp

@@ -408,21 +408,13 @@ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType&
   // Now, add isochoric and total stress, elasticity tensors
   switch (stM.isoType) {
     case ConstitutiveModelType::stIso_lin: {
-      if (ustruct) {
-        throw std::runtime_error("[compute_pk2cc] Linear isotropic material model not valid for ustruct physics.");
-      }
-
       double g1 = stM.C10;    // mu
       S += g1*Idm;
       return; 
     } break;
 
     // St.Venant-Kirchhoff
     case ConstitutiveModelType::stIso_StVK: {
-      if (ustruct) {
-        throw std::runtime_error("[compute_pk2cc] St.Venant-Kirchhoff material model not valid for ustruct physics.");
-      }
-      
       double g1 = stM.C10;         // lambda
       double g2 = stM.C01 * 2.0;   // 2*mu
 
@@ -432,10 +424,6 @@ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType&
 
     // modified St.Venant-Kirchhoff
     case ConstitutiveModelType::stIso_mStVK: {
-      if (ustruct) {
-        throw std::runtime_error("[compute_pk2cc] Modified St.Venant-Kirchhoff material model not valid for ustruct physics.");
-      }
-
       double g1 = stM.C10; // kappa
       double g2 = stM.C01;  // mu
 

Code/Source/solver/read_files.cpp

@@ -122,7 +122,7 @@ void face_match(ComMod& com_mod, faceType& lFa, faceType& gFa, Vector<int>& ptr)
     blk[iBlk].n = blk[iBlk].n + 1;
   }
 
-  for (int a = 0; a < gFa.nNo; a++) {
+  for (int a = 0; a < lFa.nNo; a++) {
     auto coord = lFa.x.col(a);
     int iBlk = find_blk(nsd, nBlkd, nFlt, xMin, dx, coord);
     auto minS = std::numeric_limits<double>::max();
@@ -1744,9 +1744,8 @@ void read_files(Simulation* simulation, const std::string& file_name)
 
     if (eq.phys == EquationType::phys_heatF) {   
       auto& eq1_params = simulation->parameters.equation_parameters[0];
-      auto& general_params = simulation->parameters.general_simulation_parameters;
       auto eq1_type = eq1_params->type.value();
-      if ((eq1_type != "fluid") && (eq1_type != "FSI") && (!general_params.use_precomputed_solution.value())) {
+      if ((eq1_type != "fluid") && (eq1_type != "FSI") && (!com_mod.usePrecomp)) {
         throw std::runtime_error("heatF equation has to be specified after fluid/FSI equation");
       }
     }     

Code/Source/solver/set_bc.cpp

@@ -84,8 +84,7 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod)
   // if (ALL(cplBC.fa.bGrp .EQ. cplBC_Dir)) RETURN
 
   // Determine current physics
-  auto& cDmn = com_mod.cDmn;
-  auto cPhys = eq.dmn[cDmn].phys;
+  auto cPhys = eq.phys;
 
   // Mechanical configuration in which to compute flowrate
   auto cfg_o = MechanicalConfigurationType::reference;
@@ -686,8 +685,7 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod)
   auto& eq = com_mod.eq[iEq];
 
   // Determine current physics
-  auto& cDmn = com_mod.cDmn;
-  auto cPhys = eq.dmn[cDmn].phys;
+  auto cPhys = eq.phys;
 
   // Configuration in which to compute flowrate
   auto cfg_o = MechanicalConfigurationType::reference;

Code/Source/solver/ustruct.cpp

@@ -210,6 +210,21 @@ void b_ustruct_3d(const ComMod& com_mod, const int eNoN, const double w, const V
   }
 }
 
+/// @brief Check is a constitutive model is valid for the ustruct equation.
+//
+bool constitutive_model_is_valid(consts::ConstitutiveModelType model)
+{ 
+  using namespace consts;
+    
+  static std::set<ConstitutiveModelType> unsupported_models { 
+    ConstitutiveModelType::stIso_lin,
+    ConstitutiveModelType::stIso_StVK,
+    ConstitutiveModelType::stIso_mStVK
+  };
+
+  return unsupported_models.count(model) == 0;
+}   
+
 /// @brief Reproduces Fortran CONSTRUCT_uSOLID.
 //
 void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array<double>& Ag, 

Code/Source/solver/ustruct.h

@@ -43,6 +43,8 @@ void b_ustruct_3d(const ComMod& com_mod, const int eNoN, const double w, const V
     const Array<double>& Nx, const Array<double>& dl, const Vector<double>& hl, const Vector<double>& nV, 
     Array<double>& lR, Array3<double>& lK, Array3<double>& lKd);
 
+bool constitutive_model_is_valid(consts::ConstitutiveModelType model);
+
 void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array<double>& Ag, const Array<double>& Yg, 
     const Array<double>& Dg);
 

tests/cases/fluid/precomputed_dye_AD/solver.xml

@@ -6,9 +6,7 @@
   <Number_of_spatial_dimensions> 3 </Number_of_spatial_dimensions> 
   <Number_of_time_steps> 1 </Number_of_time_steps> 
   <Time_step_size> 0.01 </Time_step_size> 
-  <Use_precomputed_solution> true </Use_precomputed_solution>
-  <Precomputed_solution_file_path> precomputed_velocity.vtu </Precomputed_solution_file_path>
-  <Precomputed_solution_field_name> Velocity </Precomputed_solution_field_name>
+
   <Spectral_radius_of_infinite_time_step> 0.50 </Spectral_radius_of_infinite_time_step> 
   <Searched_file_name_to_trigger_stop> STOP_SIM </Searched_file_name_to_trigger_stop> 
 
@@ -25,6 +23,12 @@
   <Debug> 1 </Debug> 
 </GeneralSimulationParameters>
 
+<Precomputed_solution> 
+  <Use_precomputed_solution> true </Use_precomputed_solution>
+  <File_path> precomputed_velocity.vtu </File_path>
+  <Field_name> Velocity </Field_name>
+</Precomputed_solution> 
+
 <Add_mesh name="msh" > 
 
   <Mesh_file_path> mesh/mesh-complete.mesh.vtu </Mesh_file_path>