Merge branch 'AlgebraicBcSettingForLinearSolver' into 'master'

Add Algebraic BC in the assemble of the HeatTransportBHE Process

See merge request ogs/ogs!5025

Documentation/ProjectFile/prj/processes/process/HEAT_TRANSPORT_BHE/t_linear.md

@@ -0,0 +1 @@
+The default setting is false. When activated, non-linear iterations will be turned off. The simulation will just go through a single iteration within a time step. If it's set to true, the **\<use_algebraic_bc\>** option must be also set to true. Only when fluid properties does not change much with temperature, this feature can be activated. The effect of skipping the non-linear iteration is much faster simulation speed.

Documentation/ProjectFile/prj/processes/process/HEAT_TRANSPORT_BHE/t_use_algebraic_bc.md

@@ -0,0 +1 @@
+The default setting is false. When algebraic boundary condition is activated, the inflow temperatures at the top and bottom of the BHE are no longer imposed as a Dirichlet type boundary condition. Instead, the linear equation system will be enlarged to accommodate the numerical constrains. The effect of this feature is to accelerate BHE simulations.

Documentation/ProjectFile/prj/processes/process/HEAT_TRANSPORT_BHE/t_weighting_factor.md

@@ -0,0 +1 @@
+This value is the weighting factor used to impose an algebraic boundary condition in the **HeatTransportBHE** process. The default value is 100. This value is only effective, when the option **\<use_algebraic_bc\>** is set to true. This value is multiplied with the max value of the column wise inner product and then imposed in the additional rows at the end of the linear equation system, so that the unknown temperatures at the top and bottom of the BHE are imposed with certain desired values.

MathLib/LinAlg/LinAlg.cpp

@@ -167,6 +167,34 @@ void matMultAdd(PETScMatrix const& A, PETScVector const& v1,
                v3.getRawVector());
 }
 
