From 76950bc12a94eff559063e0504731e38bd392c93 Mon Sep 17 00:00:00 2001 From: Taeouk Kim Date: Tue, 4 Aug 2026 16:07:03 -0700 Subject: [PATCH] Initialize the 3D pressure field from a scalar value --- Code/Source/solver/ComMod.h | 6 +++++ Code/Source/solver/Parameters.cpp | 2 ++ Code/Source/solver/Parameters.h | 1 + Code/Source/solver/Simulation.cpp | 5 +++++ Code/Source/solver/baf_ini.cpp | 19 +++++++++------- Code/Source/solver/distribute.cpp | 2 ++ Code/Source/solver/initialize.cpp | 25 ++++++++++++++------- Code/Source/solver/svOneD_interface.cpp | 20 ++++++++++++++--- Code/Source/solver/svOneD_interface.h | 3 ++- Code/Source/solver/svZeroD_interface.cpp | 28 +++++++++++++++++++++++- Code/Source/solver/svZeroD_interface.h | 3 ++- 11 files changed, 92 insertions(+), 22 deletions(-) diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 2e9f60068..73c9f2b0c 100644 --- a/Code/Source/solver/ComMod.h +++ b/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 -----// diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 8f9eed3d6..23e1b80d9 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/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); diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index b18bd4661..920f92490 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1827,6 +1827,7 @@ class GeneralSimulationParameters : public ParameterLists Parameter spectral_radius_of_infinite_time_step; Parameter time_step_size; + Parameter initial_pressure_scalar; Parameter include_xml; Parameter increment_in_saving_restart_files; diff --git a/Code/Source/solver/Simulation.cpp b/Code/Source/solver/Simulation.cpp index 7e6d56523..714e6b1fa 100644 --- a/Code/Source/solver/Simulation.cpp +++ b/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(); diff --git a/Code/Source/solver/baf_ini.cpp b/Code/Source/solver/baf_ini.cpp index e83dd8a64..00d470da5 100644 --- a/Code/Source/solver/baf_ini.cpp +++ b/Code/Source/solver/baf_ini.cpp @@ -149,13 +149,16 @@ 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) { @@ -163,11 +166,11 @@ void baf_ini(Simulation* simulation, SolutionStates& solutions) } 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 diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 00986613e..51fc75626 100644 --- a/Code/Source/solver/distribute.cpp +++ b/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); diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index 7baf4f330..ae57c98c2 100644 --- a/Code/Source/solver/initialize.cpp +++ b/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) { diff --git a/Code/Source/solver/svOneD_interface.cpp b/Code/Source/solver/svOneD_interface.cpp index db1724e6e..a04e66cfe 100644 --- a/Code/Source/solver/svOneD_interface.cpp +++ b/Code/Source/solver/svOneD_interface.cpp @@ -12,6 +12,7 @@ #include #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) ----- diff --git a/Code/Source/solver/svOneD_interface.h b/Code/Source/solver/svOneD_interface.h index ce775fd57..4ef528cec 100644 --- a/Code/Source/solver/svOneD_interface.h +++ b/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. /// diff --git a/Code/Source/solver/svZeroD_interface.cpp b/Code/Source/solver/svZeroD_interface.cpp index 2a3c8da8f..dbcee84ab 100644 --- a/Code/Source/solver/svZeroD_interface.cpp +++ b/Code/Source/solver/svZeroD_interface.cpp @@ -226,7 +226,8 @@ void print_svZeroD(int* nSrfs, const std::vector& 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). diff --git a/Code/Source/solver/svZeroD_interface.h b/Code/Source/solver/svZeroD_interface.h index 1de624f85..a6f1f3629 100644 --- a/Code/Source/solver/svZeroD_interface.h +++ b/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& 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);