Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

RichardsMechanics: mass lumping #2590

Merged
merged 3 commits into from
Aug 6, 2019
Merged

RichardsMechanics: mass lumping #2590

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
944d6b8
[PL] RM; Parse mass_lumping tag
endJunction Aug 5, 2019
128af5d
[T] RM; Newton with mass lumping.
endJunction Aug 5, 2019
dc69d0e
[PL] RM; Implement mass lumping.
endJunction Aug 5, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions ...entation/ProjectFile/prj/processes/process/RICHARDS_MECHANICS/t_mass_lumping.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
A tag for enabling diagonal lumping of the mass matrix.
7 changes: 6 additions & 1 deletion ProcessLib/RichardsMechanics/CreateRichardsMechanicsProcess.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand All @@ -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;

Expand Down
12 changes: 12 additions & 0 deletions ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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();
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In future, the name of this term is better to be unified both here and where in assemble().

}

// pressure equation, pressure part.
local_Jac
.template block<pressure_size, pressure_size>(pressure_index,
Expand Down
8 changes: 6 additions & 2 deletions ProcessLib/RichardsMechanics/RichardsMechanicsProcessData.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_)},
Expand All @@ -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_)
{
}

Expand Down Expand Up @@ -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;
};

Expand Down
20 changes: 20 additions & 0 deletions ProcessLib/RichardsMechanics/Tests.cmake
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
291 changes: 291 additions & 0 deletions Tests/Data/RichardsMechanics/RichardsFlow_2d_small_masslumping.prj
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,291 @@
<?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>
Loading