Merge pull request #2132 from TomFischer/BalanceForHTProcess

Balance for HT process

Applications/ApplicationsLib/ProjectData.cpp

@@ -448,7 +448,7 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
             process = ProcessLib::HT::createHTProcess(
                 *_mesh_vec[0], std::move(jacobian_assembler),
                 _process_variables, _parameters, integration_order,
-                process_config);
+                process_config, project_directory, output_directory);
         }
         else
 #endif

ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.cpp

@@ -55,13 +55,14 @@ CalculateSurfaceFlux::CalculateSurfaceFlux(
 
 void CalculateSurfaceFlux::integrate(GlobalVector const& x,
                                      MeshLib::PropertyVector<double>& balance,
+                                     double const t,
                                      Process const& bulk_process)
 {
     DBUG("Integrate CalculateSurfaceFlux.");
 
     GlobalExecutor::executeMemberOnDereferenced(
         &CalculateSurfaceFluxLocalAssemblerInterface::integrate,
-        _local_assemblers, x, balance, bulk_process);
+        _local_assemblers, x, balance, t, bulk_process);
 }
 
 }  // namespace ProcessLib

ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.h

@@ -29,10 +29,12 @@ class CalculateSurfaceFlux final
     /// Executes for each element of the mesh the local intergration procedure.
     /// @param x The global solution the intergration values will be fetched of.
     /// @param balance The vector the integration results will be stored in.
+    /// @param t The balance will be computed at the time t.
     /// @param bulk_process Stores the variable that is used for the
     /// integration.
     void integrate(GlobalVector const& x,
                    MeshLib::PropertyVector<double>& balance,
+                   double const t,
                    Process const& bulk_process);
 
 private:

ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFluxLocalAssembler.h

@@ -31,6 +31,7 @@ class CalculateSurfaceFluxLocalAssemblerInterface
     virtual void integrate(std::size_t element_id,
                            GlobalVector const& x,
                            MeshLib::PropertyVector<double>& balance,
+                           double const t,
                            Process const& bulk_process) = 0;
 };
 
@@ -101,11 +102,13 @@ class CalculateSurfaceFluxLocalAssembler final
     /// @param balance PropertyVector the integration result will be stored
     /// into, where balance has the same number of entries like mesh elements
     /// exists.
