Issue 546 pressure initialization

Code/Source/solver/CMakeLists.txt

@@ -204,6 +204,7 @@ set(CSRCS
   stokes.h stokes.cpp
   sv_struct.h sv_struct.cpp
   svZeroD_interface.h svZeroD_interface.cpp
+  svOneD_subroutines.h svOneD_subroutines.cpp
   txt.h txt.cpp
   utils.h utils.cpp
   ustruct.h ustruct.cpp
@@ -228,6 +229,7 @@ set(CSRCS
   SPLIT.c
 
   svZeroD_interface/LPNSolverInterface.h svZeroD_interface/LPNSolverInterface.cpp
+  svOneD_interface/OneDSolverInterface.h svOneD_interface/OneDSolverInterface.cpp
 
   BoundaryCondition.h BoundaryCondition.cpp
   RobinBoundaryCondition.h RobinBoundaryCondition.cpp

Code/Source/solver/ComMod.cpp

@@ -217,3 +217,13 @@ void svZeroDSolverInterfaceType::set_data(const svZeroDSolverInterfaceParameters
 
   has_data = true;
 }
+
+void svOneDSolverInterfaceType::set_data(const svOneDSolverInterfaceParameters& params)
+{
+  if (!params.defined()) {
+    return;
+  }
+
+  solver_library = params.shared_library();
+  has_data       = true;
+}

Code/Source/solver/ComMod.h

@@ -204,6 +204,10 @@ class bcType
 
     // Coupled BC class
     CoupledBoundaryCondition coupled_bc;
+
+    // svOneD: per-face 1D solver input file path (set when Time_dependence=Coupled
+    // and svOneDSolver_interface is active).
+    std::string oned_input_file;
 };
 
 /// @brief Class storing data for B-Splines.
@@ -771,6 +775,12 @@ class cplFaceType
 
     // RCR type BC
     rcrType RCR;
+
+    // svOneD: path to the per-face 1D solver input file.
+    std::string oned_input_file;
+
+    // Whether this face uses RCR (Windkessel) boundary condition.
+    bool isRCR = false;
 };
 
 //----------------------------
@@ -807,6 +817,26 @@ class svZeroDSolverInterfaceType
     void set_data(const svZeroDSolverInterfaceParameters& params);
 };
 
+//----------------------------
+// svOneDSolverInterfaceType
+//----------------------------
+// This class stores information used to interface to the svOneDSolver.
+//
+class svOneDSolverInterfaceType
+{
+  public:
+    // Path to the 1D solver shared library (without .so/.dylib extension,
+    // or with extension if the full path is provided).
+    std::string solver_library;
+
+    // If the data has been set for the interface.  Set to true after
+    // the svOneDSolver_interface XML element has been parsed.
+    // Note: the per-face input files are stored in cplFaceType::oned_input_file.
+    bool has_data = false;
+
+    void set_data(const svOneDSolverInterfaceParameters& params);
+};
+
 /// @brief For coupled 0D-3D problems
 //
 class cplBCType
@@ -822,6 +852,9 @@ class cplBCType
     //  Whether to use svZeroD
     bool useSvZeroD = false;
 
+    //  Whether to use svOneD (svOneDSolver)
+    bool useSv1D = false;
+
     //  Whether to initialize RCR from flow data
     bool initRCR = false;
 
@@ -853,6 +886,8 @@ class cplBCType
 
     svZeroDSolverInterfaceType svzerod_solver_interface;
 
+    svOneDSolverInterfaceType sv1d_solver_interface;
+
     /// @brief The name of history file containing "X"
     std::string saveName;
     //std::string(LEN=stdL) :: saveName = "LPN.dat";
