Issue 512 3D-1D coupling

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_interface.h svOneD_interface.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

@@ -190,10 +190,10 @@ rmshType::rmshType()
 
 
 /////////////////////////////////////////////////////////
-// s v Z e r o D S o l v e r I n t e r f a c e T y p e //
+// s v Z e r o D S o l v e r I n t e r f a c e D a t a //
 /////////////////////////////////////////////////////////
 
-void svZeroDSolverInterfaceType::set_data(const svZeroDSolverInterfaceParameters& params)
+void svZeroDSolverInterfaceData::set_data(const svZeroDSolverInterfaceParameters& params)
 {
   if (!params.defined()) {
     return;
@@ -217,3 +217,13 @@ void svZeroDSolverInterfaceType::set_data(const svZeroDSolverInterfaceParameters
 
   has_data = true;
 }
+
+void svOneDSolverInterfaceData::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,15 +775,21 @@ 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;
 };
 
 //----------------------------
-// svZeroDSolverInterfaceType
+// svZeroDSolverInterfaceData
 //----------------------------
 // This class stores information used to interface to
 // the svZeroDSolver.
 //
-class svZeroDSolverInterfaceType
+class svZeroDSolverInterfaceData
 {
   public:
 
@@ -807,6 +817,32 @@ class svZeroDSolverInterfaceType
     void set_data(const svZeroDSolverInterfaceParameters& params);
 };
 
+/// \brief Stores information used to interface with svOneDSolver.
+///
+/// This class stores the global svOneDSolver interface settings read from the
+/// solver XML file. Per-face 1D input files are stored separately in
+/// cplFaceType::oned_input_file.
+class svOneDSolverInterfaceData
+{
+  public:
+    /// \brief Path to the svOneDSolver interface shared library.
+    ///
+    /// This may be provided either with the platform-specific extension
+    /// (.so/.dylib) or without it.
+    std::string solver_library;
+
+    /// \brief True if the svOneDSolver interface settings were read from XML.
+    ///
+    /// This is set to true after the svOneDSolver_interface XML element has
+    /// been parsed successfully.
+    bool has_data = false;
+
+    /// \brief Read svOneDSolver interface settings from parsed parameters.
+    ///
+    /// \param params Parsed svOneDSolver interface parameters.
+    void set_data(const svOneDSolverInterfaceParameters& params);
+};
+
 /// @brief For coupled 0D-3D problems
 //
 class cplBCType
@@ -822,6 +858,9 @@ class cplBCType
     //  Whether to use svZeroD
     bool useSvZeroD = false;
 
+    //  Whether to use svOneD (svOneDSolver)
+    bool useSvOneD = false;
+
     //  Whether to initialize RCR from flow data
     bool initRCR = false;
 
@@ -851,7 +890,11 @@ class cplBCType
     std::string commuName;
     //std::string commuName = ".CPLBC_0D_3D.tmp";
 
-    svZeroDSolverInterfaceType svzerod_solver_interface;
+    /// @brief Data structure used for coupling with svZeroD code
+    svZeroDSolverInterfaceData svzerod_solver_interface;
+
+    /// @brief Data structure used for coupling with svOneD code
+    svOneDSolverInterfaceData svOneD_solver_interface;
 
     /// @brief The name of history file containing "X"
     std::string saveName;

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_);
 
@@ -986,7 +1048,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;
 
@@ -1084,7 +1146,7 @@ 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));
     }
 
@@ -1311,3 +1373,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);
+}