Merge pull request #2590 from endJunction/RMMassLumping

RichardsMechanics: mass lumping

Documentation/ProjectFile/prj/processes/process/RICHARDS_MECHANICS/t_mass_lumping.md

+A tag for enabling diagonal lumping of the mass matrix.

ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp

@@ -188,6 +188,10 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
         ProcessLib::RichardsFlow::createRichardsFlowMaterialProperties(
             flow_material_config, material_ids, parameters);
 
+    auto const mass_lumping =
+        //! \ogs_file_param{prj__processes__process__RICHARDS_MECHANICS__mass_lumping}
+        config.getConfigParameter<bool>("mass_lumping", false);
+
     RichardsMechanicsProcessData<DisplacementDim> process_data{
         material_ids,
         std::move(flow_material),
@@ -197,7 +201,8 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
         solid_density,
         solid_bulk_modulus,
         temperature,
-        specific_body_force};
+        specific_body_force,
+        mass_lumping};
 
     SecondaryVariableCollection secondary_variables;
 

ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h

@@ -281,6 +281,13 @@ void RichardsMechanicsLocalAssembler<
             .noalias() += N_p.transpose() * S_L * rho_LR * alpha *
                           identity2.transpose() * B * w;
     }
+
+    if (_process_data.apply_mass_lumping)
+    {
+        auto Mpp = M.template block<pressure_size, pressure_size>(
+            pressure_index, pressure_index);
+        Mpp = Mpp.colwise().sum().eval().asDiagonal();
+    }
 }
 
 template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
@@ -553,6 +560,11 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
             dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
     }
 
+    if (_process_data.apply_mass_lumping)
+    {
+        storage_p = storage_p.colwise().sum().eval().asDiagonal();
+    }
+
     // pressure equation, pressure part.
     local_Jac
         .template block<pressure_size, pressure_size>(pressure_index,

ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h

@@ -48,7 +48,8 @@ struct RichardsMechanicsProcessData
         ParameterLib::Parameter<double> const& solid_bulk_modulus_,
         ParameterLib::Parameter<double> const& temperature_,
         Eigen::Matrix<double, DisplacementDim, 1>
-            specific_body_force_)
+            specific_body_force_,
+        bool const apply_mass_lumping_)
         : material_ids(material_ids_),
           flow_material{std::move(flow_material_)},
           solid_materials{std::move(solid_materials_)},
@@ -57,7 +58,8 @@ struct RichardsMechanicsProcessData
           solid_density(solid_density_),
           solid_bulk_modulus(solid_bulk_modulus_),
           temperature(temperature_),
-          specific_body_force(std::move(specific_body_force_))
+          specific_body_force(std::move(specific_body_force_)),
+          apply_mass_lumping(apply_mass_lumping_)
     {
     }
 
@@ -105,6 +107,8 @@ struct RichardsMechanicsProcessData
     MeshLib::PropertyVector<double>* element_saturation = nullptr;
     MeshLib::PropertyVector<double>* pressure_interpolated = nullptr;
 
+    bool const apply_mass_lumping;
+
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
 };
 

ProcessLib/RichardsMechanics/Tests.cmake

@@ -130,6 +130,26 @@ AddTest(
     GLOB Richards_2d_geometry_OBSERVATION_POINT_pcs_0_ts_*.vtu displacement displacement 2e-14 0
     GLOB Richards_2d_geometry_OBSERVATION_POINT_pcs_0_ts_*.vtu pressure pressure 1e-7 1e-15
 )
