From ba7ae6bd204202df982ce20747ed578860da3950 Mon Sep 17 00:00:00 2001
From: tkurz <theresa.kurz@iws.uni-stuttgart.de>
Date: Tue, 18 Jun 2019 22:39:14 +0200
Subject: [PATCH] [test][2pncminni] add a test for evaporation with salt
precipitation for the 2pncmin nonisothermal model
---
.../2pncmin/implicit/CMakeLists.txt | 38 +-
.../implicit/isothermal/CMakeLists.txt | 34 ++
.../2pncmin/implicit/{ => isothermal}/main.cc | 0
.../implicit/{ => isothermal}/params.input | 0
.../implicit/{ => isothermal}/problem.hh | 0
.../{ => isothermal}/spatialparams.hh | 0
.../implicit/nonisothermal/CMakeLists.txt | 12 +
.../2pncmin/implicit/nonisothermal/main.cc | 247 ++++++++
.../implicit/nonisothermal/params.input | 52 ++
.../2pncmin/implicit/nonisothermal/problem.hh | 568 ++++++++++++++++++
.../implicit/nonisothermal/spatialparams.hh | 181 ++++++
.../test_2pncminni_salinization-reference.vtu | 139 +++++
12 files changed, 1235 insertions(+), 36 deletions(-)
create mode 100644 test/porousmediumflow/2pncmin/implicit/isothermal/CMakeLists.txt
rename test/porousmediumflow/2pncmin/implicit/{ => isothermal}/main.cc (100%)
rename test/porousmediumflow/2pncmin/implicit/{ => isothermal}/params.input (100%)
rename test/porousmediumflow/2pncmin/implicit/{ => isothermal}/problem.hh (100%)
rename test/porousmediumflow/2pncmin/implicit/{ => isothermal}/spatialparams.hh (100%)
create mode 100644 test/porousmediumflow/2pncmin/implicit/nonisothermal/CMakeLists.txt
create mode 100644 test/porousmediumflow/2pncmin/implicit/nonisothermal/main.cc
create mode 100644 test/porousmediumflow/2pncmin/implicit/nonisothermal/params.input
create mode 100644 test/porousmediumflow/2pncmin/implicit/nonisothermal/problem.hh
create mode 100644 test/porousmediumflow/2pncmin/implicit/nonisothermal/spatialparams.hh
create mode 100644 test/references/test_2pncminni_salinization-reference.vtu
diff --git a/test/porousmediumflow/2pncmin/implicit/CMakeLists.txt b/test/porousmediumflow/2pncmin/implicit/CMakeLists.txt
index eb412f736f..bc7eb9a410 100644
--- a/test/porousmediumflow/2pncmin/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/2pncmin/implicit/CMakeLists.txt
@@ -1,36 +1,2 @@
-dune_symlink_to_source_files(FILES params.input)
-
-# isothermal tests
-dune_add_test(NAME test_2pncmin_dissolution_box
- LABELS porousmediumflow 2pncmin
- SOURCES main.cc
- COMPILE_DEFINITIONS TYPETAG=DissolutionBox
- TIMEOUT 1500
- COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
- CMD_ARGS --script fuzzy
- --files ${CMAKE_SOURCE_DIR}/test/references/test_2pncmin_dissolution_box-reference.vtu
- ${CMAKE_CURRENT_BINARY_DIR}/test_2pncmin_dissolution_box-00044.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_2pncmin_dissolution_box -ParameterFile params.input -Problem.Name test_2pncmin_dissolution_box")
-
-dune_add_test(NAME test_2pncmin_dissolution_box_restart
- LABELS porousmediumflow 2pncmin
- TARGET test_2pncmin_dissolution_box
- COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
- CMD_ARGS --script fuzzy
- --files ${CMAKE_SOURCE_DIR}/test/references/test_2pncmin_dissolution_box-reference.vtu
- ${CMAKE_CURRENT_BINARY_DIR}/test_2pncmin_dissolution_box_restart-00005.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_2pncmin_dissolution_box params.input -Problem.Name test_2pncmin_dissolution_box_restart -TimeLoop.DtInitial 50000 -Restart.Time 756290 -Restart.File test_2pncmin_dissolution_box-00039.vtu")
-
-# the restart test has to run after the test that produces the corresponding vtu file
-set_tests_properties(test_2pncmin_dissolution_box_restart PROPERTIES DEPENDS test_2pncmin_dissolution_box)
-
-dune_add_test(NAME test_2pncmin_dissolution_tpfa
- LABELS porousmediumflow 2pncmin
- SOURCES main.cc
- COMPILE_DEFINITIONS TYPETAG=DissolutionCCTpfa
- TIMEOUT 1500
- COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
- CMD_ARGS --script fuzzy
- --files ${CMAKE_SOURCE_DIR}/test/references/test_2pncmin_dissolution_tpfa-reference.vtu
- ${CMAKE_CURRENT_BINARY_DIR}/test_2pncmin_dissolution_tpfa-00043.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_2pncmin_dissolution_tpfa -ParameterFile params.input -Problem.Name test_2pncmin_dissolution_tpfa")
+add_subdirectory(isothermal)
+add_subdirectory(nonisothermal)
diff --git a/test/porousmediumflow/2pncmin/implicit/isothermal/CMakeLists.txt b/test/porousmediumflow/2pncmin/implicit/isothermal/CMakeLists.txt
new file mode 100644
index 0000000000..4b2d0333ce
--- /dev/null
+++ b/test/porousmediumflow/2pncmin/implicit/isothermal/CMakeLists.txt
@@ -0,0 +1,34 @@
+dune_symlink_to_source_files(FILES params.input)
+
+# isothermal tests
+dune_add_test(NAME test_2pncmin_dissolution_box
+ LABELS porousmediumflow 2pncmin
+ SOURCES main.cc
+ COMPILE_DEFINITIONS TYPETAG=DissolutionBox
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_2pncmin_dissolution_box-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_2pncmin_dissolution_box-00044.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_2pncmin_dissolution_box -ParameterFile params.input -Problem.Name test_2pncmin_dissolution_box")
+
+dune_add_test(NAME test_2pncmin_dissolution_box_restart
+ LABELS porousmediumflow 2pncmin
+ TARGET test_2pncmin_dissolution_box
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_2pncmin_dissolution_box-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_2pncmin_dissolution_box_restart-00005.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_2pncmin_dissolution_box params.input -Problem.Name test_2pncmin_dissolution_box_restart -TimeLoop.DtInitial 50000 -Restart.Time 756290 -Restart.File test_2pncmin_dissolution_box-00039.vtu")
+
+# the restart test has to run after the test that produces the corresponding vtu file
+set_tests_properties(test_2pncmin_dissolution_box_restart PROPERTIES DEPENDS test_2pncmin_dissolution_box)
+
+dune_add_test(NAME test_2pncmin_dissolution_tpfa
+ LABELS porousmediumflow 2pncmin
+ SOURCES main.cc
+ COMPILE_DEFINITIONS TYPETAG=DissolutionCCTpfa
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_2pncmin_dissolution_tpfa-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_2pncmin_dissolution_tpfa-00043.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_2pncmin_dissolution_tpfa -ParameterFile params.input -Problem.Name test_2pncmin_dissolution_tpfa")
diff --git a/test/porousmediumflow/2pncmin/implicit/main.cc b/test/porousmediumflow/2pncmin/implicit/isothermal/main.cc
similarity index 100%
rename from test/porousmediumflow/2pncmin/implicit/main.cc
rename to test/porousmediumflow/2pncmin/implicit/isothermal/main.cc
diff --git a/test/porousmediumflow/2pncmin/implicit/params.input b/test/porousmediumflow/2pncmin/implicit/isothermal/params.input
similarity index 100%
rename from test/porousmediumflow/2pncmin/implicit/params.input
rename to test/porousmediumflow/2pncmin/implicit/isothermal/params.input
diff --git a/test/porousmediumflow/2pncmin/implicit/problem.hh b/test/porousmediumflow/2pncmin/implicit/isothermal/problem.hh
similarity index 100%
rename from test/porousmediumflow/2pncmin/implicit/problem.hh
rename to test/porousmediumflow/2pncmin/implicit/isothermal/problem.hh
diff --git a/test/porousmediumflow/2pncmin/implicit/spatialparams.hh b/test/porousmediumflow/2pncmin/implicit/isothermal/spatialparams.hh
similarity index 100%
rename from test/porousmediumflow/2pncmin/implicit/spatialparams.hh
rename to test/porousmediumflow/2pncmin/implicit/isothermal/spatialparams.hh
diff --git a/test/porousmediumflow/2pncmin/implicit/nonisothermal/CMakeLists.txt b/test/porousmediumflow/2pncmin/implicit/nonisothermal/CMakeLists.txt
new file mode 100644
index 0000000000..082a0bbdb4
--- /dev/null
+++ b/test/porousmediumflow/2pncmin/implicit/nonisothermal/CMakeLists.txt
@@ -0,0 +1,12 @@
+dune_symlink_to_source_files(FILES params.input)
+
+# isothermal tests
+dune_add_test(NAME test_2pncminni_salinization
+ LABELS porousmediumflow 2pncmin
+ SOURCES main.cc
+ COMPILE_DEFINITIONS TYPETAG=SalinizationBox
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/test_2pncminni_salinization-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_2pncminni_salinization-00034.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_2pncminni_salinization -ParameterFile params.input -Problem.Name test_2pncminni_salinization")
diff --git a/test/porousmediumflow/2pncmin/implicit/nonisothermal/main.cc b/test/porousmediumflow/2pncmin/implicit/nonisothermal/main.cc
new file mode 100644
index 0000000000..d6ab3b32ed
--- /dev/null
+++ b/test/porousmediumflow/2pncmin/implicit/nonisothermal/main.cc
@@ -0,0 +1,247 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup TwoPNCMinTests
+ * \brief Test for the nonisothermal two-phase n-component finite volume model used to model e.g. salt precipitation.
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+#include <dune/istl/io.hh>
+
+#include "problem.hh"
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/method.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+#include <dumux/io/loadsolution.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+ "\t-ParameterFile Parameter file (Input file) \n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = Properties::TTag::TYPETAG;
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv, usage);
+
+ // try to create a grid (from the given grid file or the input file)
+ GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
+ gridManager.init();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // get some time loop parameters
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = getParam<Scalar>("Restart.Time", 0);
+
+ // the solution vector
+ using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+ SolutionVector x(fvGridGeometry->numDofs());
+ if (restartTime > 0)
+ {
+ using IOFields = GetPropType<TypeTag, Properties::IOFields>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using SolidSystem = GetPropType<TypeTag, Properties::SolidSystem>;
+ const auto fileName = getParam<std::string>("Restart.File");
+ const auto pvName = createPVNameFunction<IOFields, PrimaryVariables, ModelTraits, FluidSystem, SolidSystem>();
+ loadSolution(x, fileName, pvName, *fvGridGeometry);
+ }
+ else
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x);
+
+ // initialize the vtk output module
+ using IOFields = GetPropType<TypeTag, Properties::IOFields>;
+ VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
+ using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
+ vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
+ IOFields::initOutputModule(vtkWriter); // Add model specific output fields
+ //add specific output
+ vtkWriter.addField(problem->getPermeability(), "Permeability");
+ // update the output fields before write
+ problem->updateVtkOutput(x);
+ vtkWriter.write(restartTime);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = AMGBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+
+ // the non-linear solver
+ using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set time for problem for implicit Euler scheme
+ problem->setTime( timeLoop->time() + timeLoop->timeStepSize() );
+ problem->setTimeStepSize( timeLoop->timeStepSize() );
+
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(xOld);
+
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(x, *timeLoop);
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // update the output fields before write
+ problem->updateVtkOutput(x);
+
+ // write vtk output
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by the newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/test/porousmediumflow/2pncmin/implicit/nonisothermal/params.input b/test/porousmediumflow/2pncmin/implicit/nonisothermal/params.input
new file mode 100644
index 0000000000..67456e2cae
--- /dev/null
+++ b/test/porousmediumflow/2pncmin/implicit/nonisothermal/params.input
@@ -0,0 +1,52 @@
+# Parameter file for test case 2pncmin.
+# Everything behind a '#' is a comment.
+# Type "./test_2pncmin --help" for more information.
+
+[TimeLoop]
+TEnd = 80000 # duration of the simulation [s]
+DtInitial = 1 # initial time step size [s]
+MaxTimeStepSize = 5000 # maximum time step size
+
+[Grid]
+Positions0 = 0 0.02
+Positions1 = 0 0.015 0.02
+Cells0 = 2
+Cells1 = 2 3
+Grading0 = 1.0
+Grading1 = 1.0 -1.2
+
+[FluidSystem]
+NTemperature = 100 # [-] number of tabularization entries
+NPressure = 100 # [-] number of tabularization entries
+PressureLow = 1e4 # [Pa]low end for tabularization of fluid properties
+PressureHigh = 3e6 # [Pa]high end for tabularization of fluid properties
+TemperatureLow = 273.15 # [K]low end for tabularization of fluid properties
+TemperatureHigh = 400.00 # [K]high end for tabularization of fluid properties
+
+[Problem]
+Name = test_2pncminni # [-] name for output files
+Temperature = 299 # [K] initial temperature
+InitialPressure = 1.00e5 # [Pa] Initial reservoir pressure
+InitialGasSaturation = 0.5 # [-] initial gas saturation
+InitialSalinity = 2.00e-2 # [-]=[kgNaCl/kgSolution] Initial salt mass fraction
+
+[FreeFlow]
+RefMoleFracH2O = 6.0e-3 # [-] mole fraction of water vapor in air above the container
+BoundaryLayerThickness = 1e-2 # [m] thickness of boundary layer over which the water vapor diffuses
+MassTransferCoefficient = 1
+RefTemperature = 299 # [K] temperature in the air above container
+
+[SpatialParams]
+SolubilityLimit = 0.26 # [-] solubility limit of NaCl in brine
+referencePorosity = 0.4 # [-] initial porosity
+referencePermeability = 1e-13 # [m²] initial permeability
+VGAlpha = 0.00015 # Van Genuchten parameter
+VGn = 14.0 # Van Genuchten parameter
+IrreducibleLiqSat = 0.25 # [-] irreducible liquid saturation
+IrreducibleGasSat = 0.025 # [-] irreducible gas saturation
+
+[Vtk]
+AddVelocity = 1 # Add extra information
+
+[Output]
+PlotFluidMatrixInteractions = false
diff --git a/test/porousmediumflow/2pncmin/implicit/nonisothermal/problem.hh b/test/porousmediumflow/2pncmin/implicit/nonisothermal/problem.hh
new file mode 100644
index 0000000000..8756092803
--- /dev/null
+++ b/test/porousmediumflow/2pncmin/implicit/nonisothermal/problem.hh
@@ -0,0 +1,568 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup TwoPNCMinTests
+ * \brief Problem where water is injected in a for flushing precipitated salt clogging a gas reservoir.
+ */
+#ifndef DUMUX_SALINIZATION_PROBLEM_HH
+#define DUMUX_SALINIZATION_PROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/elementsolution.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/discretization/cctpfa.hh>
+#include <dumux/discretization/box.hh>
+#include <dumux/porousmediumflow/2pncmin/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/material/fluidsystems/brineair.hh>
+
+#include <dumux/material/components/nacl.hh>
+#include <dumux/material/components/granite.hh>
+#include <dumux/material/solidsystems/compositionalsolidphase.hh>
+
+#include "spatialparams.hh"
+
+namespace Dumux {
+/*!
+ * \ingroup TwoPNCMinTests
+ * \brief Problem where brine is evaporating at the top boundary. The system is closed the remaining boundaries.
+ */
+template <class TypeTag>
+class SalinizationProblem;
+
+namespace Properties {
+// Create new type tags
+namespace TTag {
+struct Salinization { using InheritsFrom = std::tuple<TwoPNCMinNI>; };
+struct SalinizationBox { using InheritsFrom = std::tuple<Salinization, BoxModel>; };
+struct SalinizationCCTpfa { using InheritsFrom = std::tuple<Salinization, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the grid type
+SET_TYPE_PROP(Salinization, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<double, 2> >);
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::Salinization> { using type = SalinizationProblem<TypeTag>; };
+
+// Set fluid configuration
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::Salinization>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::BrineAir<Scalar, Components::H2O<Scalar>>;
+};
+
+template<class TypeTag>
+struct SolidSystem<TypeTag, TTag::Salinization>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using ComponentOne = Components::NaCl<Scalar>;
+ using ComponentTwo = Components::Granite<Scalar>;
+ static constexpr int numInertComponents = 1;
+ using type = SolidSystems::CompositionalSolidPhase<Scalar, ComponentOne, ComponentTwo, numInertComponents>;
+};
+
+// Set the spatial parameters
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::Salinization>
+{
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = SalinizationSpatialParams<FVGridGeometry, Scalar>;
+};
+
+// Set properties here to override the default property settings
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::Salinization> { static constexpr int value = 1; }; //!< Replace gas balance by total mass balance
+template<class TypeTag>
+struct Formulation<TypeTag, TTag::Salinization>
+{ static constexpr auto value = TwoPFormulation::p0s1; };
+
+} // end namespace Properties
+
+/*!
+ * \ingroup TwoPNCMinTests
+ * \brief Problem where water is evaporating at the top of a porous media filled
+ * container saturated with brine and air, which causes precipitation at the top.
+ *
+ * Problem with a porous media in a container, which is open to the atmosphere
+ * at the top boundary. The container has dimensions of 0.2m by 0.2m. Neumann
+ * no-flow boundaries are applied at the left, right and bottom boundary.
+ * The grid is refined towards the upper boundary to capter the relevant processes.
+ *
+ * Initially the porous medium is 50 % saturated with brine. Evaporation takes
+ * place at the top boundary and hence the temperature and liquid saturation
+ * decreases first at the top, then in the whole system, whereas the sodium
+ * chloride (NaCl) concentration increases. This results in precipitaion of NaCl
+ * at the top as the solubility limit is exceeded. Due to the low liquid saturation
+ * the top after some time, top temperature rises again.
+ *
+ * The model uses mole fractions of dissolved components and volume fractions of
+ * precipitated salt as primary variables. Make sure that the according units
+ * are used in the problem set-up.
+ *
+ * This problem uses the \ref TwoPNCMinModel.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./test_2pncminni_salinization</tt>
+ */
+template <class TypeTag>
+class SalinizationProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using SolidSystem = GetPropType<TypeTag, Properties::SolidSystem>;
+
+ enum
+ {
+ // primary variable indices
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx,
+
+ // component indices
+ // TODO: using xwNaClIdx as privaridx works here, but
+ // looks like magic. Can this be done differently??
+ xwNaClIdx = FluidSystem::NaClIdx,
+ precipNaClIdx = FluidSystem::numComponents,
+
+ // Indices of the components
+ H2OIdx = FluidSystem::H2OIdx,
+ NaClIdx = FluidSystem::NaClIdx,
+ AirIdx = FluidSystem::AirIdx,
+
+ // Indices of the phases
+ liquidPhaseIdx = FluidSystem::liquidPhaseIdx,
+ gasPhaseIdx = FluidSystem::gasPhaseIdx,
+
+ // index of the solid phase
+ sPhaseIdx = SolidSystem::comp0Idx,
+
+
+ // Index of the primary component of G and L phase
+ conti0EqIdx = Indices::conti0EqIdx, //water component
+ conti1EqIdx = Indices::conti0EqIdx + 1, //air component
+ precipNaClEqIdx = Indices::conti0EqIdx + FluidSystem::numComponents,
+ energyEqIdx = Indices::energyEqIdx,
+
+ // Phase State
+ bothPhases = Indices::bothPhases,
+
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ };
+
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using GlobalPosition = typename SubControlVolume::GlobalPosition;
+ using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+public:
+ SalinizationProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ //Fluidsystem
+ nTemperature_ = getParam<int>("FluidSystem.NTemperature");
+ nPressure_ = getParam<int>("FluidSystem.NPressure");
+ pressureLow_ = getParam<Scalar>("FluidSystem.PressureLow");
+ pressureHigh_ = getParam<Scalar>("FluidSystem.PressureHigh");
+ temperatureLow_ = getParam<Scalar>("FluidSystem.TemperatureLow");
+ temperatureHigh_ = getParam<Scalar>("FluidSystem.TemperatureHigh");
+ name_ = getParam<std::string>("Problem.Name");
+
+ //problem
+ name_ = getParam<std::string>("Problem.Name");
+ temperature_ = getParam<Scalar>("Problem.Temperature");
+
+ //inital conditions
+ initPressure_ = getParam<Scalar>("Problem.InitialPressure");
+ initGasSaturation_ = getParam<Scalar>("Problem.InitialGasSaturation");
+ initSalinity_ = getParam<Scalar>("Problem.InitialSalinity");
+
+ unsigned int codim = GET_PROP_TYPE(TypeTag, FVGridGeometry)::discMethod == DiscretizationMethod::box ? dim : 0;
+ permeability_.resize(fvGridGeometry->gridView().size(codim));
+
+ FluidSystem::init(/*Tmin=*/temperatureLow_,
+ /*Tmax=*/temperatureHigh_,
+ /*nT=*/nTemperature_,
+ /*pmin=*/pressureLow_,
+ /*pmax=*/pressureHigh_,
+ /*np=*/nPressure_);
+ }
+
+ /*!
+ * \brief The current time.
+ */
+ void setTime( Scalar time )
+ {
+ time_ = time;
+ }
+
+ /*!
+ * \brief The time step size.
+ *
+ * This is used to calculate the source term.
+ */
+ void setTimeStepSize( Scalar timeStepSize )
+ {
+ timeStepSize_ = timeStepSize;
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const std::string& name() const
+ { return name_; }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature() const
+ { return temperature_; }
+
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ */
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+ {
+ BoundaryTypes bcTypes;
+
+ // default to Neumann
+ bcTypes.setAllNeumann();
+
+ return bcTypes;
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables priVars(0.0);
+ priVars.setState(bothPhases);
+
+ return priVars;
+ }
+
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * This is the method for the case where the Neumann condition is
+ * potentially solution dependent
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param elemVolVars All volume variables for the element
+ * \param elemFluxVarsCache Flux variables caches for all faces in stencil
+ * \param scvf The sub control volume face
+ *
+ * Negative values mean influx.
+ * E.g. for the mass balance that would be the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
+ */
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
+ {
+ PrimaryVariables values(0.0);
+
+ const auto& globalPos = scvf.ipGlobal();
+ const auto& volVars = elemVolVars[scvf.insideScvIdx()];
+ const Scalar hmax = this->fvGridGeometry().bBoxMax()[1];
+
+ static const Scalar temperatureRef = getParam<Scalar>("FreeFlow.RefTemperature");
+
+
+ if(globalPos[1] > hmax - eps_)
+ {
+ // get free-flow properties:
+ static const Scalar moleFracRefH2O = getParam<Scalar>("FreeFlow.RefMoleFracH2O");
+ static const Scalar boundaryLayerThickness = getParam<Scalar>("FreeFlow.BoundaryLayerThickness");
+ static const Scalar massTransferCoefficient = getParam<Scalar>("FreeFlow.MassTransferCoefficient");
+
+ // get porous medium values:
+ Scalar moleFracH2OInside = volVars.moleFraction(gasPhaseIdx, H2OIdx);
+ Scalar referencePermeability_ = getParam<Scalar>("SpatialParams.referencePermeability", 2.23e-14);
+
+ // calculate fluxes
+ // liquid phase
+ Scalar evaporationRateMole = 0;
+ if(moleFracH2OInside - moleFracRefH2O > 0)
+ {
+ evaporationRateMole = massTransferCoefficient
+ * volVars.diffusionCoefficient(gasPhaseIdx, H2OIdx)
+ * (moleFracH2OInside - moleFracRefH2O)
+ / boundaryLayerThickness
+ * volVars.molarDensity(gasPhaseIdx);
+ }
+ else
+ {
+ evaporationRateMole = massTransferCoefficient
+ * volVars.diffusionCoefficient(gasPhaseIdx, H2OIdx)
+ * (moleFracH2OInside - moleFracRefH2O)
+ / boundaryLayerThickness
+ * 1.2;
+
+ }
+
+ values[conti0EqIdx] = evaporationRateMole;
+
+ // gas phase
+ // gas flows in
+ if (volVars.pressure(gasPhaseIdx) - 1e5 > 0) {
+ values[conti1EqIdx] = (volVars.pressure(gasPhaseIdx) - 1e5)
+ /(globalPos - fvGeometry.scv(scvf.insideScvIdx()).center()).two_norm()
+ *volVars.mobility(gasPhaseIdx)
+ *referencePermeability_
+ *volVars.molarDensity(gasPhaseIdx)
+ *volVars.moleFraction(gasPhaseIdx, AirIdx);
+ }
+ //gas flows out
+ else {
+ values[conti1EqIdx] = (volVars.pressure(gasPhaseIdx) - 1e5)
+ /(globalPos - fvGeometry.scv(scvf.insideScvIdx()).center()).two_norm()
+ *volVars.mobility(gasPhaseIdx)
+ *referencePermeability_
+ *volVars.molarDensity(gasPhaseIdx) * (1-moleFracRefH2O);
+ }
+
+
+ // energy fluxes
+ values[energyEqIdx] = FluidSystem::componentEnthalpy(volVars.fluidState(), gasPhaseIdx, H2OIdx) * values[conti0EqIdx] * FluidSystem::molarMass(H2OIdx);
+
+ values[energyEqIdx] += FluidSystem::componentEnthalpy(volVars.fluidState(), gasPhaseIdx, AirIdx)* values[conti1EqIdx] * FluidSystem::molarMass(AirIdx);
+
+ values[energyEqIdx] += FluidSystem::thermalConductivity(elemVolVars[scvf.insideScvIdx()].fluidState(), gasPhaseIdx) * (volVars.temperature() - temperatureRef)/boundaryLayerThickness;
+
+ }
+ return values;
+ }
+
+
+
+
+
+
+ /*!
+ * \brief Evaluates the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables priVars(0.0);
+ priVars.setState(bothPhases);
+ Scalar density = 1000.00; //FluidSystem::density(, liquidPhaseIdx);
+
+ priVars[pressureIdx] = initPressure_ - density*9.81*globalPos[dimWorld-1];
+ priVars[switchIdx] = initGasSaturation_; // Sg primary variable
+ priVars[xwNaClIdx] = massToMoleFrac_(initSalinity_); // mole fraction
+ priVars[precipNaClIdx] = 0.0; // [kg/m^3]
+ priVars[energyEqIdx] = temperature_; // [K]
+
+ return priVars;
+ }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluates the source term for all phases within a given
+ * sub-controlvolume.
+ *
+ * This is the method for the case where the source term is
+ * potentially solution dependent and requires some quantities that
+ * are specific to the fully-implicit method.
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param elemVolVars All volume variables for the element
+ * \param scv The subcontrolvolume
+ *
+ * For this method, the \a values parameter stores the conserved quantity rate
+ * generated or annihilated per volume unit. Positive values mean
+ * that the conserved quantity is created, negative ones mean that it vanishes.
+ * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
+ */
+ NumEqVector source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ {
+ NumEqVector source(0.0);
+
+ const auto& volVars = elemVolVars[scv];
+
+ Scalar moleFracNaCl_wPhase = volVars.moleFraction(liquidPhaseIdx, NaClIdx);
+ Scalar massFracNaCl_Max_wPhase = this->spatialParams().solubilityLimit();
+ Scalar moleFracNaCl_Max_wPhase = massToMoleFrac_(massFracNaCl_Max_wPhase);
+ Scalar saltPorosity = this->spatialParams().minimalPorosity(element, scv);
+
+ // precipitation of amount of salt whic hexeeds the solubility limit
+ using std::abs;
+ Scalar precipSalt = volVars.porosity() * volVars.molarDensity(liquidPhaseIdx)
+ * volVars.saturation(liquidPhaseIdx)
+ * abs(moleFracNaCl_wPhase - moleFracNaCl_Max_wPhase);
+ if (moleFracNaCl_wPhase < moleFracNaCl_Max_wPhase)
+ precipSalt *= -1;
+
+ // make sure we don't dissolve more salt than previously precipitated
+ if (precipSalt*timeStepSize_ + volVars.solidVolumeFraction(sPhaseIdx)* volVars.solidComponentMolarDensity(sPhaseIdx)< 0)
+ precipSalt = -volVars.solidVolumeFraction(sPhaseIdx)* volVars.solidComponentMolarDensity(sPhaseIdx)/timeStepSize_;
+
+ // make sure there is still pore space available for precipitation
+ if (volVars.solidVolumeFraction(sPhaseIdx) >= this->spatialParams().referencePorosity(element, scv) - saltPorosity && precipSalt > 0)
+ precipSalt = 0;
+
+ source[conti0EqIdx + NaClIdx] += -precipSalt;
+ source[precipNaClEqIdx] += precipSalt;
+ return source;
+
+ }
+
+ /*!
+ * \brief Adds additional VTK output data to the VTKWriter.
+ *
+ * Function is called by the output module on every write.
+ */
+
+ const std::vector<Scalar>& getPermeability()
+ {
+ return permeability_;
+ }
+
+ void updateVtkOutput(const SolutionVector& curSol)
+ {
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ const auto elemSol = elementSolution(element, curSol, this->fvGridGeometry());
+
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ VolumeVariables volVars;
+ volVars.update(elemSol, *this, element, scv);
+ const auto dofIdxGlobal = scv.dofIndex();
+ permeability_[dofIdxGlobal] = volVars.permeability();
+ }
+ }
+ }
+
+ Scalar extrusionFactorAtPos(const GlobalPosition &globalPos) const
+ {
+ // As a default, i.e. if the user's problem does not overload
+ // any extrusion factor method, return 1.0
+ return 0.054977871437821;
+ }
+
+
+
+private:
+
+ /*!
+ * \brief Returns the molality of NaCl (mol NaCl / kg water) for a given mole fraction.
+ *
+ * \param XwNaCl The XwNaCl [kg NaCl / kg solution]
+ */
+ static Scalar massToMoleFrac_(Scalar XwNaCl)
+ {
+ const Scalar Mw = 18.015e-3; //FluidSystem::molarMass(H2OIdx); /* molecular weight of water [kg/mol] */ //TODO use correct link to FluidSyswem later
+ const Scalar Ms = 58.44e-3; //FluidSystem::molarMass(NaClIdx); /* molecular weight of NaCl [kg/mol] */
+
+ const Scalar X_NaCl = XwNaCl;
+ /* XwNaCl: conversion from mass fraction to mol fraction */
+ auto xwNaCl = -Mw * X_NaCl / ((Ms - Mw) * X_NaCl - Ms);
+ return xwNaCl;
+ }
+
+
+ std::string name_;
+
+ Scalar initPressure_;
+ Scalar initGasSaturation_;
+ Scalar initSalinity_;
+
+ Scalar temperature_;
+
+ Scalar pressureLow_, pressureHigh_;
+ Scalar temperatureLow_, temperatureHigh_;
+ int nTemperature_;
+ int nPressure_;
+
+ Scalar time_ = 0.0;
+ Scalar timeStepSize_ = 0.0;
+ static constexpr Scalar eps_ = 1e-6;
+
+ std::vector<double> permeability_;
+
+ Dumux::GnuplotInterface<double> gnuplot_;
+ Dumux::GnuplotInterface<double> gnuplot2_;
+ std::vector<double> x_;
+ std::vector<double> y_;
+ std::vector<double> y2_;
+
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/2pncmin/implicit/nonisothermal/spatialparams.hh b/test/porousmediumflow/2pncmin/implicit/nonisothermal/spatialparams.hh
new file mode 100644
index 0000000000..1a47d8bdd2
--- /dev/null
+++ b/test/porousmediumflow/2pncmin/implicit/nonisothermal/spatialparams.hh
@@ -0,0 +1,181 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup TwoPNCMinTests
+ * \brief Spatial parameters for the salinization problem, where evaporation
+ * from a porous medium saturated wit brine and air leads to precipitation of salt.
+ */
+
+#ifndef DUMUX_SALINIZATION_SPATIAL_PARAMETERS_HH
+#define DUMUX_SALINIZATION_SPATIAL_PARAMETERS_HH
+
+#include <dumux/io/gnuplotinterface.hh>
+#include <dumux/io/plotmateriallaw.hh>
+#include <dumux/material/spatialparams/fv.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+#include <dumux/material/fluidmatrixinteractions/porosityprecipitation.hh>
+#include <dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup TwoPNCMinTests
+ * \brief Spatial parameters for the salinization problem, where evaporation
+ * from a porous medium saturated wit brine and air leads to precipitation of salt.
+ */
+template<class FVGridGeometry, class Scalar>
+class SalinizationSpatialParams
+: public FVSpatialParams<FVGridGeometry, Scalar,
+ SalinizationSpatialParams<FVGridGeometry, Scalar>>
+{
+ using GridView = typename FVGridGeometry::GridView;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using ParentType = FVSpatialParams<FVGridGeometry, Scalar,
+ SalinizationSpatialParams<FVGridGeometry, Scalar>>;
+
+ using EffectiveLaw = RegularizedVanGenuchten<Scalar>;
+
+ using GlobalPosition = typename SubControlVolume::GlobalPosition;
+
+public:
+ // type used for the permeability (i.e. tensor or scalar)
+ using PermeabilityType = Scalar;
+ //! Export the material law type used
+ using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
+ using MaterialLawParams = typename MaterialLaw::Params;
+
+ SalinizationSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ solubilityLimit_ = getParam<Scalar>("SpatialParams.SolubilityLimit", 0.26);
+ referencePorosity_ = getParam<Scalar>("SpatialParams.referencePorosity", 0.11);
+ referencePermeability_ = getParam<Scalar>("SpatialParams.referencePermeability", 2.23e-14);
+ irreducibleLiqSat_ = getParam<Scalar>("SpatialParams.IrreducibleLiqSat", 0.2);
+ irreducibleGasSat_ = getParam<Scalar>("SpatialParams.IrreducibleGasSat", 1e-3);
+ vgAlpha_ = getParam<Scalar>("SpatialParams.VGAlpha", 1.5);
+ vgn_ = getParam<Scalar>("SpatialParams.VGn", 4.0);
+
+ plotFluidMatrixInteractions_ = getParam<bool>("Output.PlotFluidMatrixInteractions");
+
+ //Van Genuchen parameters
+ // residual saturations
+ materialParams_.setSwr(irreducibleLiqSat_);
+ materialParams_.setSnr(irreducibleGasSat_);
+
+ //Van Genuchen parameters
+ materialParams_.setVgAlpha(vgAlpha_ );
+ materialParams_.setVgn(vgn_);
+ }
+
+ /*!
+ * \brief Defines the minimum porosity \f$[-]\f$ distribution
+ *
+ * \param element The finite element
+ * \param scv The sub-control volume
+ */
+ Scalar minimalPorosity(const Element& element, const SubControlVolume &scv) const
+ { return 1e-5; }
+
+ /*!
+ * \brief Defines the volume fraction of the inert component
+ *
+ * \param globalPos The global position in the domain
+ * \param compIdx The index of the inert solid component
+ */
+ template<class SolidSystem>
+ Scalar inertVolumeFractionAtPos(const GlobalPosition& globalPos, int compIdx) const
+ { return 1.0-referencePorosity_; }
+
+ /*!
+ * \brief Defines the reference porosity \f$[-]\f$ distribution.
+ *
+ * This is the porosity of the porous medium without any of the
+ * considered solid phases.
+ *
+ * \param element The finite element
+ * \param scv The sub-control volume
+ */
+ Scalar referencePorosity(const Element& element, const SubControlVolume &scv) const
+ { return referencePorosity_; }
+
+ /*! Intrinsic permeability tensor K \f$[m^2]\f$ depending
+ * on the position in the domain
+ *
+ * \param element The finite volume element
+ * \param scv The sub-control volume
+ * \param elemSol The element solution
+ *
+ * Solution dependent permeability function
+ */
+ template<class ElementSolution>
+ PermeabilityType permeability(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ {
+ auto priVars = evalSolution(element, element.geometry(), elemSol, scv.center(), /*ignoreState=*/true);
+
+ Scalar sumPrecipitates = priVars[/*numComp*/3];
+
+ using std::max;
+ const auto poro = max(/*minPoro*/1e-5, referencePorosity_ - sumPrecipitates);
+ return permLaw_.evaluatePermeability(referencePermeability_, referencePorosity_, poro);
+ }
+
+ // the solubility limit of NaCl
+ Scalar solubilityLimit() const
+ { return solubilityLimit_; }
+
+ Scalar theta(const SubControlVolume &scv) const
+ { return 10.0; }
+
+ // return the brooks-corey context depending on the position
+ const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
+ { return materialParams_; }
+
+ // define which phase is to be considered as the wetting phase
+ template<class FluidSystem>
+ int wettingPhaseAtPos(const GlobalPosition& globalPos) const
+ { return FluidSystem::H2OIdx; }
+
+private:
+
+ MaterialLawParams materialParams_;
+
+ PermeabilityKozenyCarman<PermeabilityType> permLaw_;
+
+ Scalar solubilityLimit_;
+ Scalar referencePorosity_;
+ PermeabilityType referencePermeability_ = 0.0;
+ Scalar irreducibleLiqSat_;
+ Scalar irreducibleGasSat_;
+ Scalar vgAlpha_;
+ Scalar vgn_ ;
+
+ bool plotFluidMatrixInteractions_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/references/test_2pncminni_salinization-reference.vtu b/test/references/test_2pncminni_salinization-reference.vtu
new file mode 100644
index 0000000000..1d14a543b1
--- /dev/null
+++ b/test/references/test_2pncminni_salinization-reference.vtu
@@ -0,0 +1,139 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+ <UnstructuredGrid>
+ <Piece NumberOfCells="10" NumberOfPoints="18">
+ <PointData Scalars="S_liq" Vectors="velocity_liq (m/s)">
+ <DataArray type="Float32" Name="S_liq" NumberOfComponents="1" format="ascii">
+ 0.371972 0.371972 0.354829 0.354829 0.371972 0.354829 0.328148 0.328148 0.328148 0.315167 0.315167 0.315167
+ 0.300425 0.300425 0.300425 0.283783 0.283783 0.283783
+ </DataArray>
+ <DataArray type="Float32" Name="p_liq" NumberOfComponents="1" format="ascii">
+ 92441.8 92441.8 92339.1 92339.1 92441.8 92339.1 92142.7 92142.7 92142.7 92022.1 92022.1 92022.1
+ 91851.9 91851.9 91851.9 91583.7 91583.7 91583.7
+ </DataArray>
+ <DataArray type="Float32" Name="rho_liq" NumberOfComponents="1" format="ascii">
+ 1012.93 1012.93 1012.93 1012.93 1012.93 1012.93 1012.94 1012.94 1012.94 1013.29 1013.29 1013.29
+ 1021.51 1021.51 1021.51 1200.29 1200.29 1200.29
+ </DataArray>
+ <DataArray type="Float32" Name="mob_liq" NumberOfComponents="1" format="ascii">
+ 6.3638 6.3638 4.25016 4.25016 6.3638 4.25016 1.94427 1.94427 1.94427 1.19814 1.19814 1.19814
+ 0.594915 0.594915 0.594915 0.111125 0.111125 0.111125
+ </DataArray>
+ <DataArray type="Float32" Name="S_gas" NumberOfComponents="1" format="ascii">
+ 0.628028 0.628028 0.645171 0.645171 0.628028 0.645171 0.671852 0.671852 0.671852 0.684833 0.684833 0.684833
+ 0.699575 0.699575 0.699575 0.716217 0.716217 0.716217
+ </DataArray>
+ <DataArray type="Float32" Name="p_gas" NumberOfComponents="1" format="ascii">
+ 100002 100002 100002 100002 100002 100002 100002 100002 100002 100002 100002 100002
+ 100002 100002 100002 100001 100001 100001
+ </DataArray>
+ <DataArray type="Float32" Name="rho_gas" NumberOfComponents="1" format="ascii">
+ 1.21014 1.21014 1.21013 1.21013 1.21014 1.21013 1.21012 1.21012 1.21012 1.21012 1.21012 1.21012
+ 1.21014 1.21014 1.21014 1.21073 1.21073 1.21073
+ </DataArray>
+ <DataArray type="Float32" Name="mob_gas" NumberOfComponents="1" format="ascii">
+ 39281.7 39281.7 41576.5 41576.5 39281.7 41576.5 45236.5 45236.5 45236.5 47150.2 47150.2 47150.2
+ 49764.6 49764.6 49764.6 52753.5 52753.5 52753.5
+ </DataArray>
+ <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
+ 7560.19 7560.19 7662.75 7662.75 7560.19 7662.75 7859.09 7859.09 7859.09 7979.69 7979.69 7979.69
+ 8149.9 8149.9 8149.9 8417.42 8417.42 8417.42
+ </DataArray>
+ <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
+ 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
+ 0.4 0.4 0.4 0.399824 0.399824 0.399824
+ </DataArray>
+ <DataArray type="Float32" Name="x^H2O_liq" NumberOfComponents="1" format="ascii">
+ 0.993731 0.993731 0.993731 0.993731 0.993731 0.993731 0.993726 0.993726 0.993726 0.993565 0.993565 0.993565
+ 0.989756 0.989756 0.989756 0.902258 0.902258 0.902258
+ </DataArray>
+ <DataArray type="Float32" Name="x^Air_liq" NumberOfComponents="1" format="ascii">
+ 1.69669e-05 1.69669e-05 1.69666e-05 1.69666e-05 1.69669e-05 1.69666e-05 1.69657e-05 1.69657e-05 1.69657e-05 1.69654e-05 1.69654e-05 1.69654e-05
+ 1.69661e-05 1.69661e-05 1.69661e-05 1.69884e-05 1.69884e-05 1.69884e-05
+ </DataArray>
+ <DataArray type="Float32" Name="x^NaCl_liq" NumberOfComponents="1" format="ascii">
+ 0.00625221 0.00625221 0.00625235 0.00625235 0.00625221 0.00625235 0.00625709 0.00625709 0.00625709 0.0064178 0.0064178 0.0064178
+ 0.0102271 0.0102271 0.0102271 0.0977251 0.0977251 0.0977251
+ </DataArray>
+ <DataArray type="Float32" Name="rhoMolar_liq" NumberOfComponents="1" format="ascii">
+ 55447.9 55447.9 55447.9 55447.9 55447.9 55447.9 55447.8 55447.8 55447.8 55447.1 55447.1 55447.1
+ 55430 55430 55430 54642.7 54642.7 54642.7
+ </DataArray>
+ <DataArray type="Float32" Name="x^H2O_gas" NumberOfComponents="1" format="ascii">
+ 0.0149431 0.0149431 0.0149438 0.0149438 0.0149431 0.0149438 0.014946 0.014946 0.014946 0.0149445 0.0149445 0.0149445
+ 0.014888 0.014888 0.014888 0.0135726 0.0135726 0.0135726
+ </DataArray>
+ <DataArray type="Float32" Name="x^Air_gas" NumberOfComponents="1" format="ascii">
+ 0.985057 0.985057 0.985056 0.985056 0.985057 0.985056 0.985054 0.985054 0.985054 0.985056 0.985056 0.985056
+ 0.985112 0.985112 0.985112 0.986427 0.986427 0.986427
+ </DataArray>
+ <DataArray type="Float32" Name="x^NaCl_gas" NumberOfComponents="1" format="ascii">
+ -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0 -0
+ -0 -0 -0 -0 -0 -0
+ </DataArray>
+ <DataArray type="Float32" Name="rhoMolar_gas" NumberOfComponents="1" format="ascii">
+ 42.0241 42.0241 42.024 42.024 42.0241 42.024 42.0236 42.0236 42.0236 42.0234 42.0234 42.0234
+ 42.0233 42.0233 42.0233 42.0228 42.0228 42.0228
+ </DataArray>
+ <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
+ 3 3 3 3 3 3 3 3 3 3 3 3
+ 3 3 3 3 3 3
+ </DataArray>
+ <DataArray type="Float32" Name="precipitateVolumeFraction^NaCl" NumberOfComponents="1" format="ascii">
+ 0 0 0 0 0 0 0 0 0 0 0 0
+ 0 0 0 0.00017601 0.00017601 0.00017601
+ </DataArray>
+ <DataArray type="Float32" Name="T" NumberOfComponents="1" format="ascii">
+ 286.208 286.208 286.209 286.209 286.208 286.209 286.211 286.211 286.211 286.212 286.212 286.212
+ 286.213 286.213 286.213 286.214 286.214 286.214
+ </DataArray>
+ <DataArray type="Float32" Name="velocity_liq (m/s)" NumberOfComponents="3" format="ascii">
+ -6.7602e-22 2.38612e-09 0 -3.79683e-22 2.38612e-09 0 -6.3085e-22 4.64661e-09 0 -3.23156e-22 4.64661e-09 0
+ -8.3345e-23 2.38612e-09 0 -1.5462e-23 4.64661e-09 0 -3.42343e-22 8.41592e-09 0 -2.7727e-22 8.41592e-09 0
+ -2.12196e-22 8.41592e-09 0 -4.95383e-22 1.05518e-08 0 -3.18852e-22 1.05518e-08 0 -1.42321e-22 1.05518e-08 0
+ -1.73772e-22 1.1066e-08 0 -5.82074e-23 1.1066e-08 0 5.73568e-23 1.1066e-08 0 1.94852e-22 1.0953e-08 0
+ 1.33006e-22 1.0953e-08 0 7.11605e-23 1.0953e-08 0
+ </DataArray>
+ <DataArray type="Float32" Name="velocity_gas (m/s)" NumberOfComponents="3" format="ascii">
+ -9.77477e-18 -2.47197e-09 0 -5.91631e-18 -2.47197e-09 0 -6.62494e-18 -4.83592e-09 0 -6.20143e-18 -4.83592e-09 0
+ -2.05785e-18 -2.47197e-09 0 -5.77792e-18 -4.83592e-09 0 -5.75993e-18 -7.66282e-09 0 -7.17523e-18 -7.66282e-09 0
+ -8.59053e-18 -7.66282e-09 0 -6.91496e-18 2.53119e-08 0 -1.30386e-17 2.53119e-08 0 -1.91623e-17 2.53119e-08 0
+ 5.96753e-18 1.0712e-06 0 -8.72983e-18 1.0712e-06 0 -2.34272e-17 1.0712e-06 0 -3.41432e-19 2.08366e-06 0
+ -2.79445e-18 2.08366e-06 0 -5.24746e-18 2.08366e-06 0
+ </DataArray>
+ <DataArray type="Float32" Name="Permeability" NumberOfComponents="1" format="ascii">
+ 1e-13 1e-13 1e-13 1e-13 1e-13 1e-13 1e-13 1e-13 1e-13 1e-13 1e-13 1e-13
+ 9.99523e-14 9.99523e-14 9.99523e-14 9.98571e-14 9.98571e-14 9.98571e-14
+ </DataArray>
+ </PointData>
+ <CellData Scalars="process rank">
+ <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
+ 0 0 0 0 0 0 0 0 0 0
+ </DataArray>
+ </CellData>
+ <Points>
+ <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+ 0 0 0 0.01 0 0 0 0.0075 0 0.01 0.0075 0
+ 0.02 0 0 0.02 0.0075 0 0 0.015 0 0.01 0.015 0
+ 0.02 0.015 0 0 0.016978 0 0.01 0.016978 0 0.02 0.016978 0
+ 0 0.0186264 0 0.01 0.0186264 0 0.02 0.0186264 0 0 0.02 0
+ 0.01 0.02 0 0.02 0.02 0
+ </DataArray>
+ </Points>
+ <Cells>
+ <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+ 0 1 3 2 1 4 5 3 2 3 7 6
+ 3 5 8 7 6 7 10 9 7 8 11 10
+ 9 10 13 12 10 11 14 13 12 13 16 15
+ 13 14 17 16
+ </DataArray>
+ <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+ 4 8 12 16 20 24 28 32 36 40
+ </DataArray>
+ <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+ 9 9 9 9 9 9 9 9 9 9
+ </DataArray>
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+</VTKFile>
--
GitLab