[PCS] Added the computation of the secondary variables after during the un-coupling

ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h

@@ -154,9 +154,67 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
 template <typename ShapeFunction, typename IntegrationMethod,
           unsigned GlobalDim>
 void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
-    computeSecondaryVariable(std::vector<double> const& local_x)
+    computeSecondaryVariableConcrete(std::vector<double> const& local_x)
 {
-  (void) local_x;
+    auto const local_matrix_size = local_x.size();
+    assert(local_matrix_size == ShapeFunction::NPOINTS);
+
+    const auto local_p_vec =
+        MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
+
+    unsigned const n_integration_points =
+        _integration_method.getNumberOfPoints();
+
+    SpatialPosition pos;
+    pos.setElementID(_element.getID());
+    _material_properties.setMaterialID(pos);
+    const Eigen::MatrixXd& perm =
+        _material_properties.getPermeability(t, pos, _element.getDimension());
+
+    for (unsigned ip = 0; ip < n_integration_points; ip++)
+    {
+        auto const& sm = _shape_matrices[ip];
+        double p = 0.;
+        NumLib::shapeFunctionInterpolate(local_x, sm.N, p);
+        // TODO : compute _temperature from the heat transport pcs
+
+        const double rho_g =
+            _material_properties.getLiquidDensity(p, _temperature) *
+            _gravitational_acceleration;
+        // Compute viscosity:
+        const double mu = _material_properties.getViscosity(p, _temperature);
+
+        // Assemble Laplacian, K, and RHS by the gravitational term
+        if (perm.size() == 1)  // Save time for isotropic permeability.
+        {
+            //  Use scalar number for isotropic permeability
+            //  to save the computation time.
+            const double K = perm(0, 0) / mu;
+            // Compute the velocity
+            GlobalDimVectorType darcy_velocity = -K * sm.dNdx * local_p_vec;
+            // gravity term
+            if (_gravitational_axis_id >= 0)
+                darcy_velocity[GlobalDim - 1] -= K * rho_g;
+            for (unsigned d = 0; d < GlobalDim; ++d)
+            {
+                _darcy_velocities[d][ip] = darcy_velocity[d];
+            }
+        }
+        else
+        {
+            // Compute the velocity
+            GlobalDimVectorType darcy_velocity = -perm * sm.dNdx * local_p_vec / mu;
+            if (_gravitational_axis_id >= 0)
+            {
+                darcy_velocity.noalias() -=  rho_g * perm.col(GlobalDim - 1) / mu;
+            }
+            for (unsigned d = 0; d < GlobalDim; ++d)
+            {
+                _darcy_velocities[d][ip] = darcy_velocity[d];
+            }
+        }
+    }
+
 }
 
 }  // end of namespace

ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h

@@ -85,7 +85,7 @@ public:
                   std::vector<double>& local_K_data,
                   std::vector<double>& local_b_data) override;
 
-    void computeSecondaryVariable(std::vector<double> const& local_x) override;
+    void computeSecondaryVariableConcrete(std::vector<double> const& local_x) override;
 
     Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
         const unsigned integration_point) const override

ProcessLib/LiquidFlow/LiquidFlowProcess.cpp

@@ -113,5 +113,13 @@ void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
         dx_dx, M, K, b, Jac);
 }
 
+void LiquidFlowProcess::computeSecondaryVariableConcrete(GlobalVector const& x)
+{
+    DBUG("Compute the velocity for LiquidFlowProcess.");
+    GlobalExecutor::executeMemberOnDereferenced(
+            &LiquidFlowLocalAssemblerInterface::computeSecondaryVariable,
+            _local_assemblers, *_local_to_global_index_map, x);
+}
+        
 }  // end of namespace
 }  // end of namespace

ProcessLib/LiquidFlow/LiquidFlowProcess.h

@@ -69,6 +69,8 @@ public:
         double const gravitational_acceleration,
         BaseLib::ConfigTree const& config);
 
+    void computeSecondaryVariableConcrete(GlobalVector const& x) override;
+
     bool isLinear() const override { return true; }
 private:
     void initializeConcreteProcess(

ProcessLib/LocalAssemblerInterface.cpp

@@ -26,11 +26,13 @@ void LocalAssemblerInterface::assembleWithJacobian(
         "assembler.");
 }
 
-void LocalAssemblerInterface::computeSecondaryVariable(std::vector<double> const& /*local_x*/)
+void LocalAssemblerInterface::computeSecondaryVariable(std::size_t const mesh_item_id,
+                              NumLib::LocalToGlobalIndexMap const& dof_table,
+                              GlobalVector const& x)
 {
-    OGS_FATAL(
-        "computeSecondaryVariable(...) function is not implemented in the local "
-        "assembler.");
+    auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
+    auto const local_x = x.get(indices);
+    computeSecondaryVariableConcrete(local_x);
 }
 
 void LocalAssemblerInterface::preTimestep(

ProcessLib/LocalAssemblerInterface.h

@@ -45,7 +45,9 @@ public:
                                       std::vector<double>& local_b_data,
                                       std::vector<double>& local_Jac_data);
 
-    virtual void computeSecondaryVariable(std::vector<double> const& local_x);
+    virtual void computeSecondaryVariable(std::size_t const mesh_item_id,
+                              NumLib::LocalToGlobalIndexMap const& dof_table,
+                              GlobalVector const& x);
 
     virtual void preTimestep(std::size_t const mesh_item_id,
                              NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -72,6 +74,9 @@ private:
     }
 
     virtual void postTimestepConcrete(std::vector<double> const& /*local_x*/) {}
+
+    virtual void computeSecondaryVariableConcrete
+                         (std::vector<double> const& /*local_x*/) {}
 };
 
 } // namespace ProcessLib

ProcessLib/Process.cpp

@@ -242,6 +242,11 @@ void Process::postTimestep(GlobalVector const& x)
     postTimestepConcreteProcess(x);
 }
 
+void Process::computeSecondaryVariable(GlobalVector const& x)
+{
+    computeSecondaryVariableConcrete(x);
+}
+
 void Process::preIteration(const unsigned iter, const GlobalVector &x)
 {
     // In every new iteration cached values of secondary variables are expired.

ProcessLib/Process.h

@@ -59,6 +59,9 @@ public:
     void preIteration(const unsigned iter,
                       GlobalVector const& x) override final;
 
+    /// computeSecondaryVariable for the coupled equations or for output.
+    void computeSecondaryVariable(GlobalVector const& x);
+
     NumLib::IterationResult postIteration(GlobalVector const& x) override final;
 
     void initialize();
@@ -148,6 +151,8 @@ private:
     virtual void preIterationConcreteProcess(const unsigned /*iter*/,
                                              GlobalVector const& /*x*/){}
 
+    virtual void computeSecondaryVariableConcrete(GlobalVector const& /*x*/) {};
+
     virtual NumLib::IterationResult postIterationConcreteProcess(
         GlobalVector const& /*x*/)
     {

ProcessLib/UncoupledProcessesTimeLoop.cpp

@@ -505,6 +505,8 @@ bool UncoupledProcessesTimeLoop::loop()
             nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
                 x, timestep, t, delta_t, *spd, *_output);
 
+            pcs.computeSecondaryVariable(x);
+
             INFO("[time] Solving process #%u took %g s in timestep #%u.",
                  timestep, time_timestep.elapsed(), pcs_idx);