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.

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

RichardsMechanics implementation. #2189

Merged
merged 6 commits into from Aug 22, 2018

Conversation

Projects
None yet
2 participants
@endJunction
Copy link
Member

endJunction commented Aug 19, 2018

Based on #2188

Monolithic scheme. (Staggered scheme functions are prepared, but not implemented.)
Current implementation relies on the RichardsFlow material and solid material model implementations. Some of the properties are repeated therefore, and not all properties (of RichardsFlow material) are used, storage term, for example...
Tests for coupled and decoupled saturated flow are provided. Decoupled unsaturated test (with comparison to the RichardsFlow implementation) is given. There is no test for coupled unsaturated flow for now.

#ifdef OGS_BUILD_PROCESS_RICHARDSMECHANICS
if (type == "RICHARDS_MECHANICS")
{
//! \ogs_file_param{prj__processes__process__HYDRO_MECHANICS__dimension}

This comment has been minimized.

@wenqing

wenqing Aug 20, 2018

Member

...__RICHARDS_MECHANIC_dimension

{
OGS_FATAL("Not implemented.");
/*
double const Sg = 1 - saturation;

This comment has been minimized.

@wenqing

wenqing Aug 20, 2018

Member

Will these lines be used later on?


S_L = _process_data.flow_material->getSaturation(
material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
if (S_L < 0 || S_L > 1)

This comment has been minimized.

@wenqing

wenqing Aug 20, 2018

Member

This if-condition can be removed because the calculated saturation has already been restricted in the range of [0,1] in the member function of capillary models call by _process_data.flow_material->getSaturation.

double const specific_storage =
dS_L_dp_cap * (p_cap_ip * a0 - porosity) +
S_L * (porosity / K_LR + a0);
M.template block<pressure_size, pressure_size>(pressure_index,

This comment has been minimized.

@wenqing

wenqing Aug 20, 2018

Member

Alternatively, the PDE can be scaled with rho_LR, which gives more accurate solution for the problems with known fluid velocity on the boundary and with varying fluid density.

This comment has been minimized.

@wenqing

wenqing Aug 20, 2018

Member

It is also good without scaling. Then the Neumman BC means the mass per unit area across the boundary.


S_L = _process_data.flow_material->getSaturation(
material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
if (S_L < 0 || S_L > 1)

This comment has been minimized.

@wenqing

wenqing Aug 20, 2018

Member

This if-condition can be removed.

NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
double const S_L = _process_data.flow_material->getSaturation(
material_id, t, x_position, -p_cap_ip, temperature, p_cap_ip);
if (S_L < 0 || S_L > 1)

This comment has been minimized.

@wenqing

wenqing Aug 20, 2018

Member

This can be removed.

.template block<pressure_size, pressure_size>(pressure_index,
pressure_index)
.noalias() -=
N_p.transpose() * rho_LR * dspecific_storage_dp_cap * N_p * w;

This comment has been minimized.

@wenqing

wenqing Aug 20, 2018

Member

as discussed, dot p_c is missing.

dS_L_dp_cap *
(porosity / K_LR + a0 * 3 * S_L + dS_L_dp_cap * p_cap_ip * a0);

storage_p.noalias() +=

This comment has been minimized.

@wenqing

wenqing Aug 20, 2018

Member

I see that the pressure part is not scaled with rho_LR. It is OK.

@endJunction endJunction force-pushed the endJunction:RM branch 4 times, most recently from 7cb2bfa to 7d6a63e Aug 20, 2018

@endJunction endJunction removed the WIP 🏗 label Aug 21, 2018

@endJunction

This comment has been minimized.

Copy link
Member Author

endJunction commented Aug 21, 2018

TODO's after merge:

  • Use GPML (need some time with NG)
  • Swelling is not included.
  • Test of coupled unsaturated flow (Vin?)
  • Add docu pdf (@nagelt ?)
@wenqing
Copy link
Member

wenqing left a comment

Looks good except one sign.

pressure_index)
.noalias() -= N_p.transpose() * rho_LR * p_cap_dot_ip *
dspecific_storage_dp_cap * N_p * w;

This comment has been minimized.

@wenqing

wenqing Aug 21, 2018

Member

-= should be +=, because
d (alpha dot p)/dp = d (alpha)/dp dot p + alpha d(dot p)/dp. In which
d (alpha)/dp dot p = d (alpha)/dp_c dp_c/dp dot p = d (alpha)/dp_c (-1) (-dot p_c) = d (alpha)/dp_c dot p_c)

This comment has been minimized.

@endJunction

endJunction Aug 21, 2018

Author Member

What you say, is correct, though the convergence is worse. I'm looking into this, but if there is no improvement, I'll omit the term completely for now...

@endJunction endJunction force-pushed the endJunction:RM branch from 39f7477 to 5661b33 Aug 22, 2018

@endJunction endJunction force-pushed the endJunction:RM branch from 5661b33 to 86ff420 Aug 22, 2018

@endJunction endJunction force-pushed the endJunction:RM branch from 86ff420 to 73e6aa9 Aug 22, 2018

@endJunction endJunction force-pushed the endJunction:RM branch from 73e6aa9 to 22264ac Aug 22, 2018

@endJunction

This comment has been minimized.

Copy link
Member Author

endJunction commented Aug 22, 2018

Jacobian fixed and corresponds to the numerical Jacobian in RM and RF.
RF example with 2nd order integration delivers slightly different results and was changed to use the 3rd order for comparison with RM. Since density is cancelled out in RF, for comparison it is set to 1.

@endJunction endJunction merged commit ff4ed73 into ufz:master Aug 22, 2018

3 checks passed

continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/jenkins/pr-merge This commit looks good
Details
deploy/netlify Deploy preview ready!
Details

@endJunction endJunction deleted the endJunction:RM branch Aug 22, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.