From 1a73766e7cb75bb6c5164b5cfa09b27c8dce316c Mon Sep 17 00:00:00 2001 From: Yonghui <huangyh56@gmail.com> Date: Mon, 23 Jan 2017 00:52:10 +0100 Subject: [PATCH] Implement the material properties in PP model --- Applications/ApplicationsLib/ProjectData.cpp | 2 +- .../i_relative_permeability.md | 1 + .../relative_permeability/a_id.md | 1 + .../i_relative_permeability.md | 1 + .../porous_medium/t_capillary_pressure.md | 1 + .../material_property/t_gas_density.md | 1 + .../material_property/t_gas_viscosity.md | 1 + .../material_property/t_liquid_density.md | 1 + .../material_property/t_liquid_viscosity.md | 1 + .../TWOPHASE_FLOW_PP/t_mass_lumping.md | 2 +- .../TWOPHASE_FLOW_PP/t_specific_body_force.md | 2 +- .../process/TWOPHASE_FLOW_PP/t_temperature.md | 2 + .../PiecewiseLinearInterpolation.cpp | 2 +- ProcessLib/TwoPhaseFlowWithPP/CMakeLists.txt | 67 +++++-- ...teTwoPhaseFlowWithPPMaterialProperties.cpp | 139 ++++++++++++++ ...eateTwoPhaseFlowWithPPMaterialProperties.h | 30 ++++ .../CreateTwoPhaseFlowWithPPProcess.cpp | 82 +++------ .../CreateTwoPhaseFlowWithPPProcess.h | 2 +- .../TwoPhaseFlowWithPPLocalAssembler-impl.h | 170 +++++++++--------- .../TwoPhaseFlowWithPPLocalAssembler.h | 79 ++++++-- .../TwoPhaseFlowWithPPMaterialProperties.cpp | 169 +++++++++++++++++ .../TwoPhaseFlowWithPPMaterialProperties.h | 129 +++++++++++++ .../TwoPhaseFlowWithPPProcess.cpp | 6 +- .../TwoPhaseFlowWithPPProcess.h | 1 - .../TwoPhaseFlowWithPPProcessData.h | 48 ++--- 25 files changed, 729 insertions(+), 211 deletions(-) create mode 100644 Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/i_relative_permeability.md create mode 100644 Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/a_id.md create mode 100644 Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/i_relative_permeability.md create mode 100644 Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/t_capillary_pressure.md create mode 100644 Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_density.md create mode 100644 Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_viscosity.md create mode 100644 Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_density.md create mode 100644 Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_viscosity.md create mode 100644 Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_temperature.md create mode 100644 ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.cpp create mode 100644 ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.h create mode 100644 ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.cpp create mode 100644 ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.h diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp index 2436b1a083b..d3d46f9b442 100644 --- a/Applications/ApplicationsLib/ProjectData.cpp +++ b/Applications/ApplicationsLib/ProjectData.cpp @@ -437,7 +437,7 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config, else if (type == "TWOPHASE_FLOW_PP") { process = - ProcessLib::TwoPhaseFlowWithPP::CreateTwoPhaseFlowWithPPProcess( + ProcessLib::TwoPhaseFlowWithPP::createTwoPhaseFlowWithPPProcess( *_mesh_vec[0], std::move(jacobian_assembler), _process_variables, _parameters, integration_order, process_config, _curves); diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/i_relative_permeability.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/i_relative_permeability.md new file mode 100644 index 00000000000..e3e0c731598 --- /dev/null +++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/i_relative_permeability.md @@ -0,0 +1 @@ +A tag for the relative permeability model, which holds two ids: one is for the Nonwetting phase and the other is for the Wetting phase. diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/a_id.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/a_id.md new file mode 100644 index 00000000000..a64a7806df4 --- /dev/null +++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/a_id.md @@ -0,0 +1 @@ +A tag for the relative permeability model id: 0 is Nonwetting phase and 1 is Wetting phase. diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/i_relative_permeability.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/i_relative_permeability.md new file mode 100644 index 00000000000..ea0efd1bb55 --- /dev/null +++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/i_relative_permeability.md @@ -0,0 +1 @@ +Defines the relative permeability model. diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/t_capillary_pressure.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/t_capillary_pressure.md new file mode 100644 index 00000000000..19d854870d0 --- /dev/null +++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/t_capillary_pressure.md @@ -0,0 +1 @@ +Capillary pressure model. diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_density.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_density.md new file mode 100644 index 00000000000..bcf97b55a9e --- /dev/null +++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_density.md @@ -0,0 +1 @@ +Density model of the gas phase. diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_viscosity.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_viscosity.md new file mode 100644 index 00000000000..a02e98700f6 --- /dev/null +++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_viscosity.md @@ -0,0 +1 @@ +Viscosity model of the gas phase. diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_density.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_density.md new file mode 100644 index 00000000000..5c8db5dec3d --- /dev/null +++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_density.md @@ -0,0 +1 @@ +Density model of the liquid phase. diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_viscosity.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_viscosity.md new file mode 100644 index 00000000000..0451570185b --- /dev/null +++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_viscosity.md @@ -0,0 +1 @@ +Viscosity model of the liquid phase. diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_mass_lumping.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_mass_lumping.md index 0d6032b0aad..2cfa81d1102 100644 --- a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_mass_lumping.md +++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_mass_lumping.md @@ -1 +1 @@ -\copydoc ProcessLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPProcessData::_has_mass_lumping +\copydoc ProcessLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPProcessData::has_mass_lumping diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_specific_body_force.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_specific_body_force.md index 7437fbf0be7..c1658dbc4dd 100644 --- a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_specific_body_force.md +++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_specific_body_force.md @@ -1 +1 @@ -\copydoc ProcessLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPProcessData::_specific_body_force +\copydoc ProcessLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPProcessData::specific_body_force diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_temperature.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_temperature.md new file mode 100644 index 00000000000..ee031866715 --- /dev/null +++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_temperature.md @@ -0,0 +1,2 @@ +An input of a reference temperature (K) for two-phase flow process. +It is fixed to be constant for the simulation due to an isothermal assumption. diff --git a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp index 4fec226af6f..df5ac0e9b9a 100644 --- a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp +++ b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp @@ -95,7 +95,7 @@ double PiecewiseLinearInterpolation::getDerivative( interval_idx = 1; } - if (interval_idx > 2 && interval_idx < _supp_pnts.size() - 1) + if (interval_idx > 1 && interval_idx < _supp_pnts.size() - 2) { // left and right support points. double const x_ll = _supp_pnts[interval_idx - 2]; diff --git a/ProcessLib/TwoPhaseFlowWithPP/CMakeLists.txt b/ProcessLib/TwoPhaseFlowWithPP/CMakeLists.txt index b21e3b6008d..12a9b699452 100644 --- a/ProcessLib/TwoPhaseFlowWithPP/CMakeLists.txt +++ b/ProcessLib/TwoPhaseFlowWithPP/CMakeLists.txt @@ -1,26 +1,69 @@ AddTest( - NAME 2D_TwoPhase_PP_Lia_quad + NAME 2D_TwoPhase_PP_Lia_quad_1 PATH Parabolic/TwoPhaseFlowPP/Liakopoulos EXECUTABLE ogs - EXECUTABLE_ARGS Twophase_Lia_quad2_short.prj + EXECUTABLE_ARGS TwoPhase_Lia_quad_short.prj TESTER vtkdiff REQUIREMENTS NOT OGS_USE_MPI - ABSTOL 1e-8 RELTOL 1e-12 + ABSTOL 1e-2 RELTOL 1e-4 DIFF_DATA - Lia_20.vtu twophaseflow_pcs_0_ts_119_t_20.000000.vtu capillary_pressure capillary_pressure - Lia_20.vtu twophaseflow_pcs_0_ts_119_t_20.000000.vtu gas_pressure gas_pressure - Lia_20.vtu twophaseflow_pcs_0_ts_119_t_20.000000.vtu saturation saturation + h2_Liako_20.vtu twophaseflow_pcs_0_ts_218_t_20.000000.vtu saturation saturation ) AddTest( - NAME LARGE_2D_TwoPhase_PP_Lia_quad + NAME 2D_TwoPhase_PP_Lia_quad_2 PATH Parabolic/TwoPhaseFlowPP/Liakopoulos EXECUTABLE ogs - EXECUTABLE_ARGS Twophase_Lia_quad2_large.prj + EXECUTABLE_ARGS TwoPhase_Lia_quad_short.prj TESTER vtkdiff REQUIREMENTS NOT OGS_USE_MPI - ABSTOL 1e-8 RELTOL 1e-12 + ABSTOL 20 RELTOL 1e-3 DIFF_DATA - Lia_1000.vtu twophaseflow_pcs_0_ts_1180_t_1000.000000.vtu capillary_pressure capillary_pressure - Lia_1000.vtu twophaseflow_pcs_0_ts_1180_t_1000.000000.vtu gas_pressure gas_pressure - Lia_1000.vtu twophaseflow_pcs_0_ts_1180_t_1000.000000.vtu saturation saturation + h2_Liako_20.vtu twophaseflow_pcs_0_ts_218_t_20.000000.vtu capillary_pressure capillary_pressure + h2_Liako_20.vtu twophaseflow_pcs_0_ts_218_t_20.000000.vtu gas_pressure gas_pressure +) +AddTest( + NAME LARGE_2D_TwoPhase_PP_Lia_quad_1 + PATH Parabolic/TwoPhaseFlowPP/Liakopoulos + EXECUTABLE ogs + EXECUTABLE_ARGS TwoPhase_Lia_quad_large.prj + TESTER vtkdiff + REQUIREMENTS NOT OGS_USE_MPI + ABSTOL 1e-2 RELTOL 1e-3 + DIFF_DATA + h2_Liako_1198.vtu twophaseflow_pcs_0_ts_1198_t_1000.000000.vtu SATURATION1 saturation +) +AddTest( + NAME LARGE_2D_TwoPhase_PP_Lia_quad_2 + PATH Parabolic/TwoPhaseFlowPP/Liakopoulos + EXECUTABLE ogs + EXECUTABLE_ARGS TwoPhase_Lia_quad_large.prj + TESTER vtkdiff + REQUIREMENTS NOT OGS_USE_MPI + ABSTOL 20 RELTOL 1e-2 + DIFF_DATA + h2_Liako_1198.vtu twophaseflow_pcs_0_ts_1198_t_1000.000000.vtu PRESSURE1 capillary_pressure + h2_Liako_1198.vtu twophaseflow_pcs_0_ts_1198_t_1000.000000.vtu PRESSURE2 gas_pressure +) +AddTest( + NAME 1D_TwoPhase_PP_mcwt_1 + PATH Parabolic/TwoPhaseFlowPP/McWorts + EXECUTABLE ogs + EXECUTABLE_ARGS TwoPhase_mcwt_line.prj + TESTER vtkdiff + REQUIREMENTS NOT OGS_USE_MPI + ABSTOL 1e-2 RELTOL 1e-4 + DIFF_DATA + mcwt_1000.vtu twophaseflow_pcs_0_ts_519_t_1000.000000.vtu SATURATION1 saturation +) +AddTest( + NAME 1D_TwoPhase_PP_mcwt_2 + PATH Parabolic/TwoPhaseFlowPP/McWorts + EXECUTABLE ogs + EXECUTABLE_ARGS TwoPhase_mcwt_line.prj + TESTER vtkdiff + REQUIREMENTS NOT OGS_USE_MPI + ABSTOL 10 RELTOL 1e-3 + DIFF_DATA + mcwt_1000.vtu twophaseflow_pcs_0_ts_519_t_1000.000000.vtu PRESSURE1 capillary_pressure + mcwt_1000.vtu twophaseflow_pcs_0_ts_519_t_1000.000000.vtu PRESSURE2 gas_pressure ) diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.cpp b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.cpp new file mode 100644 index 00000000000..5c5f353d016 --- /dev/null +++ b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.cpp @@ -0,0 +1,139 @@ +/** + * \copyright + * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#include "CreateTwoPhaseFlowWithPPMaterialProperties.h" +#include <logog/include/logog.hpp> +#include "BaseLib/reorderVector.h" +#include "MaterialLib/Fluid/FluidProperty.h" +#include "MaterialLib/PorousMedium/Porosity/Porosity.h" +#include "MaterialLib/PorousMedium/Storage/Storage.h" +#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h" +#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h" +#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h" +#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h" +#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h" +#include "MeshLib/Mesh.h" +#include "MeshLib/PropertyVector.h" +#include "ProcessLib/Parameter/Parameter.h" +#include "ProcessLib/Parameter/SpatialPosition.h" +#include "TwoPhaseFlowWithPPMaterialProperties.h" + +namespace ProcessLib +{ +namespace TwoPhaseFlowWithPP +{ +std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties> +createTwoPhaseFlowWithPPMaterialProperties( + BaseLib::ConfigTree const& config, + boost::optional<MeshLib::PropertyVector<int> const&> + material_ids) +{ + DBUG("Reading material properties of two-phase flow process."); + + //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__fluid} + auto const& fluid_config = config.getConfigSubtree("fluid"); + + // Get fluid properties + //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__liquid_density} + auto const& rho_conf = fluid_config.getConfigSubtree("liquid_density"); + auto liquid_density = MaterialLib::Fluid::createFluidDensityModel(rho_conf); + //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__gas_density} + auto const& rho_gas_conf = fluid_config.getConfigSubtree("gas_density"); + auto gas_density = + MaterialLib::Fluid::createFluidDensityModel(rho_gas_conf); + //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__liquid_viscosity} + auto const& mu_conf = fluid_config.getConfigSubtree("liquid_viscosity"); + auto liquid_viscosity = MaterialLib::Fluid::createViscosityModel(mu_conf); + //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__gas_viscosity} + auto const& mu_gas_conf = fluid_config.getConfigSubtree("gas_viscosity"); + auto gas_viscosity = MaterialLib::Fluid::createViscosityModel(mu_gas_conf); + + // Get porous properties + std::vector<int> mat_ids; + std::vector<int> mat_krel_ids; + std::vector<Eigen::MatrixXd> intrinsic_permeability_models; + std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>> + porosity_models; + std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>> + storage_models; + std::vector< + std::unique_ptr<MaterialLib::PorousMedium::CapillaryPressureSaturation>> + capillary_pressure_models; + std::vector< + std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>> + relative_permeability_models; + + //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium} + auto const& poro_config = config.getConfigSubtree("porous_medium"); + //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium} + for (auto const& conf : poro_config.getConfigSubtreeList("porous_medium")) + { + //! \ogs_file_attr{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__id} + auto const id = conf.getConfigAttributeOptional<int>("id"); + mat_ids.push_back(*id); + + //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__permeability} + auto const& permeability_conf = conf.getConfigSubtree("permeability"); + intrinsic_permeability_models.emplace_back( + MaterialLib::PorousMedium::createPermeabilityModel( + permeability_conf)); + + //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__porosity} + auto const& porosity_conf = conf.getConfigSubtree("porosity"); + auto n = MaterialLib::PorousMedium::createPorosityModel(porosity_conf); + porosity_models.emplace_back(std::move(n)); + + //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__storage} + auto const& storage_conf = conf.getConfigSubtree("storage"); + auto beta = MaterialLib::PorousMedium::createStorageModel(storage_conf); + storage_models.emplace_back(std::move(beta)); + + auto const& capillary_pressure_conf = + //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__capillary_pressure} + conf.getConfigSubtree("capillary_pressure"); + auto pc = MaterialLib::PorousMedium::createCapillaryPressureModel( + capillary_pressure_conf); + capillary_pressure_models.emplace_back(std::move(pc)); + + auto const& krel_config = + //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__relative_permeability} + conf.getConfigSubtree("relative_permeability"); + for ( + auto const& krel_conf : + //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__relative_permeability__relative_permeability} + krel_config.getConfigSubtreeList("relative_permeability")) + { + auto const krel_id = + //! \ogs_file_attr{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__relative_permeability__relative_permeability__id} + krel_conf.getConfigAttributeOptional<int>("id"); + mat_krel_ids.push_back(*krel_id); + auto krel_n = + MaterialLib::PorousMedium::createRelativePermeabilityModel( + krel_conf); + relative_permeability_models.emplace_back(std::move(krel_n)); + } + BaseLib::reorderVector(relative_permeability_models, mat_krel_ids); + } + + BaseLib::reorderVector(intrinsic_permeability_models, mat_ids); + BaseLib::reorderVector(porosity_models, mat_ids); + BaseLib::reorderVector(storage_models, mat_ids); + + return std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>{ + new TwoPhaseFlowWithPPMaterialProperties{ + material_ids, std::move(liquid_density), + std::move(liquid_viscosity), std::move(gas_density), + std::move(gas_viscosity), intrinsic_permeability_models, + std::move(porosity_models), std::move(storage_models), + std::move(capillary_pressure_models), + std::move(relative_permeability_models)}}; +} + +} // end namespace +} // end namespace diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.h b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.h new file mode 100644 index 00000000000..a00d066f382 --- /dev/null +++ b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.h @@ -0,0 +1,30 @@ +/** + * \copyright + * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org) + * Distributed under a Modified BSD License. + * See accompanying file LICENSE.txt or + * http://www.opengeosys.org/project/license + * + */ + +#pragma once + + +#include "TwoPhaseFlowWithPPMaterialProperties.h" +namespace BaseLib +{ +class ConfigTree; +} + +namespace ProcessLib +{ +namespace TwoPhaseFlowWithPP +{ +std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties> +createTwoPhaseFlowWithPPMaterialProperties( + BaseLib::ConfigTree const& config, + boost::optional<MeshLib::PropertyVector<int> const&> + material_ids); + +} // end namespace +} // end namespace diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp index 460edafca73..03375787171 100644 --- a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp +++ b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp @@ -9,20 +9,20 @@ #include "CreateTwoPhaseFlowWithPPProcess.h" #include <cassert> -#include "MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.h" #include "MeshLib/MeshGenerators/MeshGenerator.h" #include "ProcessLib/Parameter/ConstantParameter.h" #include "ProcessLib/Utils/ParseSecondaryVariables.h" #include "ProcessLib/Utils/ProcessUtils.h" +#include "CreateTwoPhaseFlowWithPPMaterialProperties.h" +#include "TwoPhaseFlowWithPPMaterialProperties.h" #include "TwoPhaseFlowWithPPProcess.h" #include "TwoPhaseFlowWithPPProcessData.h" - namespace ProcessLib { namespace TwoPhaseFlowWithPP { -std::unique_ptr<Process> CreateTwoPhaseFlowWithPPProcess( +std::unique_ptr<Process> createTwoPhaseFlowWithPPProcess( MeshLib::Mesh& mesh, std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler, std::vector<ProcessVariable> const& variables, @@ -37,9 +37,6 @@ std::unique_ptr<Process> CreateTwoPhaseFlowWithPPProcess( config.checkConfigParameter("type", "TWOPHASE_FLOW_PP"); DBUG("Create TwoPhaseFlowProcess with PP model."); - - // Process variable. - //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__process_variables} auto const pv_config = config.getConfigSubtree("process_variables"); @@ -58,81 +55,48 @@ std::unique_ptr<Process> CreateTwoPhaseFlowWithPPProcess( ProcessLib::parseSecondaryVariables(config, secondary_variables, named_function_caller); // Specific body force - Eigen::VectorXd specific_body_force; - std::vector<double> const b = //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__specific_body_force} config.getConfigParameter<std::vector<double>>("specific_body_force"); assert(b.size() > 0 && b.size() < 4); - specific_body_force.resize(b.size()); + Eigen::VectorXd specific_body_force(b.size()); bool const has_gravity = MathLib::toVector(b).norm() > 0; if (has_gravity) std::copy_n(b.data(), b.size(), specific_body_force.data()); //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__mass_lumping} - auto mass_lump = config.getConfigParameter<bool>("mass_lumping"); + auto const mass_lumping = config.getConfigParameter<bool>("mass_lumping"); + + auto& temperature = findParameter<double>( + config, + //! \ogs_file_param_special{prj__processes__process__TWOPHASE_FLOW_PP__temperature} + "temperature", parameters, 1); //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property} auto const& mat_config = config.getConfigSubtree("material_property"); - std::unique_ptr< - MaterialLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPMaterialProperties> - material = nullptr; - + boost::optional<MeshLib::PropertyVector<int> const&> material_ids; if (mesh.getProperties().existsPropertyVector<int>("MaterialIDs")) { INFO("The twophase flow is in heterogeneous porous media."); - const bool has_material_ids = true; - auto const& mat_ids = - mesh.getProperties().getPropertyVector<int>("MaterialIDs"); - material = MaterialLib::TwoPhaseFlowWithPP:: - CreateTwoPhaseFlowMaterialProperties(mat_config, has_material_ids, - *mat_ids); - TwoPhaseFlowWithPPProcessData process_data{ - specific_body_force, - has_gravity, - mass_lump, - std::move(material), - *curves.at("curve_PC_S"), - *curves.at("curve_S_Krel_wet"), - *curves.at("curve_S_Krel_nonwet")}; - return std::unique_ptr<Process>{new TwoPhaseFlowWithPPProcess{ - mesh, std::move(jacobian_assembler), parameters, integration_order, - std::move(process_variables), std::move(process_data), - std::move(secondary_variables), std::move(named_function_caller), - mat_config, curves}}; + material_ids = + *mesh.getProperties().getPropertyVector<int>("MaterialIDs"); } else { INFO("The twophase flow is in homogeneous porous media."); + } + std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties> + material = createTwoPhaseFlowWithPPMaterialProperties(mat_config, material_ids); - MeshLib::Properties dummy_property; - - auto const& dummy_property_vector = - dummy_property.createNewPropertyVector<int>( - "MaterialIDs", MeshLib::MeshItemType::Cell, 1); + TwoPhaseFlowWithPPProcessData process_data{ + specific_body_force, has_gravity, mass_lumping, temperature, std::move(material)}; - // Since dummy_property_vector is only visible in this function, - // the following constant, has_material_ids, is employed to indicate - // that material_ids does not exist. - const bool has_material_ids = false; - material = MaterialLib::TwoPhaseFlowWithPP:: - CreateTwoPhaseFlowMaterialProperties(mat_config, has_material_ids, - *dummy_property_vector); - TwoPhaseFlowWithPPProcessData process_data{ - specific_body_force, - has_gravity, - mass_lump, - std::move(material), - *curves.at("curve_PC_S"), - *curves.at("curve_S_Krel_wet"), - *curves.at("curve_S_Krel_nonwet")}; - return std::unique_ptr<Process>{new TwoPhaseFlowWithPPProcess{ - mesh, std::move(jacobian_assembler), parameters, integration_order, - std::move(process_variables), std::move(process_data), - std::move(secondary_variables), std::move(named_function_caller), - mat_config, curves}}; - } + return std::unique_ptr<Process>{new TwoPhaseFlowWithPPProcess{ + mesh, std::move(jacobian_assembler), parameters, integration_order, + std::move(process_variables), std::move(process_data), + std::move(secondary_variables), std::move(named_function_caller), + mat_config, curves}}; } } // end of namespace diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h index 9f46501de39..d6a43e27c3b 100644 --- a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h +++ b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h @@ -16,7 +16,7 @@ namespace ProcessLib { namespace TwoPhaseFlowWithPP { -std::unique_ptr<Process> CreateTwoPhaseFlowWithPPProcess( +std::unique_ptr<Process> createTwoPhaseFlowWithPPProcess( MeshLib::Mesh& mesh, std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler, std::vector<ProcessVariable> const& variables, diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h index 939b6644688..ee734e2f5e2 100644 --- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h +++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h @@ -7,6 +7,27 @@ * */ +/** +* common nomenclature +* --------------primary variable---------------------- +* pn_int_pt pressure for nonwetting phase at each integration point +* pc_int_pt capillary pressure at each integration point +* --------------secondary variable-------------------- +* temperature capillary pressure +* Sw wetting phase saturation +* dSw_dpc derivative of wetting phase saturation with respect +* to capillary pressure +* rho_nonwet density of nonwetting phase +* drhononwet_dpn derivative of nonwetting phase density with respect +*to nonwetting phase pressure +* rho_wet density of wetting phase +* k_rel_nonwet relative permeability of nonwetting phase +* mu_nonwet viscosity of nonwetting phase +* lambda_nonwet mobility of nonwetting phase +* k_rel_wet relative permeability of wetting phase +* mu_wet viscosity of wetting phase +* lambda_wet mobility of wetting phase +*/ #pragma once #include "TwoPhaseFlowWithPPLocalAssembler.h" @@ -39,18 +60,12 @@ void TwoPhaseFlowWithPPLocalAssembler< auto local_b = MathLib::createZeroedVector<LocalVectorType>( local_b_data, local_matrix_size); - NodalMatrixType mass_operator; - mass_operator.setZero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS); - auto Mgp = local_M.template block<nonwet_pressure_size, nonwet_pressure_size>( nonwet_pressure_matrix_index, nonwet_pressure_matrix_index); auto Mgpc = local_M.template block<nonwet_pressure_size, cap_pressure_size>( nonwet_pressure_matrix_index, cap_pressure_matrix_index); - auto Mlp = local_M.template block<cap_pressure_size, nonwet_pressure_size>( - cap_pressure_matrix_index, nonwet_pressure_matrix_index); - auto Mlpc = local_M.template block<cap_pressure_size, cap_pressure_size>( cap_pressure_matrix_index, cap_pressure_matrix_index); @@ -61,9 +76,6 @@ void TwoPhaseFlowWithPPLocalAssembler< local_K.template block<nonwet_pressure_size, nonwet_pressure_size>( nonwet_pressure_matrix_index, nonwet_pressure_matrix_index); - auto Kgpc = local_K.template block<nonwet_pressure_size, cap_pressure_size>( - nonwet_pressure_matrix_index, cap_pressure_matrix_index); - auto Klp = local_K.template block<cap_pressure_size, nonwet_pressure_size>( cap_pressure_matrix_index, nonwet_pressure_matrix_index); @@ -81,105 +93,97 @@ void TwoPhaseFlowWithPPLocalAssembler< SpatialPosition pos; pos.setElementID(_element.getID()); - _process_data._material->setMaterialID(pos); - - const Eigen::MatrixXd& perm = _process_data._material->getPermeability( - t, pos, _element.getDimension()); - assert(perm.rows() == GlobalDim || perm.rows() == 1); - GlobalDimMatrixType permeability = - GlobalDimMatrixType::Zero(GlobalDim, GlobalDim); - if (perm.rows() == GlobalDim) + const int material_id = + _process_data.material->getMaterialID(pos.getElementID().get()); + + const Eigen::MatrixXd& perm = _process_data.material->getPermeability( + material_id, t, pos, _element.getDimension()); + assert(perm.rows() == _element.getDimension() || perm.rows() == 1); + GlobalDimMatrixType permeability = GlobalDimMatrixType::Zero( + _element.getDimension(), _element.getDimension()); + if (perm.rows() == _element.getDimension()) permeability = perm; else if (perm.rows() == 1) permeability.diagonal().setConstant(perm(0, 0)); - MathLib::PiecewiseLinearInterpolation const& interpolated_Pc = - _process_data._interpolated_Pc; - MathLib::PiecewiseLinearInterpolation const& interpolated_Kr_wet = - _process_data._interpolated_Kr_wet; - MathLib::PiecewiseLinearInterpolation const& interpolated_Kr_nonwet = - _process_data._interpolated_Kr_nonwet; - - // Note: currently only isothermal case is considered, so the temperature is - // assumed to be const - // the variation of temperature will be taken into account in future + for (unsigned ip = 0; ip < n_integration_points; ip++) { - auto const& sm = _shape_matrices[ip]; - double pc_int_pt = 0.; - double pg_int_pt = 0.; - NumLib::shapeFunctionInterpolate(local_x, sm.N, pg_int_pt, pc_int_pt); - - _pressure_wetting[ip] = pg_int_pt - pc_int_pt; + double pn_int_pt = 0.; + NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N, pn_int_pt, + pc_int_pt); - auto const& wp = _integration_method.getWeightedPoint(ip); - const double integration_factor = - sm.integralMeasure * sm.detJ * wp.getWeight(); + _pressure_wet[ip] = pn_int_pt - pc_int_pt; - double const rho_gas = - _process_data._material->getGasDensity(pg_int_pt, _temperature); - double const rho_w = _process_data._material->getLiquidDensity( - _pressure_wetting[ip], _temperature); + const double temperature = _process_data.temperature(t, pos)[0]; + double const rho_nonwet = + _process_data.material->getGasDensity(pn_int_pt, temperature); + double const rho_wet = _process_data.material->getLiquidDensity( + _pressure_wet[ip], temperature); - double const Sw = - (pc_int_pt < 0) ? 1.0 : interpolated_Pc.getValue(pc_int_pt); + double const Sw = _process_data.material->getSaturation( + material_id, t, pos, pn_int_pt, temperature, pc_int_pt); _saturation[ip] = Sw; - double dSwdPc = interpolated_Pc.getDerivative(pc_int_pt); - if (pc_int_pt > interpolated_Pc.getSupportMax()) - dSwdPc = - interpolated_Pc.getDerivative(interpolated_Pc.getSupportMax()); - else if (pc_int_pt < interpolated_Pc.getSupportMin()) - dSwdPc = - interpolated_Pc.getDerivative(interpolated_Pc.getSupportMin()); - double const porosity = _process_data._material->getPorosity( - t, pos, pg_int_pt, _temperature, 0); + double dSw_dpc = _process_data.material->getSaturationDerivative( + material_id, t, pos, pn_int_pt, temperature, Sw); + + double const porosity = _process_data.material->getPorosity( + material_id, t, pos, pn_int_pt, temperature, 0); // Assemble M matrix // nonwetting - double const drhogas_dpg = _process_data._material->getDerivGasDensity( - pg_int_pt, _temperature); + double const drhononwet_dpn = + _process_data.material->getGasDensityDerivative(pn_int_pt, + temperature); - mass_operator.noalias() = sm.N.transpose() * sm.N * integration_factor; + Mgp.noalias() += + porosity * (1 - Sw) * drhononwet_dpn * _ip_data[ip].massOperator; + Mgpc.noalias() += + -porosity * rho_nonwet * dSw_dpc * _ip_data[ip].massOperator; - Mgp.noalias() += porosity * (1 - Sw) * drhogas_dpg * mass_operator; - Mgpc.noalias() += -porosity * rho_gas * dSwdPc * mass_operator; - Mlp.noalias() += 0.0 * mass_operator; - Mlpc.noalias() += porosity * dSwdPc * rho_w * mass_operator; + Mlpc.noalias() += + porosity * dSw_dpc * rho_wet * _ip_data[ip].massOperator; - // Assemble M matrix // nonwet - double const k_rel_G = interpolated_Kr_nonwet.getValue(Sw); - double const mu_gas = - _process_data._material->getGasViscosity(pg_int_pt, _temperature); - double const lambda_G = k_rel_G / mu_gas; + double const k_rel_nonwet = + _process_data.material->getNonwetRelativePermeability( + t, pos, pn_int_pt, temperature, Sw); + double const mu_nonwet = + _process_data.material->getGasViscosity(pn_int_pt, temperature); + double const lambda_nonwet = k_rel_nonwet / mu_nonwet; // wet - double const k_rel_L = interpolated_Kr_wet.getValue(Sw); - double const mu_liquid = _process_data._material->getLiquidViscosity( - _pressure_wetting[ip], _temperature); - double const lambda_L = k_rel_L / mu_liquid; + double const k_rel_wet = + _process_data.material->getWetRelativePermeability( + t, pos, _pressure_wet[ip], temperature, Sw); + double const mu_wet = _process_data.material->getLiquidViscosity( + _pressure_wet[ip], temperature); + double const lambda_wet = k_rel_wet / mu_wet; - laplace_operator.noalias() = - sm.dNdx.transpose() * permeability * sm.dNdx * integration_factor; + laplace_operator.noalias() = _ip_data[ip].dNdx.transpose() * + permeability * _ip_data[ip].dNdx * + _ip_data[ip].integration_weight; - Kgp.noalias() += rho_gas * lambda_G * laplace_operator; - Kgpc.noalias() += 0.0 * laplace_operator; - Klp.noalias() += rho_w * lambda_L * laplace_operator; - Klpc.noalias() += -rho_w * lambda_L * laplace_operator; + Kgp.noalias() += rho_nonwet * lambda_nonwet * laplace_operator; - if (_process_data._has_gravity) - { - auto const& b = _process_data._specific_body_force; - Bg.noalias() += rho_gas * rho_gas * lambda_G * sm.dNdx.transpose() * - permeability * b * integration_factor; - Bl.noalias() += rho_w * rho_w * lambda_L * sm.dNdx.transpose() * - permeability * b * integration_factor; + Klp.noalias() += rho_wet * lambda_wet * laplace_operator; + Klpc.noalias() += -rho_wet * lambda_wet * laplace_operator; + if (_process_data.has_gravity) + { + auto const& b = _process_data.specific_body_force; + + NodalVectorType gravity_operator = _ip_data[ip].dNdx.transpose() * + permeability * b * + _ip_data[ip].integration_weight; + Bg.noalias() += + rho_nonwet * rho_nonwet * lambda_nonwet * gravity_operator; + Bl.noalias() += rho_wet * rho_wet * lambda_wet * gravity_operator; } // end of has gravity - } // end of GP - if (_process_data._has_mass_lumping) + } + if (_process_data.has_mass_lumping) { for (unsigned row = 0; row < Mgpc.cols(); row++) { @@ -193,8 +197,6 @@ void TwoPhaseFlowWithPPLocalAssembler< Mgp(row, column) = 0.0; Mlpc(row, row) += Mlpc(row, column); Mlpc(row, column) = 0.0; - Mlp(row, row) += Mlp(row, column); - Mlp(row, column) = 0.0; } } } diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h index f9d65553af8..2daf43b4cab 100644 --- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h +++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h @@ -20,13 +20,35 @@ #include "ProcessLib/Parameter/Parameter.h" #include "ProcessLib/Utils/InitShapeMatrices.h" -#include "MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h" +#include "TwoPhaseFlowWithPPMaterialProperties.h" #include "TwoPhaseFlowWithPPProcessData.h" namespace ProcessLib { namespace TwoPhaseFlowWithPP { +template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType, + typename NodalMatrixType> +struct IntegrationPointData final +{ + explicit IntegrationPointData( + NodalRowVectorType const& N_, GlobalDimNodalMatrixType const& dNdx_, + TwoPhaseFlowWithPPMaterialProperties& material_property_, + double const& integration_weight_, NodalMatrixType const massOperator_) + : N(N_), + dNdx(dNdx_), + mat_property(material_property_), + integration_weight(integration_weight_), + massOperator(massOperator_) + + { + } + NodalRowVectorType const N; + GlobalDimNodalMatrixType const dNdx; + TwoPhaseFlowWithPPMaterialProperties const& mat_property; + const double integration_weight; + NodalMatrixType const massOperator; +}; const unsigned NUM_NODAL_DOF = 2; class TwoPhaseFlowWithPPLocalAssemblerInterface @@ -37,7 +59,7 @@ public: virtual std::vector<double> const& getIntPtSaturation( std::vector<double>& /*cache*/) const = 0; - virtual std::vector<double> const& getIntPtWettingPressure( + virtual std::vector<double> const& getIntPtWetPressure( std::vector<double>& /*cache*/) const = 0; }; @@ -51,7 +73,10 @@ class TwoPhaseFlowWithPPLocalAssembler using LocalAssemblerTraits = ProcessLib::LocalAssemblerTraits< ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>; + using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType; + using GlobalDimNodalMatrixType = + typename ShapeMatricesType::GlobalDimNodalMatrixType; using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType; using NodalVectorType = typename ShapeMatricesType::NodalVectorType; using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType; @@ -68,15 +93,29 @@ public: TwoPhaseFlowWithPPProcessData const& process_data) : _element(element), _integration_method(integration_order), - _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType, - IntegrationMethod, GlobalDim>( - element, is_axially_symmetric, _integration_method)), _process_data(process_data), _saturation( std::vector<double>(_integration_method.getNumberOfPoints())), - _pressure_wetting( + _pressure_wet( std::vector<double>(_integration_method.getNumberOfPoints())) { + unsigned const n_integration_points = + _integration_method.getNumberOfPoints(); + _ip_data.reserve(n_integration_points); + auto const shape_matrices = + initShapeMatrices<ShapeFunction, ShapeMatricesType, + IntegrationMethod, GlobalDim>( + element, is_axially_symmetric, _integration_method); + for (unsigned ip = 0; ip < n_integration_points; ip++) + { + auto const& sm = shape_matrices[ip]; + _ip_data.emplace_back( + sm.N, sm.dNdx, *_process_data.material, + sm.integralMeasure * sm.detJ * + _integration_method.getWeightedPoint(ip).getWeight(), + sm.N.transpose() * sm.N * sm.integralMeasure * sm.detJ * + _integration_method.getWeightedPoint(ip).getWeight()); + } } void assemble(double const t, std::vector<double> const& local_x, @@ -87,7 +126,7 @@ public: Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix( const unsigned integration_point) const override { - auto const& N = _shape_matrices[integration_point].N; + auto const& N = _ip_data[integration_point].N; // assumes N is stored contiguously in memory return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size()); @@ -100,28 +139,32 @@ public: return _saturation; } - std::vector<double> const& getIntPtWettingPressure( + std::vector<double> const& getIntPtWetPressure( std::vector<double>& /*cache*/) const override { - assert(_pressure_wetting.size() > 0); - return _pressure_wetting; + assert(_pressure_wet.size() > 0); + return _pressure_wet; } private: MeshLib::Element const& _element; IntegrationMethod const _integration_method; - std::vector<ShapeMatrices, Eigen::aligned_allocator<ShapeMatrices>> - _shape_matrices; TwoPhaseFlowWithPPProcessData const& _process_data; - - // Note: currently only isothermal case is considered, so the temperature is - // assumed to be const - // the variation of temperature will be taken into account in future - double _temperature = 293.15; + std::vector< + IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType, + NodalMatrixType>, + Eigen::aligned_allocator<IntegrationPointData< + NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType>>> + _ip_data; + + // output vector for wetting phase saturation with + // respect to each integration point std::vector<double> _saturation; - std::vector<double> _pressure_wetting; + // output vector for wetting phase pressure with respect + // to each integration point + std::vector<double> _pressure_wet; static const int nonwet_pressure_coeff_index = 0; static const int cap_pressure_coeff_index = 1; diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.cpp new file mode 100644 index 00000000000..c2ed90e5a79 --- /dev/null +++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.cpp @@ -0,0 +1,169 @@ +/** +* \copyright +* Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org) +* Distributed under a Modified BSD License. +* See accompanying file LICENSE.txt or +* http://www.opengeosys.org/project/license +* +*/ + +#include "TwoPhaseFlowWithPPMaterialProperties.h" +#include <logog/include/logog.hpp> +#include "BaseLib/reorderVector.h" +#include "MaterialLib/Fluid/FluidProperty.h" +#include "MaterialLib/PorousMedium/Porosity/Porosity.h" +#include "MaterialLib/PorousMedium/Storage/Storage.h" +#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h" +#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h" +#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h" +#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h" +#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h" +#include "MeshLib/Mesh.h" +#include "MeshLib/PropertyVector.h" +#include "ProcessLib/Parameter/Parameter.h" +#include "ProcessLib/Parameter/SpatialPosition.h" +namespace ProcessLib +{ +namespace TwoPhaseFlowWithPP +{ +TwoPhaseFlowWithPPMaterialProperties::TwoPhaseFlowWithPPMaterialProperties( + boost::optional<MeshLib::PropertyVector<int> const&> const material_ids, + std::unique_ptr<MaterialLib::Fluid::FluidProperty> + liquid_density, + std::unique_ptr<MaterialLib::Fluid::FluidProperty> + liquid_viscosity, + std::unique_ptr<MaterialLib::Fluid::FluidProperty> + gas_density, + std::unique_ptr<MaterialLib::Fluid::FluidProperty> + gas_viscosity, + std::vector<Eigen::MatrixXd> + intrinsic_permeability_models, + std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>&& + porosity_models, + std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>&& + storage_models, + std::vector<std::unique_ptr< + MaterialLib::PorousMedium::CapillaryPressureSaturation>>&& + capillary_pressure_models, + std::vector< + std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&& + relative_permeability_models) + : _liquid_density(std::move(liquid_density)), + _liquid_viscosity(std::move(liquid_viscosity)), + _gas_density(std::move(gas_density)), + _gas_viscosity(std::move(gas_viscosity)), + _material_ids(material_ids), + _intrinsic_permeability_models(intrinsic_permeability_models), + _porosity_models(std::move(porosity_models)), + _storage_models(std::move(storage_models)), + _capillary_pressure_models(std::move(capillary_pressure_models)), + _relative_permeability_models(std::move(relative_permeability_models)) +{ + DBUG("Create material properties for Two-Phase flow with PP model."); +} + +int TwoPhaseFlowWithPPMaterialProperties::getMaterialID( + const std::size_t element_id) +{ + if (!_material_ids) + { + return 0; + } + + assert(element_id < _material_ids->size()); + return (*_material_ids)[element_id]; +} + +double TwoPhaseFlowWithPPMaterialProperties::getLiquidDensity( + const double p, const double T) const +{ + ArrayType vars; + vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T; + vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p; + return _liquid_density->getValue(vars); +} + +double TwoPhaseFlowWithPPMaterialProperties::getGasDensity(const double p, + const double T) const +{ + ArrayType vars; + vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T; + vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p; + return _gas_density->getValue(vars); +} + +double TwoPhaseFlowWithPPMaterialProperties::getGasDensityDerivative( + const double p, const double T) const +{ + ArrayType vars; + vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T; + vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p; + + return _gas_density->getdValue(vars, + MaterialLib::Fluid::PropertyVariableType::p); +} +double TwoPhaseFlowWithPPMaterialProperties::getLiquidViscosity( + const double p, const double T) const +{ + ArrayType vars; + vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T; + vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p; + return _liquid_viscosity->getValue(vars); +} + +double TwoPhaseFlowWithPPMaterialProperties::getGasViscosity( + const double p, const double T) const +{ + ArrayType vars; + vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T; + vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] = p; + return _gas_viscosity->getValue(vars); +} + +Eigen::MatrixXd const& TwoPhaseFlowWithPPMaterialProperties::getPermeability( + const int material_id, const double /*t*/, + const ProcessLib::SpatialPosition& /*pos*/, const int /*dim*/) const +{ + return _intrinsic_permeability_models[material_id]; +} + +double TwoPhaseFlowWithPPMaterialProperties::getPorosity( + const int material_id, const double /*t*/, + const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/, + const double T, const double porosity_variable) const +{ + return _porosity_models[material_id]->getValue(porosity_variable, T); +} + +double TwoPhaseFlowWithPPMaterialProperties::getNonwetRelativePermeability( + const double /*t*/, const ProcessLib::SpatialPosition& /*pos*/, + const double /*p*/, const double /*T*/, const double saturation) const +{ + return _relative_permeability_models[0]->getValue(saturation); +} + +double TwoPhaseFlowWithPPMaterialProperties::getWetRelativePermeability( + const double /*t*/, const ProcessLib::SpatialPosition& /*pos*/, + const double /*p*/, const double /*T*/, const double saturation) const +{ + return _relative_permeability_models[1]->getValue(saturation); +} + +double TwoPhaseFlowWithPPMaterialProperties::getSaturation( + const int material_id, const double /*t*/, + const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/, + const double /*T*/, const double pc) const +{ + return _capillary_pressure_models[material_id]->getSaturation(pc); +} +double TwoPhaseFlowWithPPMaterialProperties::getSaturationDerivative( + const int material_id, const double /*t*/, + const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/, + const double /*T*/, const double saturation) const +{ + const double dpcdsw = + _capillary_pressure_models[material_id]->getdPcdS(saturation); + return 1 / dpcdsw; +} +} // end of namespace +} // end of namespace diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.h new file mode 100644 index 00000000000..fd3b0657d01 --- /dev/null +++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.h @@ -0,0 +1,129 @@ +/** +* \copyright +* Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org) +* Distributed under a Modified BSD License. +* See accompanying file LICENSE.txt or +* http://www.opengeosys.org/project/license +* +*/ + +#pragma once + +#include <memory> +#include <vector> +#include "MaterialLib/Fluid/FluidPropertyHeaders.h" +#include "MaterialLib/PhysicalConstant.h" +#include "MaterialLib/PorousMedium/Porosity/Porosity.h" +#include "MaterialLib/PorousMedium/PorousPropertyHeaders.h" +#include "MaterialLib/PorousMedium/Storage/Storage.h" +#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h" +#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h" +#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h" +#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h" + +namespace MeshLib +{ +template <typename PROP_VAL_TYPE> +class PropertyVector; +} + +namespace ProcessLib +{ +class SpatialPosition; +namespace TwoPhaseFlowWithPP +{ + /** This class has a collection of material properties for two-phase flow with PP model + * and it makes description of the material properties for two-phase condition, + * i.e. the gas/liquid density and viscosity models, respectively, + * the relative permeability models with respect to two phases, + * the capillary pressure-saturation relationships. + * It generally provides the computation of the PDE coefficients for two-phase flow. + */ + +class TwoPhaseFlowWithPPMaterialProperties +{ +public: + using ArrayType = MaterialLib::Fluid::FluidProperty::ArrayType; + + TwoPhaseFlowWithPPMaterialProperties( + boost::optional<MeshLib::PropertyVector<int> const&> const material_ids, + std::unique_ptr<MaterialLib::Fluid::FluidProperty> + liquid_density, + std::unique_ptr<MaterialLib::Fluid::FluidProperty> + liquid_viscosity, + std::unique_ptr<MaterialLib::Fluid::FluidProperty> + gas_density, + std::unique_ptr<MaterialLib::Fluid::FluidProperty> + gas_viscosity, + std::vector<Eigen::MatrixXd> + intrinsic_permeability_models, + std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>&& + porosity_models, + std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>&& + storage_models, + std::vector<std::unique_ptr< + MaterialLib::PorousMedium::CapillaryPressureSaturation>>&& + capillary_pressure_models, + std::vector< + std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&& + relative_permeability_models); + + int getMaterialID(const std::size_t element_id); + + Eigen::MatrixXd const& getPermeability( + const int material_id, + const double t, + const ProcessLib::SpatialPosition& pos, + const int dim) const; + + double getPorosity(const int material_id, const double t, + const ProcessLib::SpatialPosition& pos, const double p, + const double T, const double porosity_variable) const; + + double getNonwetRelativePermeability(const double t, + const ProcessLib::SpatialPosition& pos, + const double p, const double T, + const double saturation) const; + double getWetRelativePermeability(const double t, + const ProcessLib::SpatialPosition& pos, + const double p, const double T, + const double saturation) const; + double getSaturation(const int material_id, const double t, + const ProcessLib::SpatialPosition& pos, const double p, + const double T, const double pc) const; + double getSaturationDerivative(const int material_id, const double t, + const ProcessLib::SpatialPosition& pos, + const double p, const double T, + const double saturation) const; + double getLiquidDensity(const double p, const double T) const; + double getGasDensity(const double p, const double T) const; + double getGasViscosity(const double p, const double T) const; + double getLiquidViscosity(const double p, const double T) const; + double getGasDensityDerivative(double const p, double const T) const; + +protected: + std::unique_ptr<MaterialLib::Fluid::FluidProperty> _liquid_density; + std::unique_ptr<MaterialLib::Fluid::FluidProperty> _liquid_viscosity; + std::unique_ptr<MaterialLib::Fluid::FluidProperty> _gas_density; + std::unique_ptr<MaterialLib::Fluid::FluidProperty> _gas_viscosity; + + /** Use two phase models for different material zones. + * Material IDs must be given as mesh element properties. + */ + boost::optional<MeshLib::PropertyVector<int> const&> const _material_ids; + + std::vector<Eigen::MatrixXd> _intrinsic_permeability_models; + std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>> + _porosity_models; + std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>> + _storage_models; + std::vector< + std::unique_ptr<MaterialLib::PorousMedium::CapillaryPressureSaturation>> + _capillary_pressure_models; + std::vector< + std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>> + _relative_permeability_models; +}; + +} // end of namespace +} // end of namespace diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp index 2cbb77f5b51..6c855e3ba4c 100644 --- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp +++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp @@ -35,7 +35,7 @@ TwoPhaseFlowWithPPProcess::TwoPhaseFlowWithPPProcess( BaseLib::ConfigTree const& /*config*/, std::map<std::string, std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const& - /*curves*/) + /*curves*/) : Process(mesh, std::move(jacobian_assembler), parameters, integration_order, std::move(process_variables), std::move(secondary_variables), std::move(named_function_caller)), @@ -62,10 +62,10 @@ void TwoPhaseFlowWithPPProcess::initializeConcreteProcess( &TwoPhaseFlowWithPPLocalAssemblerInterface::getIntPtSaturation)); _secondary_variables.addSecondaryVariable( - "pressure_wetting", 1, + "pressure_wet", 1, makeExtrapolator(getExtrapolator(), _local_assemblers, &TwoPhaseFlowWithPPLocalAssemblerInterface:: - getIntPtWettingPressure)); + getIntPtWetPressure)); } void TwoPhaseFlowWithPPProcess::assembleConcreteProcess(const double t, diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h index 42854353fe4..376d91d356f 100644 --- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h +++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h @@ -9,7 +9,6 @@ #pragma once -#include "MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h" #include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h" #include "NumLib/DOF/LocalToGlobalIndexMap.h" #include "ProcessLib/Process.h" diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h index 6da20a21694..bf722b2473d 100644 --- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h +++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h @@ -8,7 +8,7 @@ */ #pragma once - +#include "TwoPhaseFlowWithPPMaterialProperties.h" namespace MeshLib { class Element; @@ -27,29 +27,23 @@ struct TwoPhaseFlowWithPPProcessData Eigen::VectorXd const specific_body_force_, bool const has_gravity_, bool const has_mass_lumping_, - std::unique_ptr<MaterialLib::TwoPhaseFlowWithPP:: - TwoPhaseFlowWithPPMaterialProperties>&& material_, - MathLib::PiecewiseLinearInterpolation const& interpolated_Pc_, - MathLib::PiecewiseLinearInterpolation const& interpolated_Kr_wet_, - MathLib::PiecewiseLinearInterpolation const& interpolated_Kr_nonwet_) - : _specific_body_force(specific_body_force_), - _has_gravity(has_gravity_), - _has_mass_lumping(has_mass_lumping_), - _material(std::move(material_)), - _interpolated_Pc(interpolated_Pc_), - _interpolated_Kr_wet(interpolated_Kr_wet_), - _interpolated_Kr_nonwet(interpolated_Kr_nonwet_) + Parameter<double> const& temperature_, + std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>&& material_) + : specific_body_force(specific_body_force_), + has_gravity(has_gravity_), + has_mass_lumping(has_mass_lumping_), + temperature(temperature_), + material(std::move(material_)) + { } TwoPhaseFlowWithPPProcessData(TwoPhaseFlowWithPPProcessData&& other) - : _specific_body_force(other._specific_body_force), - _has_gravity(other._has_gravity), - _has_mass_lumping(other._has_mass_lumping), - _material(std::move(other._material)), - _interpolated_Pc(other._interpolated_Pc), - _interpolated_Kr_wet(other._interpolated_Kr_wet), - _interpolated_Kr_nonwet(other._interpolated_Kr_nonwet) + : specific_body_force(other.specific_body_force), + has_gravity(other.has_gravity), + has_mass_lumping(other.has_mass_lumping), + temperature(other.temperature), + material(std::move(other.material)) { } @@ -66,18 +60,14 @@ struct TwoPhaseFlowWithPPProcessData //! Specific body forces applied to solid and fluid. //! It is usually used to apply gravitational forces. //! A vector of displacement dimension's length. - Eigen::VectorXd const _specific_body_force; + Eigen::VectorXd const specific_body_force; - bool const _has_gravity; + bool const has_gravity; //! Enables lumping of the mass matrix. - bool const _has_mass_lumping; - std::unique_ptr< - MaterialLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPMaterialProperties> - _material; - MathLib::PiecewiseLinearInterpolation const& _interpolated_Pc; - MathLib::PiecewiseLinearInterpolation const& _interpolated_Kr_wet; - MathLib::PiecewiseLinearInterpolation const& _interpolated_Kr_nonwet; + bool const has_mass_lumping; + Parameter<double> const& temperature; + std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties> material; }; } // namespace TwoPhaseFlowWithPP -- GitLab