+    /// @param t The integration is performed at the time t.
     /// @param bulk_process Reference to the bulk process is used to compute the
     /// flux.
     void integrate(std::size_t element_id,
                    GlobalVector const& x,
                    MeshLib::PropertyVector<double>& balance,
+                   double const t,
                    Process const& bulk_process) override
     {
         auto surface_element_normal =
@@ -118,16 +121,16 @@ class CalculateSurfaceFluxLocalAssembler final
 
         std::size_t const n_integration_points =
             _integration_method.getNumberOfPoints();
-        // balance[id_of_element] =
+        // balance[id_of_element] +=
         //   int_{\Gamma_e} \alpha \cdot flux \cdot normal \d \Gamma
         for (std::size_t ip(0); ip < n_integration_points; ip++)
         {
             auto const& wp = _integration_method.getWeightedPoint(ip);
 
             auto const bulk_element_point = MeshLib::getBulkElementPoint(
                 bulk_process.getMesh(), _bulk_element_id, _bulk_face_id, wp);
-            auto const bulk_flux =
-                bulk_process.getFlux(_bulk_element_id, bulk_element_point, x);
+            auto const bulk_flux = bulk_process.getFlux(
+                _bulk_element_id, bulk_element_point, t, x);
 
             for (int component_id(0);
                  component_id < balance.getNumberOfComponents();

ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h

@@ -104,10 +104,9 @@ class LocalAssemblerData : public GroundwaterFlowLocalAssemblerInterface
 
     /// Computes the flux in the point \c p_local_coords that is given in local
     /// coordinates using the values from \c local_x.
-    // TODO add time dependency
-    std::vector<double> getFlux(
-        MathLib::Point3d const& p_local_coords,
-        std::vector<double> const& local_x) const override
+    Eigen::Vector3d getFlux(MathLib::Point3d const& p_local_coords,
+                            double const t,
+                            std::vector<double> const& local_x) const override
     {
         // eval dNdx and invJ at p
         using FemType =
@@ -123,19 +122,16 @@ class LocalAssemblerData : public GroundwaterFlowLocalAssemblerInterface
         // here, which is not affected by axial symmetry.
         fe.computeShapeFunctions(p_local_coords.getCoords(), shape_matrices,
                                  GlobalDim, false);
-        std::vector<double> flux;
-        flux.resize(3);
 
         // fetch hydraulic conductivity
         SpatialPosition pos;
         pos.setElementID(_element.getID());
-        // TODO remove follwing line if time dependency is implemented
-        double const t = 0.0;
         auto const k = _process_data.hydraulic_conductivity(t, pos)[0];
 
-        Eigen::Map<Eigen::RowVectorXd>(flux.data(), flux.size()) =
+        Eigen::Vector3d flux;
+        flux.head<GlobalDim>() =
             -k * shape_matrices.dNdx *
-            Eigen::Map<const Eigen::VectorXd>(local_x.data(), local_x.size());
+            Eigen::Map<const NodalVectorType>(local_x.data(), local_x.size());
 
         return flux;
     }

ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h

@@ -47,20 +47,23 @@ class GroundwaterFlowProcess final : public Process
     bool isLinear() const override { return true; }
     //! @}
 
-    std::vector<double> getFlux(std::size_t element_id,
-                                MathLib::Point3d const& p,
-                                GlobalVector const& x) const override
+    Eigen::Vector3d getFlux(std::size_t element_id,
+                            MathLib::Point3d const& p,
+                            double const t,
+                            GlobalVector const& x) const override
     {
         // fetch local_x from primary variable
         std::vector<GlobalIndexType> indices_cache;
         auto const r_c_indices = NumLib::getRowColumnIndices(
             element_id, *_local_to_global_index_map, indices_cache);
         std::vector<double> local_x(x.get(r_c_indices.rows));
 
-        return _local_assemblers[element_id]->getFlux(p, local_x);
+        return _local_assemblers[element_id]->getFlux(p, t, local_x);
     }
 
     void postTimestepConcreteProcess(GlobalVector const& x,
+                                     const double t,
+                                     const double /*delta_t*/,
                                      int const process_id) override
     {
         //For this single process, process_id is always zero.
@@ -87,7 +90,7 @@ class GroundwaterFlowProcess final : public Process
                 _balance_mesh->getProperties()
                     .template getPropertyVector<double>(_balance_pv_name);
 
-            balance.integrate(x, *balance_pv, *this);
+            balance.integrate(x, *balance_pv, t, *this);
             // post: surface_mesh has vectorial element property
 
             // TODO output, if output classes are ready this has to be

ProcessLib/HT/CreateHTProcess.cpp

@@ -12,9 +12,12 @@
 #include "MaterialLib/Fluid/FluidProperties/CreateFluidProperties.h"
 #include "MaterialLib/PorousMedium/CreatePorousMediaProperties.h"
 
+#include "MeshLib/IO/readMeshFromFile.h"
+
 #include "ProcessLib/Output/CreateSecondaryVariables.h"
 #include "ProcessLib/Parameter/ConstantParameter.h"
 #include "ProcessLib/Utils/ProcessUtils.h"
+#include "ProcessLib/CalculateSurfaceFlux/ParseCalculateSurfaceFluxData.h"
 
 #include "HTProcess.h"
 #include "HTMaterialProperties.h"
@@ -30,7 +33,9 @@ std::unique_ptr<Process> createHTProcess(
     std::vector<ProcessVariable> const& variables,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     unsigned const integration_order,
-    BaseLib::ConfigTree const& config)
+    BaseLib::ConfigTree const& config,
+    std::string const& project_directory,
+    std::string const& output_directory)
 {
     //! \ogs_file_param{prj__processes__process__type}
     config.checkConfigParameter("type", "HT");
@@ -194,6 +199,33 @@ std::unique_ptr<Process> createHTProcess(
         DBUG("Use \'%s\' as Biot's constant.", biot_constant->name.c_str());
     }
 
+    // for the balance
+    std::string mesh_name; // surface mesh the balance will computed on
+    std::string balance_pv_name;
+    std::string balance_out_fname;
+    std::unique_ptr<MeshLib::Mesh> surface_mesh;
+    ProcessLib::parseCalculateSurfaceFluxData(
+        config, mesh_name, balance_pv_name, balance_out_fname);
+
+    if (!mesh_name.empty())  // balance is optional
+    {
+        mesh_name = BaseLib::copyPathToFileName(mesh_name, project_directory);
+
+        balance_out_fname =
+            BaseLib::copyPathToFileName(balance_out_fname, output_directory);
+
+        surface_mesh.reset(MeshLib::IO::readMeshFromFile(mesh_name));
+
+        DBUG(
+            "read balance meta data:\n\tbalance mesh:\"%s\"\n\tproperty name: "
+            "\"%s\"\n\toutput to: \"%s\"",
+            mesh_name.c_str(), balance_pv_name.c_str(),
+            balance_out_fname.c_str());
+
+        // Surface mesh and bulk mesh must have equal axial symmetry flags!
+        surface_mesh->setAxiallySymmetric(mesh.isAxiallySymmetric());
+    }
+
     std::unique_ptr<HTMaterialProperties> material_properties =
         std::make_unique<HTMaterialProperties>(
             std::move(porous_media_properties),
@@ -223,7 +255,8 @@ std::unique_ptr<Process> createHTProcess(
         mesh, std::move(jacobian_assembler), parameters, integration_order,
         std::move(process_variables), std::move(material_properties),
         std::move(secondary_variables), std::move(named_function_caller),
-        use_monolithic_scheme);
+        use_monolithic_scheme, std::move(surface_mesh),
+        std::move(balance_pv_name), std::move(balance_out_fname));
 }
 
 }  // namespace HT

ProcessLib/HT/CreateHTProcess.h

@@ -22,7 +22,9 @@ std::unique_ptr<Process> createHTProcess(
     std::vector<ProcessVariable> const& variables,
     std::vector<std::unique_ptr<ParameterBase>> const& parameters,
     unsigned const integration_order,
-    BaseLib::ConfigTree const& config);
+    BaseLib::ConfigTree const& config,
+    std::string const& project_directory,
+    std::string const& output_directory);
 
 }  // namespace HT
 }  // namespace ProcessLib

