Initialize the 3D pressure field from a scalar value

Code/Source/solver/ComMod.h

@@ -1663,6 +1663,10 @@ class ComMod {
 
     /// @brief Whether to use precomputed state-variable solutions
     bool usePrecomp = false;
+
+    /// @brief Whether a scalar initial pressure is provided in solver.xml
+    bool have_initial_pressure_scalar = false;
+
     //----- int members -----//
 
     /// @brief Current domain
@@ -1747,6 +1751,8 @@ class ComMod {
     /// @brief Time
     double time = 0.0;
 
+    /// @brief Scalar pressure used to initialize the 3D pressure field
+    double initial_pressure_scalar = 0.0;
 
     //----- string members -----//
 

Code/Source/solver/Parameters.cpp

@@ -2846,6 +2846,8 @@ GeneralSimulationParameters::GeneralSimulationParameters() {
                 start_saving_after_time_step);
   set_parameter("Starting time step", 0, !required, starting_time_step);
 
+  set_parameter("Initial_pressure_scalar", 0.0, !required,
+                initial_pressure_scalar);
   set_parameter("Time_step_size", 0.0, required, time_step_size);
   set_parameter("Verbose", false, !required, verbose);
   set_parameter("Warning", false, !required, warning);

Code/Source/solver/Parameters.h

@@ -1827,6 +1827,7 @@ class GeneralSimulationParameters : public ParameterLists
 
     Parameter<double> spectral_radius_of_infinite_time_step;
     Parameter<double> time_step_size;
+    Parameter<double> initial_pressure_scalar;
 
     Parameter<std::string> include_xml;
     Parameter<int> increment_in_saving_restart_files;

Code/Source/solver/Simulation.cpp

@@ -61,6 +61,11 @@ void Simulation::set_module_parameters()
   com_mod.startTS = general.starting_time_step.value();
   com_mod.dt = general.time_step_size.value();
 
+  com_mod.have_initial_pressure_scalar = general.initial_pressure_scalar.defined();
+  if (com_mod.have_initial_pressure_scalar) {
+    com_mod.initial_pressure_scalar = general.initial_pressure_scalar.value();
+  }
+
   com_mod.stopTrigName = general.searched_file_name_to_trigger_stop.value();
   com_mod.ichckIEN = general.check_ien_order.value();
   com_mod.saveVTK = general.save_results_to_vtk_format.value();

Code/Source/solver/baf_ini.cpp

@@ -149,25 +149,28 @@ void baf_ini(Simulation* simulation, SolutionStates& solutions)
       }
     }
 
+    // Create temporary SolutionStates from the initialized 3D state.
+    SolutionStates init_solutions;
+    init_solutions.old.get_acceleration() = Ao;
+    init_solutions.old.get_displacement() = Do;
+    init_solutions.old.get_velocity() = Yo;
+    init_solutions.current.get_acceleration() = Ao;
+    init_solutions.current.get_displacement() = Do;
+    init_solutions.current.get_velocity() = Yo;
     if (!com_mod.stFileFlag) {
-      // Create temporary SolutionStates for set_bc calls
-      SolutionStates temp_solutions;
-      temp_solutions.old.get_acceleration() = Ao;
-      temp_solutions.old.get_displacement() = Do;
-      temp_solutions.old.get_velocity() = Yo;
-      set_bc::rcr_init(com_mod, cm_mod, temp_solutions);
+      set_bc::rcr_init(com_mod, cm_mod, init_solutions);
     }
 
     if (com_mod.cplBC.useGenBC) {
       set_bc::genBC_Integ_X(com_mod, cm_mod, "I");
     }
 
     if (com_mod.cplBC.useSvZeroD) {
-      svZeroD::init_svZeroD(com_mod, cm_mod);
+      svZeroD::init_svZeroD(com_mod, cm_mod, init_solutions);
     }
 
     if (com_mod.cplBC.useSvOneD) {
-      svOneD::init_svOneD(com_mod, cm_mod);
+      svOneD::init_svOneD(com_mod, cm_mod, init_solutions);
     }
 
     // Initialize cap integration for Coupled boundary conditions

Code/Source/solver/distribute.cpp

@@ -331,6 +331,8 @@ void distribute(Simulation* simulation)
     cm.bcast(cm_mod, &com_mod.nEq);
     cm.bcast(cm_mod, &com_mod.dt);
     cm.bcast(cm_mod, &com_mod.precompDt);
+    cm.bcast(cm_mod, &com_mod.have_initial_pressure_scalar);
+    cm.bcast(cm_mod, &com_mod.initial_pressure_scalar);
 
     cm.bcast(cm_mod, &com_mod.zeroAve);
     cm.bcast(cm_mod, &com_mod.cmmInit);

Code/Source/solver/initialize.cpp