+AddTest(
+    NAME RichardsMechanics_RichardsFlow_2d_small_masslumping
+    PATH RichardsMechanics
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS RichardsFlow_2d_small_masslumping.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS NOT OGS_USE_MPI
+    DIFF_DATA
+    GLOB RichardsFlow_2d_small_masslumping_pcs_0_ts_*.vtu displacement displacement 2e-14 0
+    GLOB RichardsFlow_2d_small_masslumping_pcs_0_ts_*.vtu sigma sigma 1e-8 0
+    GLOB RichardsFlow_2d_small_masslumping_pcs_0_ts_*.vtu epsilon epsilon 1e-15 0
+    GLOB RichardsFlow_2d_small_masslumping_pcs_0_ts_*.vtu pressure pressure 1e-7 1e-15
+    GLOB RichardsFlow_2d_small_masslumping_pcs_0_ts_*.vtu pressure_interpolated pressure_interpolated 1e-7 1e-15
+    GLOB RichardsFlow_2d_small_masslumping_pcs_0_ts_*.vtu saturation saturation 1e-11 1e-15
+    GLOB RichardsFlow_2d_small_masslumping_pcs_0_ts_*.vtu saturation_avg saturation_avg 1e-11 1e-15
+    GLOB RichardsFlow_2d_small_masslumping_pcs_0_ts_*.vtu velocity velocity 1e-15 1e-15
+    GLOB RichardsFlow_2d_small_masslumping_pcs_0_ts_*.vtu HydraulicFlow HydraulicFlow 1e-13 0
+    GLOB RichardsFlow_2d_small_masslumping_pcs_0_ts_*.vtu NodalForces NodalForces 1e-9 0
+)
 AddTest(
     NAME RichardsMechanics_RichardsFlow_2d_quasinewton
     PATH RichardsMechanics

Tests/Data/RichardsMechanics/RichardsFlow_2d_small_masslumping.prj

+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <mesh>Richards_2d.vtu</mesh>
+    <geometry>Richards_2d.gml</geometry>
+    <processes>
+        <process>
+            <name>RM</name>
+            <type>RICHARDS_MECHANICS</type>
+            <integration_order>3</integration_order>
+            <dimension>2</dimension>
+            <constitutive_relation>
+                <type>LinearElasticIsotropic</type>
+                <youngs_modulus>E</youngs_modulus>
+                <poissons_ratio>nu</poissons_ratio>
+            </constitutive_relation>
+            <solid_bulk_modulus>K_SR</solid_bulk_modulus>
+            <fluid_bulk_modulus>K_LR</fluid_bulk_modulus>
+            <biot_coefficient>alpha</biot_coefficient>
+            <solid_density>rho_sr</solid_density>
+            <process_variables>
+                <displacement>displacement</displacement>
+                <pressure>pressure</pressure>
+            </process_variables>
+            <material_property>
+                <fluid>
+                    <density>
+                        <type>Constant</type>
+                        <value> 1 </value>
+                    </density>
+                    <viscosity>
+                        <type>Constant</type>
+                        <value> 1.e-3 </value>
+                    </viscosity>
+                </fluid>
+                <porous_medium>
+                    <porous_medium id="0">
+                        <permeability>
+                            <permeability_tensor_entries>k</permeability_tensor_entries>
+                            <type>Constant</type>
+                        </permeability>
+                        <porosity>
+                            <type>Constant</type>
+                            <porosity_parameter>phi</porosity_parameter>
+                        </porosity>
+                        <storage>
+                            <type>Constant</type>
+                            <value> 0.0 </value>
+                        </storage>
+                        <capillary_pressure>
+                            <type>vanGenuchten</type>
+                            <m>0.789029535864979</m>
+                            <pc_max>36333.30</pc_max>
+                            <pd>3633.33</pd>
+                            <smax>0.95</smax>
+                            <sr>0.1689</sr>
+                        </capillary_pressure>
+                        <relative_permeability>
+                            <type>WettingPhaseVanGenuchten</type>
+                            <krel_min>1e-12</krel_min>
+                            <m>0.789029535864979</m>
+                            <smax>0.95</smax>
+                            <sr>0.1689</sr>
+                        </relative_permeability>
+                    </porous_medium>
+                </porous_medium>
+            </material_property>
+            <secondary_variables>
+                <secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
+                <secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
+                <secondary_variable type="static" internal_name="velocity" output_name="velocity"/>
+                <secondary_variable type="static" internal_name="saturation" output_name="saturation"/>
+            </secondary_variables>
+            <specific_body_force>0 -9.81</specific_body_force>
+            <mass_lumping>true</mass_lumping>
+            <temperature> temp </temperature>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="RM">
+                <nonlinear_solver>nonlinear_solver</nonlinear_solver>
+                <convergence_criterion>
+                    <type>PerComponentDeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstols>1e-7 1e-11 1e-10</abstols>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial>0</t_initial>
+                    <t_end>2000</t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>1</repeat>
+                            <delta_t>1</delta_t>
+                        </pair>
+                        <pair>
+                            <repeat>1</repeat>
+                            <delta_t>2</delta_t>
+                        </pair>
+                        <pair>
+                            <repeat>1</repeat>
+                            <delta_t>5</delta_t>
+                        </pair>
+                        <pair>
+                            <repeat>1</repeat>
+                            <delta_t>10</delta_t>
+                        </pair>
+                        <pair>
+                            <repeat>50</repeat>
+                            <delta_t>20</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>RichardsFlow_2d_small_masslumping</prefix>
+            <timesteps>
+                <pair>
+                    <repeat>9</repeat>
+                    <each_steps>1</each_steps>
+                </pair>
+                <pair>
+                    <repeat>100</repeat>
+                    <each_steps>10</each_steps>
+                </pair>
+            </timesteps>
+            <output_iteration_results>false</output_iteration_results>
+            <variables>
+                <variable>displacement</variable>
+                <variable>pressure</variable>
+                <variable>sigma</variable>
+                <variable>epsilon</variable>
+                <variable>velocity</variable>
+                <variable>saturation</variable>
+            </variables>
+        </output>
+    </time_loop>
+    <parameters>
+        <!-- Mechanics -->
+        <parameter>
+            <name>K_SR</name>
+            <type>Constant</type>
+            <value>160e300</value>
+        </parameter>
+        <parameter>
+            <name>E</name>
+            <type>Constant</type>
+            <value>5e9</value>
+        </parameter>
+        <parameter>
+            <name>nu</name>
+            <type>Constant</type>
+            <value>.1</value>
+        </parameter>
+        <!-- Model parameters -->
+        <parameter>
+            <name>k</name>
+            <type>Constant</type>
+            <value>4.46e-13</value>
+        </parameter>
+        <parameter>
+            <name>alpha</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>phi</name>
+            <type>Constant</type>
+            <value>0.38</value>
+        </parameter>
+        <parameter>
+            <name>rho_sr</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>K_LR</name>
+            <type>Constant</type>
+            <value>2.2e300</value>
+        </parameter>
+        <parameter>
+            <name>temp</name>
+            <type>Constant</type>
+            <value>293.15</value>
+        </parameter>
+        <parameter>
+            <name>displacement0</name>
+            <type>Constant</type>
+            <values>0 0</values>
+        </parameter>
+        <parameter>
+            <name>dirichlet0</name>
+            <type>Constant</type>
+            <values>0</values>
+        </parameter>
+        <parameter>
+            <name>p0</name>
+            <type>Constant</type>
+            <value>-5000</value>
+        </parameter>
+        <parameter>
+            <name>p_Dirichlet_bottom</name>
+            <type>Constant</type>
+            <value>1000</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>displacement</name>
+            <components>2</components>
+            <order>2</order>
+            <initial_condition>displacement0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>Richards_2d_geometry</geometrical_set>
+                    <geometry>LEFT</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>Richards_2d_geometry</geometrical_set>
+                    <geometry>RIGHT</geometry>
+                    <type>Dirichlet</type>
+                    <component>0</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>Richards_2d_geometry</geometrical_set>
+                    <geometry>BOTTOM</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>Richards_2d_geometry</geometrical_set>
+                    <geometry>TOP</geometry>
+                    <type>Dirichlet</type>
+                    <component>1</component>
+                    <parameter>dirichlet0</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+        <process_variable>
+            <name>pressure</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>p0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <geometrical_set>Richards_2d_geometry</geometrical_set>
+                    <geometry>BOTTOM</geometry>
+                    <type>Dirichlet</type>
+                    <parameter>p_Dirichlet_bottom</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <geometrical_set>Richards_2d_geometry</geometrical_set>
+                    <geometry>TOP</geometry>
+                    <type>Dirichlet</type>
+                    <parameter>p0</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>nonlinear_solver</name>
+            <type>Newton</type>
+            <max_iter>100</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+            <!--damping>0.9</damping-->
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <eigen>
+                <solver_type>SparseLU</solver_type>
+                <scaling>true</scaling>
+                <!--
+                <solver_type>PardisoLU</solver_type>
+                -->
+            </eigen>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>