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

Create the two phase flow process #1530

Merged
merged 2 commits into from
Nov 20, 2016

Conversation

Yonghui56
Copy link
Contributor

In this pull request, a volume based two phase flow process is created. Nonwetting phase pressure and capillary pressure are selected as the primary variables. The Liakopoulos experiment is selected as a benchmark for testing.
Note that currently the capillary pressure - saturation, relative permeability --saturation relationships are delivered by curves. Later on these will be updated based on the material properties work @wenqing

@@ -371,6 +371,16 @@ if(NOT OGS_USE_MPI)
DIFF_DATA
h_us_quad_1000.vtu richards_pcs_0_ts_100_t_100.000000.vtu PRESSURE1 pressure
)
AddTest(
Copy link
Member

Choose a reason for hiding this comment

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

tabs to be replaced with whitespaces. Please change your editor's default choice, since I have found this error already multiple times.

namespace TwoPhaseFlowWithPP
{
/** Create a porosity model
* @param config ConfigTree object has a tag of <porosity>
Copy link
Member

Choose a reason for hiding this comment

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

Is this indeed creating a porosity model?

class TwoPhaseFlowWithPPMaterialProperties
{
public:
typedef MaterialLib::Fluid::FluidProperty::ArrayType ArrayType;
Copy link
Member

Choose a reason for hiding this comment

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

using is preferred over typedef

void setMaterialID(const ProcessLib::SpatialPosition& pos);

/**
* \brief Compute the coefficient of the mass term by
Copy link
Member

Choose a reason for hiding this comment

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

Is this correct documentation for getPermeability?

std::unique_ptr<MaterialLib::Fluid::FluidProperty> _gas_density;
std::unique_ptr<MaterialLib::Fluid::FluidProperty> _gas_viscosity;

/** Use porous medium models for different material zones.
Copy link
Member

Choose a reason for hiding this comment

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

The "porous medium models" is from different class. Need to change this to two phase models. (Not only here, but in all similar places.)

// Process variable.
auto process_variables = findProcessVariables(
variables, config,
{//! \ogs_file_param_special{process__TWOPHASE_FLOW__WITHPP__process_variables__process_variable}
Copy link
Member

@endJunction endJunction Nov 10, 2016

Choose a reason for hiding this comment

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

One underscore too much in TWOPHASE_FLOW__WITHPP.
Actually the name is wrong. Please use "TWOPHASE_FLOW_PP" here and below.

Eigen::VectorXd specific_body_force;

std::vector<double> const b =
//! \ogs_file_param_special{process__TWOPHASE_FLOW_specific_body_force}
Copy link
Member

Choose a reason for hiding this comment

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

needs an underscore before 'specific'

#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
#include "NumLib/Function/Interpolation.h"
#include "TwoPhaseFlowWithPPProcessData.h"
// using namespace Eigen;
Copy link
Member

Choose a reason for hiding this comment

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

rm

#ifndef OGS_TWOPHASEFLOWWITHPPLOCALASSEMBLER_IMPL_H
#define OGS_TWOPHASEFLOWWITHPPLOCALASSEMBLER_IMPL_H

#include <iostream>
Copy link
Member

Choose a reason for hiding this comment

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

is this needed?


// Note: currently only isothermal case is considered, so the temperature is
// assumed to be const
// the variation of temperatura will be taken into account in future
Copy link
Member

Choose a reason for hiding this comment

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

formatting. typo in temperature

sm.dNdx.transpose() * permeability * sm.dNdx * integration_factor;
Klpc.noalias() +=
K_mat_coeff(cap_pressure_coeff_index, cap_pressure_coeff_index) *
sm.dNdx.transpose() * permeability * sm.dNdx * integration_factor;
Copy link
Member

Choose a reason for hiding this comment

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

common expressions should be calculated once.

@endJunction
Copy link
Member

@Yonghui56 You have to answer the comments or change the code. So long I'll set this PR to "needs an update". When ready you can set it to "please review".

@ogsbot
Copy link
Member

ogsbot commented Nov 16, 2016

Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1721/

@ogsbot
Copy link
Member

ogsbot commented Nov 16, 2016

Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1722/

@@ -65,6 +64,9 @@ void TwoPhaseFlowWithPPLocalAssembler<
auto Mlpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
cap_pressure_matrix_index, cap_pressure_matrix_index);

NodalMatrixType laplace_factor;
laplace_factor.setZero(nonwet_pressure_size, cap_pressure_size);
Copy link
Member

Choose a reason for hiding this comment

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

I think the NodalMatrixType is not of size nonwet_pressure_size times cap_pressure_size, isn't it?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@endJunction NodalMatrixType holds the size of [ShapeFunction::NPOINTS] times [ShapeFunction::NPOINTS], while nonwet_pressure_size is equal to ShapeFunction::NPOINTS, so is cap_pressure_size.

// nonwetting
double const drhogas_dpg = _process_data._material->getDerivGasDensity(
pg_int_pt, _temperature);
M_mat_coeff(nonwet_pressure_coeff_index, nonwet_pressure_coeff_index) =
Copy link
Member

@wenqing wenqing Nov 17, 2016

Choose a reason for hiding this comment

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

I think matrix type of M_mat_coeff is not needed. You can put the coefficient directly to the matrix calculation like

Mgp.noalias() +=    (porosity * (1 - Sw) * drhogas_dp ) *
+                         sm.N.transpose() * sm.N * integration_factor;

double const lambda_L = k_rel_L / mu_liquid;
K_mat_coeff(cap_pressure_coeff_index, nonwet_pressure_coeff_index) =
rho_w * lambda_L;
K_mat_coeff(cap_pressure_coeff_index, cap_pressure_coeff_index) =
Copy link
Member

Choose a reason for hiding this comment

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

Similarly, to have a matrix of K_mat_coeff is not necessary.

std::move(secondary_variables), std::move(named_function_caller)),
_process_data(std::move(process_data))
{
DBUG("Create Two-Phase flow process.");
Copy link
Member

Choose a reason for hiding this comment

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

Create two-thase flow process with primary variables of capillary pressure and gas pressure, or simplyCreate TwoPhaseFlowWithPPProcess``

@@ -0,0 +1,99 @@
/**
Copy link
Member

Choose a reason for hiding this comment

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

Since files in TwoPhaseModels are for a set of process related parameters, I would like to move them to ProcessLib/TwoPhaseFlowWithPP.

In MaterialLib, I think, it should only has physical chemical models or parameters of fluids, solids and porous medium.

Copy link
Member

@wenqing wenqing left a comment

Choose a reason for hiding this comment

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

Looks good

EXECUTABLE ogs
EXECUTABLE_ARGS Twophase_Lia_quad2.prj
TESTER vtkdiff
ABSTOL 1e-1 RELTOL 1e-1
Copy link
Member

Choose a reason for hiding this comment

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

Why is the tolerance so loose? On my computer I can set 1e-9 and 1e-12 for the absolute and relative tolerances.

Copy link
Member

@endJunction endJunction Nov 17, 2016

Choose a reason for hiding this comment

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

You have also the gas_pressure with 1e-8 and 1e-13, and saturation with 1e-14 and 1e-14
absolute and relative tolerances (on my computer -- might be slightly looser on other machines), which you can add to the ctest; Just another lines with

Lia_20.vtu twophaseflow_pcs_0_ts_200_t_20.000000.vtu gas_pressure gas_pressure

taking 1e-8 and 1e-12 as the tolerances.

@endJunction
Copy link
Member

Squash after adding additional ctest vtkdiff lines; I'll merge then.

@wenqing
Copy link
Member

wenqing commented Nov 17, 2016

Calculation of velocity and saturation as secondary variables can be added later on.

@Yonghui56 Yonghui56 force-pushed the pr_twophasewithpp_curve_test branch 5 times, most recently from bdf8afe to 538306d Compare November 18, 2016 13:25
@Yonghui56 Yonghui56 force-pushed the pr_twophasewithpp_curve_test branch 2 times, most recently from 61cd686 to 657a074 Compare November 18, 2016 14:43
@endJunction
Copy link
Member

Good! When the tests are finished I'll merge it.

@endJunction endJunction merged commit d15070c into ufz:master Nov 20, 2016
@ogsbot
Copy link
Member

ogsbot commented Jun 19, 2020

OpenGeoSys development has been moved to GitLab.

See this pull request on GitLab.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants