OGS
ComputeMicroPorosity.h
Go to the documentation of this file.
1
10#pragma once
11
12#include <Eigen/LU>
13#include <cassert>
14#include <iosfwd>
15
20
22{
23template <int DisplacementDim>
42
43template <int DisplacementDim>
44std::ostream& operator<<(std::ostream& os,
46{
47 return os << "phi_m: " << state.phi_m << ", "
48 << "e_sw: " << state.e_sw << ", "
49 << "p_L_m: " << state.p_L_m << ", "
50 << "sigma_sw: " << state.sigma_sw.transpose();
51}
52
62
63template <int DisplacementDim>
66 I_2_C_el_inverse,
67 double const rho_LR_m, // for simplification equal to rho_LR_M
68 double const mu_LR,
69 MicroPorosityParameters const& micro_porosity_parameters,
70 double const alpha_B, double const phi, double const p_L,
71 double const p_L_m_prev,
72 MaterialPropertyLib::VariableArray const& /*variables_prev*/,
73 double const S_L_m_prev, double const phi_m_prev,
74 ParameterLib::SpatialPosition const pos, double const t, double const dt,
75 MaterialPropertyLib::Property const& saturation_micro,
76 MaterialPropertyLib::Property const& swelling_stress_rate)
77{
78 namespace MPL = MaterialPropertyLib;
79 static constexpr int kelvin_vector_size =
81 // phi_m, e_sw, p_L_m, sigma_sw
82 static constexpr int nls_size = 1 + 1 + 1 + kelvin_vector_size;
83
84 static constexpr int i_phi_m = 0;
85 static constexpr int i_e_sw = 1;
86 static constexpr int i_p_L_m = 2;
87 static constexpr int i_sigma_sw = 3;
88
89 using ResidualVectorType = Eigen::Matrix<double, nls_size, 1>;
90 using JacobianMatrix =
91 Eigen::Matrix<double, nls_size, nls_size, Eigen::RowMajor>;
92
93 Eigen::FullPivLU<Eigen::Matrix<double, nls_size, nls_size, Eigen::RowMajor>>
94 linear_solver;
95
96 JacobianMatrix jacobian;
97
98 // Agglomerated solution vector construction.
99 ResidualVectorType solution = ResidualVectorType::Zero();
100
101 if (p_L >= 0 && p_L_m_prev >= 0)
102 {
103 return {solution[i_phi_m], solution[i_e_sw], solution[i_p_L_m],
104 solution.template segment<kelvin_vector_size>(i_sigma_sw)};
105 }
106
107 double const alpha_bar =
108 micro_porosity_parameters.mass_exchange_coefficient;
109
110 auto const update_residual = [&](ResidualVectorType& residual)
111 {
112 double const delta_phi_m = solution[i_phi_m];
113 double const delta_e_sw = solution[i_e_sw];
114 auto const& delta_sigma_sw =
115 solution.template segment<kelvin_vector_size>(i_sigma_sw);
116 double const delta_p_L_m = solution[i_p_L_m];
117
118 double const phi_m = phi_m_prev + delta_phi_m;
119 double const p_L_m = p_L_m_prev + delta_p_L_m;
120 MPL::VariableArray variables_prev;
121 variables_prev.capillary_pressure = -p_L_m_prev;
122 variables_prev.liquid_saturation = S_L_m_prev;
123
124 MPL::VariableArray variables;
125 variables.capillary_pressure = -p_L_m;
126
127 double const S_L_m =
128 saturation_micro.template value<double>(variables, pos, t, dt);
129 variables.liquid_saturation = S_L_m;
130 double const delta_S_L_m = S_L_m - S_L_m_prev;
131
132 auto const sigma_sw_dot =
133 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
134 swelling_stress_rate.template value<Eigen::Matrix3d>(
135 variables, variables_prev, pos, t, dt));
136
137 residual[i_phi_m] = delta_phi_m - (alpha_B - phi) * delta_e_sw;
138 residual[i_e_sw] = delta_e_sw + I_2_C_el_inverse.dot(delta_sigma_sw);
139 residual.template segment<kelvin_vector_size>(i_sigma_sw).noalias() =
140 delta_sigma_sw - sigma_sw_dot * dt;
141
142 residual[i_p_L_m] =
143 rho_LR_m *
144 (phi_m * delta_S_L_m - (alpha_B - phi) * S_L_m * delta_e_sw) +
145 phi_m * S_L_m * rho_LR_m * delta_e_sw -
146 micro_porosity_parameters.mass_exchange_coefficient / mu_LR *
147 (p_L - p_L_m) * dt;
148 };
149
150 auto const update_jacobian = [&](JacobianMatrix& jacobian)
151 {
152 jacobian = JacobianMatrix::Identity();
153
154 double const delta_phi_m = solution[i_phi_m];
155 double const delta_e_sw = solution[i_e_sw];
156 double const delta_p_L_m = solution[i_p_L_m];
157
158 double const phi_m = phi_m_prev + delta_phi_m;
159 double const p_L_m = p_L_m_prev + delta_p_L_m;
160 MPL::VariableArray variables_prev;
161 variables_prev.capillary_pressure = -p_L_m_prev;
162 MPL::VariableArray variables;
163 variables.capillary_pressure = -p_L_m;
164
165 double const S_L_m =
166 saturation_micro.template value<double>(variables, pos, t, dt);
167 variables_prev.liquid_saturation = S_L_m_prev;
168 variables.liquid_saturation = S_L_m;
169 double const delta_S_L_m = S_L_m - S_L_m_prev;
170
171 double const dS_L_m_dp_cap_m = saturation_micro.template dValue<double>(
172 variables, MPL::Variable::capillary_pressure, pos, t, dt);
173 auto const dsigma_sw_dS_L_m =
174 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
175 swelling_stress_rate.template dValue<Eigen::Matrix3d>(
176 variables, variables_prev, MPL::Variable::liquid_saturation,
177 pos, t, dt));
178
179 jacobian(i_phi_m, i_e_sw) = -(alpha_B - phi);
180
181 jacobian.template block<1, kelvin_vector_size>(i_e_sw, i_sigma_sw) =
182 I_2_C_el_inverse.transpose();
183
184 jacobian.template block<kelvin_vector_size, 1>(i_sigma_sw, i_p_L_m) =
185 -dsigma_sw_dS_L_m * dS_L_m_dp_cap_m;
186
187 jacobian(i_p_L_m, i_phi_m) =
188 rho_LR_m * (delta_S_L_m + S_L_m * delta_e_sw);
189
190 jacobian(i_p_L_m, i_e_sw) = -rho_LR_m * S_L_m * (alpha_B - phi - phi_m);
191
192 jacobian(i_p_L_m, i_p_L_m) =
193 alpha_bar / mu_LR * dt -
194 rho_LR_m * (phi_m - (alpha_B - phi - phi_m) * delta_e_sw) *
195 dS_L_m_dp_cap_m;
196 };
197
198 auto const update_solution = [&](ResidualVectorType const& increment)
199 { solution += increment; };
200
201 auto newton_solver =
202 NumLib::NewtonRaphson<decltype(linear_solver), JacobianMatrix,
203 decltype(update_jacobian), ResidualVectorType,
204 decltype(update_residual),
205 decltype(update_solution)>(
206 linear_solver, update_jacobian, update_residual, update_solution,
207 micro_porosity_parameters.nonlinear_solver_parameters);
208
209 auto const success_iterations = newton_solver.solve(jacobian);
210
211 if (!success_iterations)
212 {
213 OGS_FATAL(
214 "Could not find solution for local double structure nonlinear "
215 "problem.");
216 }
217
218 return {solution[i_phi_m], solution[i_e_sw], solution[i_p_L_m],
219 solution.template segment<kelvin_vector_size>(i_sigma_sw)};
220}
221} // namespace ProcessLib::RichardsMechanics
#define OGS_FATAL(...)
Definition Error.h:26
std::optional< int > solve(JacobianMatrix &jacobian) const
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
std::ostream & operator<<(std::ostream &os, MicroPorosityStateSpace< DisplacementDim > const &state)
MicroPorosityStateSpace< DisplacementDim > computeMicroPorosity(MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &I_2_C_el_inverse, double const rho_LR_m, double const mu_LR, MicroPorosityParameters const &micro_porosity_parameters, double const alpha_B, double const phi, double const p_L, double const p_L_m_prev, MaterialPropertyLib::VariableArray const &, 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 nonlinear_solver_parameters
MicroPorosityStateSpace & operator+=(MicroPorosityStateSpace const &state)
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > sigma_sw