Commit 30f93664 authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov

Merge branch 'RM_UseDoublePorosityModel' into 'master'

RM use double porosity model

See merge request ogs/ogs!3441
parents d9825ecb 10e698bf
An optional tag enables double structure porosity model. It requires
micro-saturation property to be available in the material properties.
\copydoc ProcessLib::RichardsMechanics::MicroPorosityParameters::mass_exchange_coefficient
......@@ -49,20 +49,30 @@ std::ostream& operator<<(std::ostream& os,
<< "sigma_sw: " << state.sigma_sw.transpose();
}
struct MicroPorosityParameters
{
NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters;
/// Coefficient \f$\bar\alpha\f$ of the micro structure mass exchange. If
/// this is given then the micro_saturation MPL property must be provided
/// too.
double mass_exchange_coefficient;
};
template <int DisplacementDim>
MicroPorosityStateSpace<DisplacementDim> computeMicroPorosity(
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const&
I_2_C_el_inverse,
double const rho_LR_m, // for simplification equal to rho_LR_M
double const mu_LR, double const alpha_bar, double const alpha_B,
double const phi_M, double const p_L, double const p_L_m_prev,
double const mu_LR,
MicroPorosityParameters const& micro_porosity_parameters,
double const alpha_B, double const phi_M, double const p_L,
double const p_L_m_prev,
MaterialPropertyLib::VariableArray const& /*variables_prev*/,
double const S_L_m_prev, double const phi_m_prev,
ParameterLib::SpatialPosition const pos, double const t, double const dt,
MaterialPropertyLib::Property const& saturation_micro,
MaterialPropertyLib::Property const& swelling_stress_rate,
NumLib::NewtonRaphsonSolverParameters const& nonlinear_solver_parameters)
MaterialPropertyLib::Property const& swelling_stress_rate)
{
namespace MPL = MaterialPropertyLib;
static constexpr int kelvin_vector_size =
......@@ -93,6 +103,9 @@ MicroPorosityStateSpace<DisplacementDim> computeMicroPorosity(
solution.template segment<kelvin_vector_size>(i_sigma_sw)};
}
double const alpha_bar =
micro_porosity_parameters.mass_exchange_coefficient;
auto const update_residual = [&](ResidualVectorType& residual) {
double const delta_phi_m = solution[i_phi_m];
double const delta_e_sw = solution[i_e_sw];
......@@ -131,7 +144,8 @@ MicroPorosityStateSpace<DisplacementDim> computeMicroPorosity(
rho_LR_m *
(phi_m * delta_S_L_m - (alpha_B - phi_M) * S_L_m * delta_e_sw) +
phi_m * S_L_m * rho_LR_m * delta_e_sw -
alpha_bar / mu_LR * (p_L - p_L_m) * dt;
micro_porosity_parameters.mass_exchange_coefficient / mu_LR *
(p_L - p_L_m) * dt;
};
auto const update_jacobian = [&](JacobianMatrix& jacobian) {
......@@ -172,9 +186,8 @@ MicroPorosityStateSpace<DisplacementDim> computeMicroPorosity(
jacobian.template block<kelvin_vector_size, 1>(
i_sigma_sw, i_p_L_m) = -dsigma_sw_dS_L_m * dS_L_m_dp_cap_m;
jacobian(i_p_L_m, i_phi_m) =
rho_LR_m * (delta_S_L_m + phi_m * delta_e_sw);
rho_LR_m * (delta_S_L_m + S_L_m * delta_e_sw);
jacobian(i_p_L_m, i_e_sw) = -rho_LR_m * S_L_m * (alpha_B - phi_M - phi_m);
......@@ -194,7 +207,7 @@ MicroPorosityStateSpace<DisplacementDim> computeMicroPorosity(
decltype(update_residual),
decltype(update_solution)>(
linear_solver, update_jacobian, update_residual, update_solution,
nonlinear_solver_parameters);
micro_porosity_parameters.nonlinear_solver_parameters);
auto const success_iterations = newton_solver.solve(jacobian);
......
......@@ -17,6 +17,7 @@
#include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
#include "MaterialLib/SolidModels/MechanicsBase.h"
#include "NumLib/CreateNewtonRaphsonSolverParameters.h"
#include "ParameterLib/Utils.h"
#include "ProcessLib/Output/CreateSecondaryVariables.h"
#include "ProcessLib/Utils/ProcessUtils.h"
......@@ -169,6 +170,20 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
MathLib::KelvinVector::KelvinVectorDimensions<DisplacementDim>::value,
&mesh);
std::optional<MicroPorosityParameters> micro_porosity_parameters;
if (auto const micro_porosity_config =
//! \ogs_file_param{prj__processes__process__RICHARDS_MECHANICS__micro_porosity}
config.getConfigSubtreeOptional("micro_porosity"))
{
micro_porosity_parameters = MicroPorosityParameters{
NumLib::createNewtonRaphsonSolverParameters(
//! \ogs_file_param{prj__processes__process__RICHARDS_MECHANICS__micro_porosity__nonlinear_solver}
micro_porosity_config->getConfigSubtree("nonlinear_solver")),
//! \ogs_file_param{prj__processes__process__RICHARDS_MECHANICS__micro_porosity__mass_exchange_coefficient}
micro_porosity_config->getConfigParameter<double>(
"mass_exchange_coefficient")};
}
auto const mass_lumping =
//! \ogs_file_param{prj__processes__process__RICHARDS_MECHANICS__mass_lumping}
config.getConfigParameter<bool>("mass_lumping", false);
......@@ -183,6 +198,7 @@ std::unique_ptr<Process> createRichardsMechanicsProcess(
std::move(solid_constitutive_relations),
initial_stress,
specific_body_force,
micro_porosity_parameters,
mass_lumping,
explicit_hm_coupling_in_unsaturated_zone};
......
......@@ -60,8 +60,12 @@ struct IntegrationPointData final
typename ShapeMatricesTypePressure::GlobalDimVectorType v_darcy;
double liquid_pressure_m = std::numeric_limits<double>::quiet_NaN();
double liquid_pressure_m_prev = std::numeric_limits<double>::quiet_NaN();
double saturation = std::numeric_limits<double>::quiet_NaN();
double saturation_prev = std::numeric_limits<double>::quiet_NaN();
double saturation_m = std::numeric_limits<double>::quiet_NaN();
double saturation_m_prev = std::numeric_limits<double>::quiet_NaN();
double porosity = std::numeric_limits<double>::quiet_NaN();
double porosity_prev = std::numeric_limits<double>::quiet_NaN();
double transport_porosity = std::numeric_limits<double>::quiet_NaN();
......@@ -85,8 +89,10 @@ struct IntegrationPointData final
sigma_eff_prev = sigma_eff;
sigma_sw_prev = sigma_sw;
saturation_prev = saturation;
saturation_m_prev = saturation_m;
porosity_prev = porosity;
transport_porosity_prev = transport_porosity;
liquid_pressure_m_prev = liquid_pressure_m;
material_state_variables->pushBackState();
}
......
......@@ -64,6 +64,22 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const = 0;
virtual std::vector<double> getMicroSaturation() const = 0;
virtual std::vector<double> const& getIntPtMicroSaturation(
const double t,
std::vector<GlobalVector*> const& x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const = 0;
virtual std::vector<double> getMicroPressure() const = 0;
virtual std::vector<double> const& getIntPtMicroPressure(
const double t,
std::vector<GlobalVector*> const& x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const = 0;
virtual std::vector<double> getPorosity() const = 0;
virtual std::vector<double> const& getIntPtPorosity(
......
......@@ -13,6 +13,7 @@
#include <cassert>
#include "ComputeMicroPorosity.h"
#include "IntegrationPointData.h"
#include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
......@@ -219,6 +220,8 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
p_cap_ip;
variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
_ip_data[ip].liquid_pressure_m_prev = -p_cap_ip;
_ip_data[ip].liquid_pressure_m = -p_cap_ip;
auto const temperature =
medium->property(MPL::PropertyType::reference_temperature)
......@@ -229,6 +232,17 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
medium->property(MPL::PropertyType::saturation)
.template value<double>(variables, x_position, t, dt);
if (medium->hasProperty(MPL::PropertyType::saturation_micro))
{
MPL::VariableArray vars;
vars[static_cast<int>(MPL::Variable::capillary_pressure)] =
p_cap_ip;
double const S_L_m =
medium->property(MPL::PropertyType::saturation_micro)
.template value<double>(vars, x_position, t, dt);
_ip_data[ip].saturation_m_prev = S_L_m;
}
// Set eps_m_prev from potentially non-zero eps and sigma_sw from
// restart.
auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
......@@ -801,8 +815,9 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
// Swelling and possibly volumetric strain rate update.
auto& sigma_sw = _ip_data[ip].sigma_sw;
auto const& sigma_sw_prev = _ip_data[ip].sigma_sw_prev;
if (!medium->hasProperty(MPL::PropertyType::saturation_micro))
{
auto const& sigma_sw_prev = _ip_data[ip].sigma_sw_prev;
// If there is swelling, compute it. Update volumetric strain rate,
// s.t. it corresponds to the mechanical part only.
......@@ -828,33 +843,79 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
MPL::Variable::volumetric_strain)]) +=
identity2.transpose() * C_el.inverse() * sigma_sw_prev;
}
}
if (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
{
variables_prev[static_cast<int>(
MPL::Variable::transport_porosity)] =
_ip_data[ip].transport_porosity_prev;
auto const mu =
liquid_phase.property(MPL::PropertyType::viscosity)
.template value<double>(variables, x_position, t, dt);
_ip_data[ip].transport_porosity =
solid_phase.property(MPL::PropertyType::transport_porosity)
.template value<double>(variables, variables_prev,
x_position, t, dt);
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
_ip_data[ip].transport_porosity;
}
else
{
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
phi;
// TODO (naumov) saturation_micro must be always defined together with
// the micro_porosity_parameters.
if (medium->hasProperty(MPL::PropertyType::saturation_micro))
{
double const phi_M_prev = _ip_data[ip].transport_porosity_prev;
double const phi_prev = _ip_data[ip].porosity_prev;
double const phi_m_prev = phi_prev - phi_M_prev;
double const p_L_m_prev = _ip_data[ip].liquid_pressure_m_prev;
auto const S_L_m_prev = _ip_data[ip].saturation_m_prev;
auto const [delta_phi_m, delta_e_sw, delta_p_L_m, delta_sigma_sw] =
computeMicroPorosity<DisplacementDim>(
identity2.transpose() * C_el.inverse(), rho_LR, mu,
*_process_data.micro_porosity_parameters, alpha, phi,
-p_cap_ip, p_L_m_prev, variables_prev, S_L_m_prev,
phi_m_prev, x_position, t, dt,
medium->property(MPL::PropertyType::saturation_micro),
solid_phase.property(
MPL::PropertyType::swelling_stress_rate));
auto& phi_M = _ip_data[ip].transport_porosity;
phi_M = phi - (phi_m_prev + delta_phi_m);
variables_prev[static_cast<int>(
MPL::Variable::transport_porosity)] = phi_M_prev;
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
phi_M;
auto& p_L_m = _ip_data[ip].liquid_pressure_m;
p_L_m = p_L_m_prev + delta_p_L_m;
{ // Update micro saturation.
MPL::VariableArray variables_prev;
variables_prev[static_cast<int>(
MPL::Variable::capillary_pressure)] = -p_L_m_prev;
MPL::VariableArray variables;
variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
-p_L_m;
_ip_data[ip].saturation_m =
medium->property(MPL::PropertyType::saturation_micro)
.template value<double>(variables, x_position, t, dt);
}
sigma_sw = sigma_sw_prev + delta_sigma_sw;
}
if (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
{
variables_prev[static_cast<int>(
MPL::Variable::transport_porosity)] =
_ip_data[ip].transport_porosity_prev;
_ip_data[ip].transport_porosity =
solid_phase.property(MPL::PropertyType::transport_porosity)
.template value<double>(variables, variables_prev,
x_position, t, dt);
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
_ip_data[ip].transport_porosity;
}
else
{
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
phi;
}
double const k_rel =
medium->property(MPL::PropertyType::relative_permeability)
.template value<double>(variables, x_position, t, dt);
auto const mu =
liquid_phase.property(MPL::PropertyType::viscosity)
.template value<double>(variables, x_position, t, dt);
// Set mechanical variables for the intrinsic permeability model
// For stress dependent permeability.
......@@ -934,7 +995,17 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
.noalias() +=
N_u_op.transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
// For the swelling stress with double structure model the corresponding
// Jacobian u-p entry would be required, but it does not improve
// convergence and sometimes worsens it:
// if (medium->hasProperty(MPL::PropertyType::saturation_micro))
// {
// -B.transpose() *
// dsigma_sw_dS_L_m* dS_L_m_dp_cap_m*(p_L_m - p_L_m_prev) /
// p_cap_dot_ip / dt* N_p* w;
// }
if (!medium->hasProperty(MPL::PropertyType::saturation_micro) &&
solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
{
using DimMatrix = Eigen::Matrix<double, 3, 3>;
auto const dsigma_sw_dS_L =
......@@ -1036,6 +1107,31 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
local_rhs.template segment<pressure_size>(pressure_index).noalias() +=
dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
if (medium->hasProperty(MPL::PropertyType::saturation_micro))
{
double const alpha_bar = _process_data.micro_porosity_parameters
->mass_exchange_coefficient;
auto const p_L_m = _ip_data[ip].liquid_pressure_m;
local_rhs.template segment<pressure_size>(pressure_index)
.noalias() -=
N_p.transpose() * alpha_bar / mu * (-p_cap_ip - p_L_m) * w;
local_Jac
.template block<pressure_size, pressure_size>(pressure_index,
pressure_index)
.noalias() += N_p.transpose() * alpha_bar / mu * N_p * w;
if (p_cap_dot_ip != 0)
{
double const p_L_m_prev = _ip_data[ip].liquid_pressure_m_prev;
local_Jac
.template block<pressure_size, pressure_size>(
pressure_index, pressure_index)
.noalias() += N_p.transpose() * alpha_bar / mu *
(p_L_m - p_L_m_prev) / (dt * p_cap_dot_ip) *
N_p * w;
}
}
}
if (_process_data.apply_mass_lumping)
......@@ -1246,6 +1342,58 @@ std::vector<double> const& RichardsMechanicsLocalAssembler<
_ip_data, &IpData::saturation, cache);
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
std::vector<double> RichardsMechanicsLocalAssembler<
ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
DisplacementDim>::getMicroSaturation() const
{
std::vector<double> result;
getIntPtMicroSaturation(0, {}, {}, result);
return result;
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
std::vector<double> const& RichardsMechanicsLocalAssembler<
ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
DisplacementDim>::
getIntPtMicroSaturation(
const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& cache) const
{
return ProcessLib::getIntegrationPointScalarData(
_ip_data, &IpData::saturation_m, cache);
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
std::vector<double> RichardsMechanicsLocalAssembler<
ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
DisplacementDim>::getMicroPressure() const
{
std::vector<double> result;
getIntPtMicroPressure(0, {}, {}, result);
return result;
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
std::vector<double> const& RichardsMechanicsLocalAssembler<
ShapeFunctionDisplacement, ShapeFunctionPressure, IntegrationMethod,
DisplacementDim>::
getIntPtMicroPressure(
const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& cache) const
{
return ProcessLib::getIntegrationPointScalarData(
_ip_data, &IpData::liquid_pressure_m, cache);
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
std::vector<double> RichardsMechanicsLocalAssembler<
......
......@@ -189,6 +189,20 @@ public:
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const override;
std::vector<double> getMicroSaturation() const override;
std::vector<double> const& getIntPtMicroSaturation(
const double t,
std::vector<GlobalVector*> const& x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const override;
std::vector<double> getMicroPressure() const override;
std::vector<double> const& getIntPtMicroPressure(
const double t,
std::vector<GlobalVector*> const& x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const override;
std::vector<double> getPorosity() const override;
std::vector<double> const& getIntPtPorosity(
const double t,
......
......@@ -314,6 +314,12 @@ void RichardsMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
add_secondary_variable("saturation", 1,
&LocalAssemblerIF::getIntPtSaturation);
add_secondary_variable("micro_saturation", 1,
&LocalAssemblerIF::getIntPtMicroSaturation);
add_secondary_variable("micro_pressure", 1,
&LocalAssemblerIF::getIntPtMicroPressure);
add_secondary_variable("porosity", 1,
&LocalAssemblerIF::getIntPtPorosity);
......
......@@ -10,14 +10,13 @@
#pragma once
#include "ParameterLib/Parameter.h"
#include <Eigen/Dense>
#include <memory>
#include <utility>
#include <Eigen/Dense>
#include "ComputeMicroPorosity.h"
#include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
#include "ParameterLib/Parameter.h"
namespace MaterialLib
{
......@@ -55,6 +54,8 @@ struct RichardsMechanicsProcessData
/// A vector of displacement dimension's length.
Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
std::optional<MicroPorosityParameters> micro_porosity_parameters;
bool const apply_mass_lumping;
/// If set, improves convergence of the global Newton method in unsaturated
......
......@@ -9,6 +9,7 @@ if (NOT OGS_USE_MPI)
OgsTest(PROJECTFILE RichardsMechanics/RichardsFlow_2d_small.prj RUNTIME 9)
OgsTest(PROJECTFILE RichardsMechanics/RichardsFlow_2d_small_masslumping.prj RUNTIME 10)
OgsTest(PROJECTFILE RichardsMechanics/RichardsFlow_2d_quasinewton.prj RUNTIME 80)
OgsTest(PROJECTFILE RichardsMechanics/double_porosity_swelling.prj)
OgsTest(PROJECTFILE RichardsMechanics/deformation_dependent_porosity.prj RUNTIME 8)
OgsTest(PROJECTFILE RichardsMechanics/deformation_dependent_porosity_swelling.prj RUNTIME 11)
OgsTest(PROJECTFILE RichardsMechanics/orthotropic_power_law_permeability_xyz.prj RUNTIME 80)
......
<?xml version="1.0" encoding="ISO-8859-1"?>
<OpenGeoSysProject>
<mesh>square_1x1_quad_1e2.vtu</mesh>
<geometry>square_1x1.gml</geometry>
<processes>
<process>
<name>RM</name>
<type>RICHARDS_MECHANICS</type>
<integration_order>2</integration_order>
<dimension>2</dimension>
<micro_porosity>
<mass_exchange_coefficient>1e-12</mass_exchange_coefficient>
<nonlinear_solver>
<maximum_iterations>100</maximum_iterations>
<residuum_tolerance>1e-8</residuum_tolerance>
<increment_tolerance>1e-20</increment_tolerance>
</nonlinear_solver>
</micro_porosity>
<constitutive_relation>
<type>LinearElasticIsotropic</type>
<youngs_modulus>E</youngs_modulus>
<poissons_ratio>nu</poissons_ratio>
</constitutive_relation>
<process_variables>
<pressure>pressure</pressure>
<displacement>displacement</displacement>
</process_variables>
<secondary_variables>
<secondary_variable internal_name="sigma" output_name="sigma"/>
<secondary_variable internal_name="swelling_stress" output_name="swelling_stress"/>
<secondary_variable internal_name="epsilon" output_name="epsilon"/>
<secondary_variable internal_name="velocity" output_name="velocity"/>
<secondary_variable internal_name="saturation" output_name="saturation"/>
<secondary_variable internal_name="porosity" output_name="porosity"/>
<secondary_variable internal_name="transport_porosity" output_name="transport_porosity"/>
<secondary_variable internal_name="dry_density_solid" output_name="dry_density_solid"/>
<secondary_variable internal_name="micro_saturation" output_name="micro_saturation"/>
<secondary_variable internal_name="micro_pressure" output_name="micro_pressure"/>
</secondary_variables>
<specific_body_force>0 0</specific_body_force>
<initial_stress>sigma0</initial_stress>
<mass_lumping>true</mass_lumping>
</process>
</processes>
<media>
<medium>
<phases>
<phase>
<type>AqueousLiquid</type>
<properties>
<property>
<name>viscosity</name>
<type>Constant</type>
<value>1e-3</value>
</property>
<property>
<name>density</name>
<type>Constant</type>
<value>1e3</value>
</property>
</properties>
</phase>
<phase>
<type>Solid</type>
<properties>
<property>
<name>density</name>
<type>Constant</type>
<value>2e3</value>
</property>
<property>
<name>biot_coefficient</name>
<type>Constant</type>
<value>0.6</value>
</property>
<property>
<name>permeability</name>
<type>Constant</type>
<value>1e-21</value>
</property>
<property>
<name>porosity</name>
<type>PorosityFromMassBalance</type>
<initial_porosity>phi0</initial_porosity>
<minimal_porosity>0</minimal_porosity>
<maximal_porosity>1</maximal_porosity>
</property>