+void linearSysNormalize(PETScMatrix const& /*A*/, PETScMatrix& /*new_A*/,
+                        PETScVector const& /*b*/, PETScVector& /*new_b*/)
+{
+    // The following block is deactivated, because there is no tests yet for the
+    // normalization operation in PETSc. This will be a task for later.
+    /*
+    assert(&A != &new_A);
+    assert(&b != &new_b);
+
+    PetscInt n_rows(0);
+    PetscInt n_cols(0);
+    MatGetSize(A.getRawMatrix(), &n_rows, &n_cols);
+    // only when A matrix is not square
+    if (n_rows != n_cols)
+    {
+        // new_b = A^T * b
+        MatMultTranspose(A.getRawMatrix(), b.getRawVector(),
+                         new_b.getRawVector());
+        // new_A = A^T * A
+        MatTranspose(A.getRawMatrix(), MAT_INITIAL_MATRIX,
+                     &(new_A.getRawMatrix()));
+    }
+    */
+    OGS_FATAL(
+        "Normalization operation is not implemented yet for PETSc library! "
+        "Program terminated.");
+}
+
 void finalizeAssembly(PETScMatrix& A)
 {
     A.finalizeAssembly(MAT_FINAL_ASSEMBLY);
@@ -313,6 +341,25 @@ void matMultAdd(EigenMatrix const& A, EigenVector const& v1,
         v2.getRawVector() + A.getRawMatrix() * v1.getRawVector();
 }
 
+void linearSysNormalize(EigenMatrix const& A, EigenMatrix& new_A,
+                        EigenVector const& b, EigenVector& new_b)
+{
+    // make sure that new_A and new_b are not the same memory
+    assert(&A != &new_A);
+    assert(&b != &new_b);
+
+    if (A.getRawMatrix().rows() == A.getRawMatrix().cols())
+    {
+        WARN(
+            "The number of rows and columns are the same for the LHS matrix."
+            "Are you sure you still need to normalize the LHS matrix and RHS "
+            "vector? ");
+    }
+
+    new_b.getRawVector() = A.getRawMatrix().transpose() * b.getRawVector();
+    new_A.getRawMatrix() = A.getRawMatrix().transpose() * A.getRawMatrix();
+}
+
 void finalizeAssembly(EigenMatrix& x)
 {
     x.getRawMatrix().makeCompressed();

MathLib/LinAlg/LinAlg.h

@@ -196,6 +196,11 @@ void matMult(PETScMatrix const& A, PETScVector const& x, PETScVector& y);
 void matMultAdd(PETScMatrix const& A, PETScVector const& v1,
                 PETScVector const& v2, PETScVector& v3);
 
+// new_A = A^T * A
+// new_b = A^T * b
+void linearSysNormalize(PETScMatrix const& A, PETScMatrix& new_A,
+                        PETScVector const& b, PETScVector& new_b);
+
 void finalizeAssembly(PETScMatrix& A);
 void finalizeAssembly(PETScVector& x);
 
@@ -264,7 +269,10 @@ void matMult(EigenMatrix const& A, EigenVector const& x, EigenVector& y);
 // v3 = A*v1 + v2
 void matMultAdd(EigenMatrix const& A, EigenVector const& v1,
                 EigenVector const& v2, EigenVector& v3);
-
+// new_A = A^T * A
+// new_b = A^T * b
+void linearSysNormalize(EigenMatrix const& A, EigenMatrix& new_A,
+                        EigenVector const& b, EigenVector& new_b);
 void finalizeAssembly(EigenMatrix& x);
 void finalizeAssembly(EigenVector& A);
 

NumLib/ODESolver/EquationSystem.h

@@ -39,6 +39,15 @@ class EquationSystem
      */
     virtual bool isLinear() const = 0;
 
+    /*! Check whether normalization of A and rhs is required.
+     *
+     * \remark
+     * In some processes, a normalization operation is required, to calculate
+     * A^T * A, and overwrite A; also calculate A^T * rhs and overwrite rhs.
+     * This parameter reflect whether such operation is required.
+     */
+    virtual bool requiresNormalization() const = 0;
+
     /*! Prepares a new iteration in the solution process of this equation.
      *
      * \param iter the current iteration number, starting from 1.

NumLib/ODESolver/MatrixTranslator.cpp

@@ -43,6 +43,25 @@ void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
     NumLib::GlobalVectorProvider::provider.releaseVector(tmp);
 }
 
+void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
+    normalizeAandRhs(GlobalMatrix& A, GlobalVector& b) const
+{
+    namespace LinAlg = MathLib::LinAlg;
+
+    // check whether A is square?
+
+    GlobalMatrix new_A(A);
+    GlobalVector new_b(b);
+    LinAlg::copy(A, new_A);
+    LinAlg::copy(b, new_b);
+    // rhs = A^T * rhs
+    // A = A^T * A
+    LinAlg::linearSysNormalize(A, new_A, b, new_b);
+
+    LinAlg::copy(new_A, A);
+    LinAlg::copy(new_b, b);
+}
+
 void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
     computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
                     GlobalVector const& b, double const dt,

NumLib/ODESolver/MatrixTranslator.h

@@ -49,6 +49,10 @@ class MatrixTranslator<ODESystemTag::FirstOrderImplicitQuasilinear>
                             const GlobalVector& b, const GlobalVector& x_prev,
                             GlobalVector& rhs) const = 0;
 
+    //! Computes \f$ A = A^T \cdot A \f$, and
+    //! also \f$ rhs = A^T \cdot rhs \f$.
+    virtual void normalizeAandRhs(GlobalMatrix& A, GlobalVector& b) const = 0;
+
     /*! Computes \c res from \c M, \c K, \c b, \f$ \hat x \f$ and \f$ x_N \f$.
      * You might also want read the remarks on
      * \ref concept_time_discretization "time discretization".
@@ -105,6 +109,10 @@ class MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>
                     const GlobalVector& b, const GlobalVector& x_prev,
                     GlobalVector& rhs) const override;
 
+    //! Computes \f$ A = A^T \cdot A \f$, and
+    //! also \f$ rhs = A^T \cdot rhs \f$.
+    void normalizeAandRhs(GlobalMatrix& A, GlobalVector& b) const override;
+
     //! Computes \f$ r = M \cdot \hat x + K \cdot x_C - b \f$.
     void computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
                          GlobalVector const& b, double dt,

NumLib/ODESolver/NonlinearSolver.cpp

@@ -167,6 +167,11 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
         sys.assemble(x_new, x_prev, process_id);
         sys.getA(A);
         sys.getRhs(*x_prev[process_id], rhs);
+
+        // Normalize the linear equation system, if required
+        if (sys.requiresNormalization())
+            sys.getAandRhsNormalized(A, rhs);
+
         INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
 
         // Subtract non-equilibrium initial residuum if set

NumLib/ODESolver/NonlinearSystem.h

@@ -121,6 +121,12 @@ class NonlinearSystem<NonlinearSolverTag::Picard> : public EquationSystem
     virtual void getRhs(GlobalVector const& x_prev,
                         GlobalVector& rhs) const = 0;
 
+    //! Writes the A_transposed times A into \c A
+    //! and also writes A_transposed times rhs into \c rhs
+    //! \pre getA() and getRhs must have been called before.
+    virtual void getAandRhsNormalized(GlobalMatrix& A,
+                                      GlobalVector& rhs) const = 0;
+
     //! Pre-compute known solutions and possibly store them internally.
     virtual void computeKnownSolutions(GlobalVector const& x,
                                        int const process_id) = 0;

NumLib/ODESolver/TimeDiscretizedODESystem.h

@@ -120,6 +120,11 @@ class TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
 
     bool isLinear() const override { return _ode.isLinear(); }
 
+    bool requiresNormalization() const override
+    {
+        return _ode.requiresNormalization();
+    }
+
     void preIteration(const unsigned iter, GlobalVector const& x) override
     {
         _ode.preIteration(iter, x);
@@ -207,6 +212,11 @@ class TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
         _mat_trans->computeRhs(*_M, *_K, *_b, x_prev, rhs);
     }
 
+    void getAandRhsNormalized(GlobalMatrix& A, GlobalVector& rhs) const override
+    {
+        _mat_trans->normalizeAandRhs(A, rhs);
+    }
+
     void computeKnownSolutions(GlobalVector const& x,
                                int const process_id) override;
 
@@ -217,6 +227,11 @@ class TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
 
     bool isLinear() const override { return _ode.isLinear(); }
 
+    bool requiresNormalization() const override
+    {
+        return _ode.requiresNormalization();
+    }
+
     void preIteration(const unsigned iter, GlobalVector const& x) override
     {
         _ode.preIteration(iter, x);

ProcessLib/HeatTransportBHE/BHE/BHECommon.h

@@ -43,6 +43,11 @@ struct BHECommon
     GroutParameters const grout;
     FlowAndTemperatureControl const flowAndTemperatureControl;
     bool const use_python_bcs;
+    constexpr bool isPowerBC() const
+    {
+        return std::visit([](auto const& ftc) { return ftc.is_power_bc; },
+                          flowAndTemperatureControl);
+    }
 };
 }  // end of namespace BHE
 }  // end of namespace HeatTransportBHE

ProcessLib/HeatTransportBHE/BHE/CreateFlowAndTemperatureControl.cpp

@@ -67,6 +67,7 @@ FlowAndTemperatureControl createFlowAndTemperatureControl(
 
         //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__flow_and_temperature_control__FixedPowerConstantFlow__flow_rate}
         auto const flow_rate = config.getConfigParameter<double>("flow_rate");
+
         return FixedPowerConstantFlow{flow_rate, power,
                                       refrigerant.specific_heat_capacity,
                                       refrigerant.density};

ProcessLib/HeatTransportBHE/BHE/FlowAndTemperatureControl.h

@@ -35,6 +35,7 @@ struct TemperatureCurveConstantFlow
     }
     double flow_rate;
     MathLib::PiecewiseLinearInterpolation const& temperature_curve;
+    static constexpr bool is_power_bc = false;
 };
 
 struct TemperatureCurveFlowCurve
@@ -47,6 +48,7 @@ struct TemperatureCurveFlowCurve
     }
     MathLib::PiecewiseLinearInterpolation const& flow_rate_curve;
     MathLib::PiecewiseLinearInterpolation const& temperature_curve;
+    static constexpr bool is_power_bc = false;
 };
 
 struct FixedPowerConstantFlow
@@ -60,6 +62,7 @@ struct FixedPowerConstantFlow
     double power;  // Value is expected to be in Watt.
     double heat_capacity;
     double density;
+    static constexpr bool is_power_bc = true;
 };
 
 struct FixedPowerFlowCurve
@@ -74,6 +77,7 @@ struct FixedPowerFlowCurve
     double power;  // Value is expected to be in Watt.
     double heat_capacity;
     double density;
+    static constexpr bool is_power_bc = true;
 };
 
 struct PowerCurveConstantFlow
@@ -92,6 +96,7 @@ struct PowerCurveConstantFlow
     double flow_rate;
     double heat_capacity;
     double density;
+    static constexpr bool is_power_bc = true;
 };
 
 struct PowerCurveFlowCurve
@@ -111,6 +116,7 @@ struct PowerCurveFlowCurve
 
     double heat_capacity;
     double density;
+    static constexpr bool is_power_bc = true;
 };
 
 struct BuildingPowerCurveConstantFlow
@@ -134,6 +140,7 @@ struct BuildingPowerCurveConstantFlow
     double flow_rate;
     double heat_capacity;
     double density;
+    static constexpr bool is_power_bc = true;
 };
 
 using FlowAndTemperatureControl = std::variant<TemperatureCurveConstantFlow,

ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp

@@ -104,6 +104,39 @@ std::unique_ptr<Process> createHeatTransportBHEProcess(
         //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__use_server_communication}
         config.getConfigParameter<bool>("use_server_communication", false);
 
+    auto const using_algebraic_bc =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__use_algebraic_bc}
+        config.getConfigParameter<bool>("use_algebraic_bc", false);
+
+    auto const weighting_factor =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__weighting_factor}
+        config.getConfigParameter<float>("weighting_factor", 100.0);
+
+    auto const is_linear =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__linear}
+        config.getConfigParameter<bool>("linear", false);
+    if (is_linear)
+    {
+        if (!using_algebraic_bc)
+        {
+            OGS_FATAL(
+                "You specified that the process simulated by OGS is linear. "
+                "For the Heat-Transport-BHE process this can only be done "
+                "together with setting the use_algebraic_bc option to true.")
+        }
+        else
+        {
+            WARN(
+                "You specified that the process simulated by OGS is linear. "
+                "With that optimization the process will be assembled only "
+                "once and the non-linear solver will do only one iteration per "
+                "time step. No non-linearities will be resolved and OGS will "
+                "not detect if there are any non-linearities. It is your "
+                "responsibility to ensure that the assembled equation systems "
+                "are linear, indeed! There is no safety net!");
+        }
+    }
+
     for (
         auto const& bhe_config :
         //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger}
@@ -218,7 +251,8 @@ std::unique_ptr<Process> createHeatTransportBHEProcess(
 
     HeatTransportBHEProcessData process_data(
         std::move(media_map), std::move(bhes), py_object, using_tespy,
-        using_server_communication);
+        using_server_communication,
+        {using_algebraic_bc, weighting_factor, is_linear});
 
     SecondaryVariableCollection secondary_variables;