ProcessLib/HT/HTFEM.h

@@ -88,6 +88,71 @@ class HTFEM : public HTLocalAssemblerInterface
         return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
     }
 
+    /// Computes the flux in the point \c pnt_local_coords that is given in
+    /// local coordinates using the values from \c local_x.
+    Eigen::Vector3d getFlux(
+        MathLib::Point3d const& pnt_local_coords,
+        double const t,
+        std::vector<double> const& local_x) const override
+    {
+        // eval dNdx and invJ at given point
+        using FemType =
+            NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
+
+        FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
+            &this->_element));
+
+        typename ShapeMatricesType::ShapeMatrices shape_matrices(
+            ShapeFunction::DIM, GlobalDim, ShapeFunction::NPOINTS);
+
+        // Note: Axial symmetry is set to false here, because we only need dNdx
+        // here, which is not affected by axial symmetry.
+        fe.computeShapeFunctions(pnt_local_coords.getCoords(), shape_matrices,
+                                 GlobalDim, false);
+
+        // fetch permeability, viscosity, density
+        SpatialPosition pos;
+        pos.setElementID(this->_element.getID());
+
+        MaterialLib::Fluid::FluidProperty::ArrayType vars;
+
+        // local_x contains the local temperature and pressure values
+        NumLib::shapeFunctionInterpolate(
+            local_x, shape_matrices.N,
+            vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)],
+            vars[static_cast<int>(
+                MaterialLib::Fluid::PropertyVariableType::p)]);
+
+        auto const K =
+            this->_material_properties.porous_media_properties
+                .getIntrinsicPermeability(t, pos)
+                .getValue(t, pos, 0.0,
+                          vars[static_cast<int>(
+                              MaterialLib::Fluid::PropertyVariableType::T)]);
+
+        auto const mu = this->_material_properties.fluid_properties->getValue(
+            MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
+        GlobalDimMatrixType const K_over_mu = K / mu;
+
+        auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
+            &local_x[local_x.size()/2], ShapeFunction::NPOINTS);
+        GlobalDimVectorType q =
+            -K_over_mu * shape_matrices.dNdx * p_nodal_values;
+
+        if (this->_material_properties.has_gravity)
+        {
+            auto const rho_w =
+                this->_material_properties.fluid_properties->getValue(
+                    MaterialLib::Fluid::FluidPropertyType::Density, vars);
+            auto const b = this->_material_properties.specific_body_force;
+            q += K_over_mu * rho_w * b;
+        }
+
+        Eigen::Vector3d flux;
+        flux.head<GlobalDim>() = q;
+        return flux;
+    }
+
 protected:
     MeshLib::Element const& _element;
     HTMaterialProperties const& _material_properties;

ProcessLib/HT/HTLocalAssemblerInterface.h

@@ -59,14 +59,20 @@ class HTLocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
         NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
         std::vector<double>& /*cache*/) const = 0;
 
+    virtual Eigen::Vector3d getFlux(
+        MathLib::Point3d const& pnt_local_coords,
+        double const t,
+        std::vector<double> const& local_x) const = 0;
+
 protected:
-    // TODO: remove _coupled_solutions or move integration point data from local
-    // assembler class to a new class to make local assembler unique for each
-    // process.
-    /** Pointer to CoupledSolutionsForStaggeredScheme that is set in a member of
-     *  Process class, setCoupledTermForTheStaggeredSchemeToLocalAssemblers.
-     *  It is used for calculate the secondary variables like velocity for
-     *  coupled processes.
+    // TODO: remove _coupled_solutions or move integration point data from
+    // local assembler class to a new class to make local assembler unique
+    // for each process.
+    /** Pointer to CoupledSolutionsForStaggeredScheme that is set in a
+     * member of Process class,
+     * setCoupledTermForTheStaggeredSchemeToLocalAssemblers. It is used for
+     * calculate the secondary variables like velocity for coupled
+     * processes.
      */
     CoupledSolutionsForStaggeredScheme* _coupled_solutions;
 };