@@ -976,14 +976,23 @@ void zero_init(Simulation* simulation, SolutionStates& solutions)
   }
 
   if (com_mod.Pinit.size() != 0) {
-     #ifdef debug_zero_init
-     dmsg << "Initialize Yo to provided P solution";
-     #endif
-     for (int a = 0; a < com_mod.tnNo; a++) {
-       for (int i = 0; i < nsd; i++) {
-         Yo(nsd,a) = com_mod.Pinit(a);
-       }
-     }
+    #ifdef debug_zero_init
+    dmsg << "Initialize Yo to provided P solution";
+    #endif
+    for (int a = 0; a < com_mod.tnNo; a++) {
+      Yo(nsd,a) = com_mod.Pinit(a);
+    }
+
+  } else if (com_mod.have_initial_pressure_scalar) {
+    #ifdef debug_zero_init
+    dmsg << "Initialize Yo to scalar initial pressure from solver.xml";
+    #endif
+    for (int a = 0; a < com_mod.tnNo; a++) {
+      Yo(nsd,a) = com_mod.initial_pressure_scalar;
+    }
+
+    std::cout << "Initialize 3D pressure field to scalar value: "
+              << com_mod.initial_pressure_scalar << std::endl;
   }
 
   if (com_mod.Dinit.size() != 0) {

Code/Source/solver/svOneD_interface.cpp

@@ -12,6 +12,7 @@
 #include <filesystem>
 
 #include "ComMod.h"
+#include "all_fun.h"
 #include "consts.h"
 #include "utils.h"
 #include "svOneD_interface/OneDSolverInterface.h"
@@ -112,7 +113,8 @@ static std::string resolve_lib_path(const std::string& lib_base)
 // ---------------------------------------------------------------------------
 // init_svOneD
 // ---------------------------------------------------------------------------
-void init_svOneD(ComMod& com_mod, const CmMod& cm_mod)
+void init_svOneD(ComMod& com_mod, const CmMod& cm_mod,
+                 const SolutionStates& solutions)
 {
   using namespace consts;
 
@@ -152,6 +154,18 @@ void init_svOneD(ComMod& com_mod, const CmMod& cm_mod)
       st.ramp_steps         = bc.coupled_bc.get_oned_ramp_steps();
       st.ramp_ref_pressure  = bc.coupled_bc.get_oned_ramp_ref_pressure();
       st.relax_factor       = bc.coupled_bc.get_oned_relax_factor();
+      const auto& face = com_mod.msh[bc.iM].fa[bc.iFa];
+      const double area = face.area;
+      if (area > 0.0) {
+        const auto& Yo = solutions.old.get_velocity();
+        const double P_init =
+            all_fun::integ(com_mod, cm_mod, face, Yo, com_mod.nsd,
+                           solutions, std::nullopt, false,
+                           MechanicalConfigurationType::reference) / area;
+        st.P_prev_sent_old = P_init;
+        st.P_prev_sent_new = P_init;
+        st.P_neu_prev = P_init;
+      }
       oned_models.push_back(std::move(st));
     }
   }
@@ -219,8 +233,8 @@ void init_svOneD(ComMod& com_mod, const CmMod& cm_mod)
     st.solution.resize(system_size, 0.0);
     shared_lib_instance->return_solution(problem_id, st.solution.data(), system_size);
 
-    // Initial coupled value = 0; first calc_svOneD call sets the real value.
-    eq.bc[st.iBc].coupled_bc.set_pressure(0.0);
+    // Keep the BC pressure state consistent with the seeded pressure history.
+    eq.bc[st.iBc].coupled_bc.set_pressure(st.P_neu_prev);
   }
 
   // ----- Broadcast metadata for all models (Phase 2: batch exchange) -----

Code/Source/solver/svOneD_interface.h

@@ -63,7 +63,8 @@ namespace svOneD {
 
 /// @brief Initialize the 1D solver and populate the initial cplBC state.
 /// Called once from baf_ini() after the BC data structures are set up.
-void init_svOneD(ComMod& com_mod, const CmMod& cm_mod);
+void init_svOneD(ComMod& com_mod, const CmMod& cm_mod,
+                 const SolutionStates& solutions);
 
 /// @brief Advance the 1D solver by one time step and update the coupled BC value.
 ///

Code/Source/solver/svZeroD_interface.cpp

@@ -226,7 +226,8 @@ void print_svZeroD(int* nSrfs, const std::vector<int>& surfID, double Q[], doubl
 // init_svZeroD
 //--------------
 //
-void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod) 
+void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod,
+                  const SolutionStates& solutions)
 {
   using namespace consts;
 
@@ -359,6 +360,31 @@ void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod)
     }
   }
 
+  for (int s = 0; s < numCoupledSrfs; ++s) {
+    bcType* bc = nullptr;
+    if (!nth_coupled_bc(com_mod, s, &bc)) {
+      throw std::runtime_error(
+          "ERROR: [init_svZeroD] Internal error resolving Coupled BC "
+          "for initialization history.");
+    }
+
+    // For cap coupling in parallel, flow integration requires all ranks.
+    if (cm.seq() || !bc->coupled_bc.has_cap()) {
+      bc->coupled_bc.compute_flowrates(com_mod, cm_mod, solutions);
+    }
+
+    bc->coupled_bc.compute_pressures(com_mod, cm_mod, solutions);
+
+    const double Q_init = bc->coupled_bc.get_Qo();
+    const double P_init = bc->coupled_bc.get_Po();
+
+    bc->coupled_bc.set_Q_prev_sent(Q_init);
+    bc->coupled_bc.set_Q_input_prev(Q_init, Q_init);
+    bc->coupled_bc.set_P_prev_sent(P_init, P_init);
+    bc->coupled_bc.set_P_neu_prev(P_init);
+    bc->coupled_bc.set_pressure(P_init);
+  }
+
   // Broadcast initial values to follower processes
   if (!cm.seq()) {
     // Coupled BCs - broadcast Neumann pressures (one scalar bcast per BC).

Code/Source/solver/svZeroD_interface.h

@@ -17,7 +17,8 @@ void get_coupled_QP(ComMod& com_mod, double QCoupled[], double QnCoupled[], doub
 
 void print_svZeroD(int* nSrfs, const std::vector<int>& surfID, double Q[], double P[]);
 
-void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod);
+void init_svZeroD(ComMod& com_mod, const CmMod& cm_mod,
+                  const SolutionStates& solutions);
 
 void calc_svZeroD(ComMod& com_mod, const CmMod& cm_mod, char BCFlag);