@@ -1633,6 +1668,9 @@ class ComMod {
     bool usePrecomp = false;
     //----- int members -----//
 
+    /// @brief scalar initial pressure is explicitly provided in the input file
+    bool have_initial_pressure = false;
+
     /// @brief Current domain
     int cDmn = 0;
 
@@ -1715,6 +1753,9 @@ class ComMod {
     /// @brief Time
     double time = 0.0;
 
+    /// @brief initial scalar pressure value
+    double initial_pressure = 0.0;
+
 
     //----- string members -----//
 

Code/Source/solver/CoupledBoundaryCondition.cpp

@@ -22,6 +22,17 @@ CoupledBoundaryCondition::CoupledBoundaryCondition(const CoupledBoundaryConditio
     , bc_type_(other.bc_type_)
     , block_name_(other.block_name_)
     , face_name_(other.face_name_)
+    , oned_input_file_(other.oned_input_file_)
+    , oned_ramp_steps_(other.oned_ramp_steps_)
+    , oned_ramp_ref_pressure_(other.oned_ramp_ref_pressure_)
+    , oned_relax_factor_(other.oned_relax_factor_)
+    , ramp_step_count_(other.ramp_step_count_)
+    , P_prev_sent_old_(other.P_prev_sent_old_)
+    , P_prev_sent_new_(other.P_prev_sent_new_)
+    , Q_prev_sent_(other.Q_prev_sent_)
+    , P_neu_prev_(other.P_neu_prev_)
+    , Q_input_prev_old_(other.Q_input_prev_old_)
+    , Q_input_prev_new_(other.Q_input_prev_new_)
     , Qo_(other.Qo_)
     , Qn_(other.Qn_)
     , Po_(other.Po_)
@@ -53,6 +64,17 @@ CoupledBoundaryCondition& CoupledBoundaryCondition::operator=(const CoupledBound
         bc_type_ = other.bc_type_;
         block_name_ = other.block_name_;
         face_name_ = other.face_name_;
+        oned_input_file_ = other.oned_input_file_;
+        oned_ramp_steps_ = other.oned_ramp_steps_;
+        oned_ramp_ref_pressure_ = other.oned_ramp_ref_pressure_;
+        oned_relax_factor_ = other.oned_relax_factor_;
+        ramp_step_count_ = other.ramp_step_count_;
+        P_prev_sent_old_ = other.P_prev_sent_old_;
+        P_prev_sent_new_ = other.P_prev_sent_new_;
+        Q_prev_sent_ = other.Q_prev_sent_;
+        P_neu_prev_ = other.P_neu_prev_;
+        Q_input_prev_old_ = other.Q_input_prev_old_;
+        Q_input_prev_new_ = other.Q_input_prev_new_;
         Qo_ = other.Qo_;
         Qn_ = other.Qn_;
         Po_ = other.Po_;
@@ -83,6 +105,17 @@ CoupledBoundaryCondition::CoupledBoundaryCondition(CoupledBoundaryCondition&& ot
     , bc_type_(other.bc_type_)
     , block_name_(std::move(other.block_name_))
     , face_name_(std::move(other.face_name_))
+    , oned_input_file_(std::move(other.oned_input_file_))
+    , oned_ramp_steps_(other.oned_ramp_steps_)
+    , oned_ramp_ref_pressure_(other.oned_ramp_ref_pressure_)
+    , oned_relax_factor_(other.oned_relax_factor_)
+    , ramp_step_count_(other.ramp_step_count_)
+    , P_prev_sent_old_(other.P_prev_sent_old_)
+    , P_prev_sent_new_(other.P_prev_sent_new_)
+    , Q_prev_sent_(other.Q_prev_sent_)
+    , P_neu_prev_(other.P_neu_prev_)
+    , Q_input_prev_old_(other.Q_input_prev_old_)
+    , Q_input_prev_new_(other.Q_input_prev_new_)
     , Qo_(other.Qo_)
     , Qn_(other.Qn_)
     , Po_(other.Po_)
@@ -129,6 +162,17 @@ CoupledBoundaryCondition& CoupledBoundaryCondition::operator=(CoupledBoundaryCon
         bc_type_ = other.bc_type_;
         block_name_ = std::move(other.block_name_);
         face_name_ = std::move(other.face_name_);
+        oned_input_file_ = std::move(other.oned_input_file_);
+        oned_ramp_steps_ = other.oned_ramp_steps_;
+        oned_ramp_ref_pressure_ = other.oned_ramp_ref_pressure_;
+        oned_relax_factor_ = other.oned_relax_factor_;
+        ramp_step_count_ = other.ramp_step_count_;
+        P_prev_sent_old_ = other.P_prev_sent_old_;
+        P_prev_sent_new_ = other.P_prev_sent_new_;
+        Q_prev_sent_ = other.Q_prev_sent_;
+        P_neu_prev_ = other.P_neu_prev_;
+        Q_input_prev_old_ = other.Q_input_prev_old_;
+        Q_input_prev_new_ = other.Q_input_prev_new_;
         Qo_ = other.Qo_;
         Qn_ = other.Qn_;
         Po_ = other.Po_;
@@ -208,6 +252,16 @@ const std::string& CoupledBoundaryCondition::get_block_name() const
     return block_name_;
 }
 
+const std::string& CoupledBoundaryCondition::get_oned_input_file() const
+{
+    return oned_input_file_;
+}
+
+void CoupledBoundaryCondition::set_oned_input_file(const std::string& path)
+{
+    oned_input_file_ = path;
+}
+
 void CoupledBoundaryCondition::set_solution_ids(int flow_id, int pressure_id, double in_out_sign)
 {
     flow_sol_id_ = flow_id;
@@ -401,6 +455,14 @@ void CoupledBoundaryCondition::distribute(const ComMod& com_mod, const CmMod& cm
     // Distribute block name
     cm.bcast(cm_mod, block_name_);
 
+    // Distribute 1D input file path
+    cm.bcast(cm_mod, oned_input_file_);
+
+    // Distribute 1D ramp and relaxation parameters
+    cm.bcast(cm_mod, &oned_ramp_steps_);
+    cm.bcast(cm_mod, &oned_ramp_ref_pressure_);
+    cm.bcast(cm_mod, &oned_relax_factor_);
+
     // Distribute face name
     cm.bcast(cm_mod, face_name_);
 
@@ -875,8 +937,27 @@ CappingSurface::CappingSurface(const CappingSurface& other)
     , valM_(other.valM_)
     , normals_(other.normals_)
 {
-    if (other.face_) {
-        face_ = std::make_unique<faceType>(*other.face_);
+    if (other.face_ != nullptr) {
+        try {
+            face_ = std::make_unique<faceType>(*other.face_);
+        } catch (const std::exception& e) {
+            throw CappingSurfaceCopyException("[CappingSurface::copy constructor] Failed to copy face_: " +
+                                             std::string(e.what()));
+        }
+        if (face_ != nullptr) {
+            if (face_->nNo > 0 && face_->gN.size() != face_->nNo) {
+                throw CappingSurfaceCopyException("[CappingSurface::copy constructor] Invalid face_: gN.size()=" +
+                                                std::to_string(face_->gN.size()) + " != nNo=" +
+                                                std::to_string(face_->nNo));
+            }
+            if (face_->nEl > 0 && face_->IEN.ncols() != face_->nEl) {
+                throw CappingSurfaceCopyException("[CappingSurface::copy constructor] Invalid face_: IEN.ncols()=" +
+                                                std::to_string(face_->IEN.ncols()) + " != nEl=" +
+                                                std::to_string(face_->nEl));
+            }
+        }
+    } else {
+        face_.reset();
     }
 }
 
@@ -886,8 +967,13 @@ CappingSurface& CappingSurface::operator=(const CappingSurface& other)
         global_node_ids_ = other.global_node_ids_;
         valM_ = other.valM_;
         normals_ = other.normals_;
-        if (other.face_) {
-            face_ = std::make_unique<faceType>(*other.face_);
+        if (other.face_ != nullptr) {
+            try {
+                face_ = std::make_unique<faceType>(*other.face_);
+            } catch (const std::exception& e) {
+                throw CappingSurfaceCopyException("[CappingSurface::operator=] Failed to copy face_: " +
+                                                  std::string(e.what()));
+            }
         } else {
             face_.reset();
         }
@@ -986,7 +1072,7 @@ void CappingSurface::init_cap_face_quadrature(const ComMod& com_mod)
 
     try {
         if (nsd != cap_nsd_) {
-            throw CappingSurfaceBaseException("[CappingSurface::init_cap_face_quadrature] Cap surface requires nsd=3.");
+            throw CappingSurfaceGeometryException("[CappingSurface::init_cap_face_quadrature] Cap surface requires nsd=3.");
         }
         face_->nG = 1;
 
@@ -1034,6 +1120,10 @@ Array<double> CappingSurface::update_element_position_global(int e, consts::Mech
 {
     using namespace consts;
 
+    if (mesh_x.nrows() < cap_nsd_ || mesh_Do.nrows() < cap_nsd_ || mesh_Dn.nrows() < cap_nsd_) {
+        throw CappingSurfaceGeometryException(
+            "[CappingSurface::update_element_position_global] Mesh arrays must have at least 3 rows.");
+    }
     Array<double> xl(cap_nsd_, face_->eNoN);
 
     for (int a = 0; a < face_->eNoN; a++) {
@@ -1063,6 +1153,13 @@ Array<double> CappingSurface::update_element_position_global(int e, consts::Mech
 /// @return The Jacobian and normal vector.
 std::pair<double, Vector<double>> CappingSurface::compute_jacobian_and_normal(const Array<double>& xl, int e, int g) const
 {
+    if (xl.nrows() != cap_nsd_ || xl.ncols() != face_->eNoN) {
+        throw CappingSurfaceGeometryException("[CappingSurface::compute_jacobian_and_normal] xl has wrong dimensions: " +
+                                              std::to_string(xl.nrows()) + "x" + std::to_string(xl.ncols()) +
+                                              " (expected " + std::to_string(cap_nsd_) + "x" +
+                                              std::to_string(face_->eNoN) + ").");
+    }
+
     // Get the shape function derivatives for the Gauss point.
     Array<double> Nx_g = face_->Nx.rslice(g);
     Array<double> xXi(cap_nsd_, cap_insd_);
@@ -1084,14 +1181,19 @@ std::pair<double, Vector<double>> CappingSurface::compute_jacobian_and_normal(co
     Jac = utils::norm(n);
 
     if (utils::is_zero(Jac)) {
-        throw CappingSurfaceBaseException("[CappingSurface::compute_jacobian_and_normal] Zero Jacobian at Gauss point " +
+        throw CappingSurfaceGeometryException("[CappingSurface::compute_jacobian_and_normal] Zero Jacobian at Gauss point " +
                                               std::to_string(g));
     }
 
     n = n / Jac;
 
     // Check if the initial normals are provided and if they are valid.
     if (normals_.ncols() > 0 && normals_.nrows() == cap_nsd_) {
+        if (e < 0 || e >= normals_.ncols()) {
+            throw CappingSurfaceGeometryException("[CappingSurface::compute_jacobian_and_normal] Element index e=" +
+                                                std::to_string(e) + " is out of bounds for normals_ (ncols=" +
+                                                std::to_string(normals_.ncols()) + ").");
+        }
 
         Vector<double> n0(cap_nsd_);
         for (int i = 0; i < cap_nsd_; i++) {
@@ -1175,7 +1277,19 @@ void CappingSurface::compute_valM(consts::MechanicalConfigurationType cfg, const
             auto [Jac, n] = compute_jacobian_and_normal(xl, e, g);
             for (int a = 0; a < face_->eNoN; a++) {
                 int gnNo_idx = face_->IEN(a, e);
-                int cap_a = gnNo_to_cap_local.at(gnNo_idx);
+                auto it = gnNo_to_cap_local.find(gnNo_idx);
+                if (it == gnNo_to_cap_local.end()) {
+                    throw CappingSurfaceAssemblyException("[CappingSurface::compute_valM] IEN entry (element " +
+                                                          std::to_string(e) + ", node " + std::to_string(a) +
+                                                          ") contains invalid gnNo index " + std::to_string(gnNo_idx) +
+                                                          " not found in cap face nodes.");
+                }
+                int cap_a = it->second;
+                if (cap_a < 0 || cap_a >= cap_nNo) {
+                    throw CappingSurfaceAssemblyException("[CappingSurface::compute_valM] Invalid cap face-local index cap_a=" +
+                                                          std::to_string(cap_a) + " (cap_nNo=" + std::to_string(cap_nNo) +
+                                                          ")");
+                }
                 for (int i = 0; i < cap_nsd_; i++) {
                     valM_(i, cap_a) += face_->N(a, g) * face_->w(g) * Jac * n(i);
                 }
@@ -1311,3 +1425,23 @@ void CoupledBoundaryCondition::bcast_coupled_neumann_pressure(const CmMod& cm_mo
     cm.bcast(cm_mod, &pr);
     set_pressure(pr);
 }
+
+void CoupledBoundaryCondition::bcast_coupled_dir_flowrate(const CmMod& cm_mod, cmType& cm)
+{
+    if (cm.seq()) {
+        return;
+    }
+    using namespace consts;
+    if (get_bc_type() != BoundaryConditionType::bType_Dir) {
+        return;
+    }
+    double Qo = 0.0;
+    double Qn = 0.0;
+    if (cm.mas(cm_mod)) {
+        Qo = get_Qo();
+        Qn = get_Qn();
+    }
+    cm.bcast(cm_mod, &Qo);
+    cm.bcast(cm_mod, &Qn);
+    set_flowrates(Qo, Qn);
+}