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 two phase flow process with total mass density model(WIP) #1613

Merged
merged 3 commits into from Feb 6, 2017

Conversation

Projects
None yet
3 participants
@Yonghui56
Contributor

Yonghui56 commented Dec 14, 2016

As titled, this pr creates a process aiming to handle two phase multi-component process. The wetting phase pressure and total mass density are selected as primary variables. Local constitutive relations are based on Henry law and Ideal gas law and local Newton is applied for the nonlinearity. Regularized van Genuchten PC-S model is complemented to the original vG model.
(@endJunction could you please help me to add a WIP label)

@Yonghui56 Yonghui56 force-pushed the Yonghui56:pr_2phase_totalmassdensity_pl_wip branch 2 times, most recently from 06bb8b8 to 95020cf Dec 14, 2016

@endJunction

This comment has been minimized.

Member

endJunction commented Dec 14, 2016

Nice, first time it compiles on all platforms 🎆

@Yonghui56 Yonghui56 force-pushed the Yonghui56:pr_2phase_totalmassdensity_pl_wip branch from 0012e6e to b57b470 Jan 4, 2017

@endJunction

This comment has been minimized.

Member

endJunction commented Jan 6, 2017

This needs a rebase due to changes in the test system.

@endJunction

This comment has been minimized.

Member

endJunction commented Jan 9, 2017

Please update the copyright years to 2012-2017.

@Yonghui56 Yonghui56 force-pushed the Yonghui56:pr_2phase_totalmassdensity_pl_wip branch 3 times, most recently from 9877a9c to b61b9d5 Jan 16, 2017

*/
#ifndef OGS_CREATETWOPHASEFLOWPRHOMATERIALPROPERTIES_H
#define OGS_CREATETWOPHASEFLOWPRHOMATERIALPROPERTIES_H

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

Please replace the include guards with #pragma once

@@ -74,6 +74,8 @@ const double O2 = 0.032; ///< kg \per{mol}
* - Source: http://www.engineeringtoolbox.com/molecular-mass-air-d_679.html
*/
const double Air = 0.02897; ///< kg \per{mol}
const double H2 = 0.002;

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

From https://pubchem.ncbi.nlm.nih.gov/compound/Hydrogen it is 2.016 g/mol.
What is your reference?

@@ -44,14 +44,15 @@ class BrookCoreyCapillaryPressureSaturation final
/**
* @param pb Entry pressure, \f$ p_b \f$
* @param Sr Residual saturation, \f$ S_r \f$
* @param Sg_r Residual saturation of gas phase, \f$ Sg_r \f$

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

Formatting: Please use no tabs.

DBUG(
"Using value %g for nonwetting phase residual saturation in "
"capillary pressure "
"model.",

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

Formatting

const double Se = (S - _saturation_r) / (_saturation_max - _saturation_r);
const double pc = _pb * std::pow(std::pow(Se, (-1.0 / _m)) - 1.0, 1.0 - _m);
return MathLib::limitValueInInterval(pc, _minor_offset, _pc_max);
if (_has_Regularized)

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

variable name should not contain capital letters.

}
else if (Sg < 0.0)
{
return getPc_bar_vG_Sg(0.0) + get_dPCdS_vG_bar(0.0) * (Sg - 0.0);

This comment has been minimized.

@endJunction
* derivative dPCdS based on regularized van Genuchten capillary
* pressure-saturation Model
*/
double get_dPCdS_vG_bar(double Sg) const

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

function names not consistent: get_dx but getx.

double get_dPCdS_vG_bar(double Sg) const
{
double S_bar = 0.0;
S_bar = getS_bar(Sg);

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

You can directly initialize S_bar with the required value. Make it const afterwards.

@@ -334,4 +334,4 @@ createNonlinearSolver(GlobalLinearSolver& linear_solver,
}
OGS_FATAL("Unsupported nonlinear solver type");
}
}
}

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

Please avoid unrelated whitespace changes.

@@ -119,7 +119,8 @@ TESProcess::TESProcess(
config.getConfigParameterOptional<bool>("output_element_matrices"))
{
DBUG("output_element_matrices: %s", (*param) ? "true" : "false");
bool test = false;
test= *param;

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

Why is this needed?

auto const& capillary_pressure_conf =
//! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PRHO__material_property__porous_medium__porous_medium__capillary_pressure}
conf.getConfigSubtree("capillary_pressure");
auto pc = MaterialLib::PorousMedium::createCapillaryPressureModel(

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

Can the results of createXxModel be const?

This comment has been minimized.

@Yonghui56

Yonghui56 Jan 27, 2017

Contributor

No. The details maybe @wenqing can answer

This comment has been minimized.

@wenqing

wenqing Feb 1, 2017

Member

Because of std::move(pc) in line 102, it can not be const.

//! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PRHO__specific_body_force}
config.getConfigParameter<std::vector<double>>("specific_body_force");
assert(b.size() > 0 && b.size() < 4);
specific_body_force.resize(b.size());

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

You can create the specific_body_force vector with correct size. See documentation of the available constructors http://eigen.tuxfamily.org/dox/group__QuickRefPage.html

This comment has been minimized.

@Yonghui56

Yonghui56 Jan 27, 2017

Contributor

As far as I know, the correct size should be consistent to the dimension of the mesh. At this point, the correct size of mesh cannot be obtained.

This comment has been minimized.

@endJunction

endJunction Jan 27, 2017

Member

The size of specific body force is given by the config parameter (b here). So you can write at line 64

Eigen::VectorXd const specific_body_force(b.size());

and remove the definition above.

Update: Maybe const is not going to work...

This comment has been minimized.

@Yonghui56

Yonghui56 Jan 27, 2017

Contributor

Yes.

std::copy_n(b.data(), b.size(), specific_body_force.data());
//! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PRHO__mass_lumping}
auto mass_lump = config.getConfigParameter<bool>("mass_lumping");

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

mass_lump -> mass_lumping.

config,
//! \ogs_file_param_special{prj__processes__process__TWOPHASE_FLOW_PRHO__diffusion_coeff_componentb}
"diffusion_coeff_componentb", parameters, 1);
auto& diff_coeff_a = findParameter<double>(

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

Sometimes you write ...coeff_x and sometimes ...coeffx. Be nice to your readers, please!

This comment has been minimized.

@endJunction

endJunction Jan 31, 2017

Member

still, some of the identifiers are coeff_x and some coeffx.

_process_data._material->setMaterialID(pos);
const Eigen::MatrixXd& perm = _process_data._material->getPermeability(
t, pos, _element.getDimension());

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

Just curious: Should it be element's dimension or the GlobalDim value?

This comment has been minimized.

@Yonghui56

Yonghui56 Jan 27, 2017

Contributor

I think it should be element's dimension, regarding this, @wenqing has made a comment in Liquid Flow process.

This comment has been minimized.

@endJunction

endJunction Jan 27, 2017

Member

Could you tell us where the comment is, please?

This comment has been minimized.

_process_data._material->getLiquidDensity(pl_int_pt, _temperature);
double& Sw = _ip_data[ip]._sw;
double const X_h2_nonwet = 1.0; // TODO

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

Please elaborate what is todo.

NumLib::shapeFunctionInterpolate(local_x, sm.N, pl_int_pt,
totalrho_int_pt);
auto const& wp = _integration_method.getWeightedPoint(ip);
const double integration_factor =

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

You already use integration point data structure, why not put all other constant per-integration-point-values there?

double const porosity = _process_data._material->getPorosity(
t, pos, pl_int_pt, _temperature, 0);
mass_operator.noalias() = sm.N.transpose() * sm.N * integration_factor;

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

(cont'd from integration point data comment)
... especialy those which are expensive to compute...

double const mu_gas = _process_data._material->getGasViscosity(
_pressure_nonwetting[ip], _temperature);
double const lambda_G = k_rel_G / mu_gas;
double const diffusion_coeff_componenth2 =

This comment has been minimized.

@endJunction

endJunction Jan 20, 2017

Member

is it ..._component2 or is it ..._component_h2?

* John Wiley & Sons, 2012.
*/
double const Henry_const = 7.65e-6; // mol/Pa./m3
double const& molar_mass_h2 = MaterialLib::PhysicalConstant::MolarMass::H2;

This comment has been minimized.

@wenqing

wenqing Feb 1, 2017

Member

to using molar_mass_h2 = MaterialLib::PhysicalConstant::MolarMass::H2;

This comment has been minimized.

@wenqing

wenqing Feb 1, 2017

Member

name starts with _

*/
double const Henry_const = 7.65e-6; // mol/Pa./m3
double const& molar_mass_h2 = MaterialLib::PhysicalConstant::MolarMass::H2;
double const& R = MaterialLib::PhysicalConstant::IdealGasConstant;

This comment has been minimized.

@wenqing

wenqing Feb 1, 2017

Member

--> using R = MaterialLib::PhysicalConstant::IdealGasConstant;

This comment has been minimized.

@endJunction

endJunction Feb 1, 2017

Member

@wenqing I think one can write

using MaterialLib::PhysicalConstant::IdealGasConstant;

and then the IdealGasConstant is available in the current scope.

This comment has been minimized.

@wenqing

wenqing Feb 1, 2017

Member

@endJunction Nice, better in in cpp file only.

@Yonghui56 Yonghui56 force-pushed the Yonghui56:pr_2phase_totalmassdensity_pl_wip branch 2 times, most recently from f3d06c1 to f7a3a5f Feb 1, 2017

@wenqing

wenqing approved these changes Feb 1, 2017

Looks good

* De Nevers N. Physical and chemical equilibrium for chemical engineers[M].
* John Wiley & Sons, 2012.
*/
double const HenryConstH2 = 7.65e-6; /// mol/Pa./m3

This comment has been minimized.

@endJunction

endJunction Feb 2, 2017

Member

Either write HenryConstantH2 or drop the HenryConst (it is already indicated by the namespace).

This comment has been minimized.

@Yonghui56

Yonghui56 Feb 2, 2017

Contributor

HenryConstantH2 is preferred.

"Using value %g for nonwetting phase residual saturation in "
"capillary pressure model.",
(*Sg_res));
Sg_r = *Sg_res;

This comment has been minimized.

@endJunction

endJunction Feb 2, 2017

Member

What is the difference between Sg_r sg_res and Sg_res?
And later there is an s_gr....

This comment has been minimized.

@Yonghui56

Yonghui56 Feb 2, 2017

Contributor

actually Sg_res is a pointer. I change Sg_res to Sg_r_ptr. Others have been changed to Sg_r uniformly

/// Regularized van Genuchten capillary pressure-saturation Model
double getSBar(double Sg) const;
/// van Genuchten capillary pressure-saturation Model
double getPcvGSg(double Sg) const;

This comment has been minimized.

@endJunction
@endJunction

Few (I hope now) last changes needed.
When you have finished all the changes, please create few distinct commits, then we can merge it.

@Yonghui56 Yonghui56 force-pushed the Yonghui56:pr_2phase_totalmassdensity_pl_wip branch from 9077ce7 to cf467ed Feb 2, 2017

@endJunction

This comment has been minimized.

Member

endJunction commented Feb 2, 2017

  • In ufz/ogs-data git diff --check origin/master..origin/twophase_prho_flow fails. You need to pay attention to trailing whitespaces (and possibly tabs).
  • Members and their initialization order do not match: TwoPhaseFlowWithPrhoLocalAssembler.h:30:11: warning: field 'mat_property' will be initialized after field 'sw'
  • Unused variables sw_pre and rho_m_pre and possibly others should be removed.
  • I could not find the ctest(s); Am I missing something? Update: Large one ok. Small is needed.
  • I didn't find the documentation on docs.opengeosys.org (it might be I'm looking in wrong place.)

@Yonghui56 Yonghui56 force-pushed the Yonghui56:pr_2phase_totalmassdensity_pl_wip branch 9 times, most recently from adcaaec to 1caf985 Feb 2, 2017

@Yonghui56 Yonghui56 force-pushed the Yonghui56:pr_2phase_totalmassdensity_pl_wip branch from 1caf985 to bbdee08 Feb 5, 2017

@Yonghui56 Yonghui56 force-pushed the Yonghui56:pr_2phase_totalmassdensity_pl_wip branch from 9190853 to e5b257f Feb 5, 2017

@endJunction

This comment has been minimized.

Member

endJunction commented Feb 6, 2017

To be merged after lunch.

@wenqing

This comment has been minimized.

Member

wenqing commented Feb 6, 2017

@endJunction endJunction merged commit d4acfe7 into ufz:master Feb 6, 2017

3 checks passed

continuous-integration/jenkins/pr-merge This commit looks good
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
ufz/ogs/PR-1613 Build #107 succeeded in 33 min
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment