Commit f9e31fee authored by Maziar Veyskarami's avatar Maziar Veyskarami Committed by Timo Koch
Browse files

[test][porenetworkflow] add 1p non-creeping flow test

parent 05668397
add_subdirectory(nonisothermal)
add_subdirectory(noncreepingflow)
add_input_file_links()
dune_symlink_to_source_files(FILES grids)
......
add_subdirectory(nonisothermal)
add_input_file_links()
dune_add_test(NAME test_pnm_1p_noncreeping_flow
SOURCES main.cc
LABELS porenetworkflow
COMPILE_DEFINITIONS ISOTHERMAL=1
CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )"
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_pnm_1p_noncreepingflow-reference.vtp
${CMAKE_CURRENT_BINARY_DIR}/test_pnm_1p_noncreepingflow-00000.vtp
--command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_1p_noncreeping_flow")
// -*- 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
*
* \brief test for the pore-network model with non-creeping flow
*/
#include <config.h>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/porenetwork/common/pnmvtkoutputmodule.hh>
#include <dumux/porenetwork/common/boundaryflux.hh>
#include <dumux/io/grid/porenetwork/gridmanager.hh>
#include "properties.hh"
int main(int argc, char** argv)
{
using namespace Dumux;
using TypeTag = Properties::TTag::PNMOnePNonCreepingProblem;
// 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 the command line arguments and input file
////////////////////////////////////////////////////////////
// parse command line arguments
Parameters::init(argc, argv);
//////////////////////////////////////////////////////////////////////
// try to create a grid (from the given grid file or the input file)
/////////////////////////////////////////////////////////////////////
using GridManager = Dumux::PoreNetwork::GridManager<3>;
GridManager gridManager;
gridManager.init();
// we compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
auto gridData = gridManager.getGridData();
// create the finite volume grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
gridGeometry->update(*gridData);
// the spatial parameters
using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
auto spatialParams = std::make_shared<SpatialParams>(gridGeometry);
// the problem (boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry, spatialParams);
// the solution vector
using GridView = typename GridGeometry::GridView;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(leafGridView.size(GridView::dimension));
problem->applyInitialSolution(x);
auto xOld = x;
// the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
// the assembler
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
// helper class to calculate the boundary flux
const auto boundaryFlux = Dumux::PoreNetwork::BoundaryFlux(*gridVariables, assembler->localResidual(), x);
// the linear solver
using LinearSolver = ILU0RestartedGMResBackend;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
nonLinearSolver.solve(x);
// output result to vtk
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
Dumux::PoreNetwork::VtkOutputModule<GridVariables, GetPropType<TypeTag, Properties::FluxVariables>, SolutionVector>
vtkWriter(*gridVariables, x, problem->name());
IOFields::initOutputModule(vtkWriter);
vtkWriter.write(0.0);
if (mpiHelper.rank() == 0)
{
std::cout << "cumulative outflux is: "
<< boundaryFlux.getFlux("max", 0, true) << std::endl;
DumuxMessage::print();
Parameters::print();
}
return 0;
}
add_input_file_links()
dune_add_test(NAME test_pnm_1pni_noncreeping_flow
SOURCES main.cc
LABELS porenetworkflow
COMPILE_DEFINITIONS ISOTHERMAL=0
CMAKE_GUARD "( dune-foamgrid_FOUND AND HAVE_UMFPACK )"
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_pnm_1pni_noncreepingflow-reference.vtp
${CMAKE_CURRENT_BINARY_DIR}/test_pnm_1pni_noncreepingflow-00014.vtp
--command "${CMAKE_CURRENT_BINARY_DIR}/test_pnm_1pni_noncreeping_flow")
// -*- 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
*
* \brief test for the non-isothermal pore-network model with non-creeping flow
*/
#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 <dumux/assembly/fvassembler.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/porenetwork/common/pnmvtkoutputmodule.hh>
#include <dumux/porenetwork/common/boundaryflux.hh>
#include <dumux/io/grid/porenetwork/gridmanager.hh>
#include "../properties.hh"
int main(int argc, char** argv)
{
using namespace Dumux;
using TypeTag = Properties::TTag::PNMOnePNonCreepingProblem;
// 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 the command line arguments and input file
////////////////////////////////////////////////////////////
// parse command line arguments
Parameters::init(argc, argv);
//////////////////////////////////////////////////////////////////////
// try to create a grid (from the given grid file or the input file)
/////////////////////////////////////////////////////////////////////
using GridManager = Dumux::PoreNetwork::GridManager<3>;
GridManager gridManager;
gridManager.init();
// we compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
auto gridData = gridManager.getGridData();
// create the finite volume grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
gridGeometry->update(*gridData);
// the spatial parameters
using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
auto spatialParams = std::make_shared<SpatialParams>(gridGeometry);
// the problem (boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry, spatialParams);
// the solution vector
using GridView = typename GridGeometry::GridView;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(leafGridView.size(GridView::dimension));
problem->applyInitialSolution(x);
auto xOld = x;
// the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
// 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");
// intialize the vtk output module
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
Dumux::PoreNetwork::VtkOutputModule<GridVariables, GetPropType<TypeTag, Properties::FluxVariables>, SolutionVector>
vtkWriter(*gridVariables, x, problem->name());
IOFields::initOutputModule(vtkWriter); //! Add model specific output fields
vtkWriter.write(0.0);
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, 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, gridGeometry, gridVariables, timeLoop, xOld);
// the linear solver
using LinearSolver = ILU0BiCGSTABBackend;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
// time loop
timeLoop->start(); do
{
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
// solve the non-linear system
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();
// write vtk output
vtkWriter.write(timeLoop->time());
// report statistics of this time step
timeLoop->reportTimeStep();
// set new dt as suggested by 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;
}
[TimeLoop]
DtInitial = 1e-3 # [s]
TEnd = 1 # [s]
[Grid]
File = ./../../grids/pnm_converted.dgf
PoreGeometry = Cube
ThroatCrossSectionShape = Square
[Problem]
Name = test_pnm_1pni_noncreepingflow # name passed to the output routines
Verbose = true
InletPressure = 1e5
OutletPressure = 0.9e5
EnableGravity = false
[Vtk]
AddVelocity = true
[Component]
SolidHeatCapacity = 1.0 # for compatibility
SolidDensity = 1.0 # for compatibility
SolidThermalConductivity = 1.0 # for compatibility
[Grid]
File = ./../grids/pnm_converted.dgf
PoreGeometry = Cube
ThroatCrossSectionShape = Square
[Problem]
Name = test_pnm_1p_noncreepingflow # name passed to the output routines
InletPressure = 1e5
OutletPressure = 0.9e5
EnableGravity = false
[Vtk]
AddVelocity = true
// -*- 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
*
* \brief The properties for the one-phase pore-network model with non-creeping flow.
*/
#ifndef DUMUX_PNM1P_NONCREEPING_PROPERTIES_HH
#define DUMUX_PNM1P_NONCREEPING_PROPERTIES_HH
#include <dune/foamgrid/foamgrid.hh>
#include <dumux/common/properties.hh>
#include <dumux/porenetwork/1p/model.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh>
#include "spatialparams.hh"
#include "../properties.hh"
//////////
// Specify the properties
//////////
namespace Dumux::Properties {
namespace TTag {
struct PNMOnePNonCreepingProblem { using InheritsFrom = std::tuple<PNMOnePProblem>; };
}
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::PNMOnePNonCreepingProblem>
{
private:
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = Dumux::PoreNetwork::NonCreepingSpatialParams<GridGeometry, Scalar>;
};
//! The advection type
template<class TypeTag>
struct AdvectionType<TypeTag, TTag::PNMOnePNonCreepingProblem>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using TransmissibilityLaw = Dumux::PoreNetwork::TransmissibilityPatzekSilin<Scalar, false/*considerPoreBodyResistance*/>;
public:
using type = Dumux::PoreNetwork::NonCreepingFlow<Scalar, TransmissibilityLaw>;
};
} //end namespace Dumux::Properties
#endif
// -*- 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 SpatialParameters
* \brief Spatial parameters for a pore-network model with non-creeping flow.
*/
#ifndef DUMUX_PNM_NONCREEPING_SPATIAL_PARAMS_1P_HH
#define DUMUX_PNM_NONCREEPING_SPATIAL_PARAMS_1P_HH
#include <dumux/material/spatialparams/porenetwork/porenetworkbase.hh>
#include <dumux/porenetwork/common/poreproperties.hh>
#include <dumux/porenetwork/common/throatproperties.hh>
namespace Dumux::PoreNetwork {
template<class GridGeometry, class Scalar>
class NonCreepingSpatialParams : public BaseSpatialParams<GridGeometry, Scalar, NonCreepingSpatialParams<GridGeometry, Scalar>>
{
using ParentType = BaseSpatialParams<GridGeometry, Scalar, NonCreepingSpatialParams<GridGeometry, Scalar>>;
using GridView = typename GridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename GridGeometry::SubControlVolume;
public:
using PermeabilityType = Scalar;
using ParentType::ParentType;
template<class ElementSolutionVector>
auto poreShapeFactor(const Element& element,
const SubControlVolume& scv,
const ElementSolutionVector& elemSol) const
{
return Throat::shapeFactorCircle<Scalar>();
}
template<class ElementSolutionVector>
Scalar poreCrossSectionalArea(const Element& element,
const SubControlVolume& scv,
const ElementSolutionVector& elemSol) const
{
const auto poreRadius = this->gridGeometry().poreInscribedRadius(scv.dofIndex());
return poreRadius*poreRadius*M_PI;
}
template<class ElementSolutionVector>
Scalar poreLength(const Element& element,
const SubControlVolume& scv,
const ElementSolutionVector& elemSol) const
{
return this->gridGeometry().poreInscribedRadius(scv.dofIndex());
}
// dimensionless kinetic-energy coeffiecient which for non-creeping flow is equal to 1.0
template<class ElementSolutionVector>
Scalar kineticEnergyCoefficient(const Element& element,
const SubControlVolume& scv,
const ElementSolutionVector& elemSol) const
{ return 1.0; }
// dimensionless momentum coeffiecient which for non-creeping flow is equal to 1.0
template<class ElementSolutionVector>
Scalar momentumCoefficient(const Element& element,
const SubControlVolume& scv,
const ElementSolutionVector& elemSol) const
{ return 1.0; }
};
} // end namespace Dumux::PoreNetwork
#endif
<?xml version="1.0"?>
<VTKFile type="PolyData" version="0.1" byte_order="LittleEndian">
<PolyData>
<Piece NumberOfLines="137" NumberOfPoints="90">
<PointData Scalars="p">
<DataArray type="Float32" Name="p" NumberOfComponents="1" format="ascii">
98363.1 100000 98013 100000 97885.8 100000 100000 100000 97566.5 98397.5 97540.4 100000
97803.4 97137.9 97137.9 97137.9 97919.9 95980.5 95980.5 97137.9 97137.9 95980.5 96677.3 96426.2
96097.4 96658.6 96658.6 96039.1 97280.2 97141.1 97788 97185 97137.9 97134.6 96811.5 97072.7
97493.4 97597 97559.6 96812.3 95893.8 95937 96658.6 95659.2 93643.4 96079.3 97117.8 97157.5
97030.4 96995.5 96835.2 96648.9 97117.3 95717.9 97197.5 96783 96491.7 94963.4 95186.9 95495.2
96730.2 95718.8 96065.5 92382.2 95387.1 94739.4 93529.7 95880.2 94623.4 93148.8 91202.5 96379.5
96405.5 96397.9 96402.3 96096.2 96888.6 96648.9 96648.9 95955.3 95319 97117.3 95880.2 90000
90000 93817.6 93148.8 93148.8 90000 93529.7
</DataArray>
<DataArray type="Float32" Name="poreInscribedRadius" NumberOfComponents="1" format="ascii">
0.00010842 0.0002263 0.00018872 8.6711e-05 0.00014246 0.00017024 0.00014006 5.8137e-05 0.00017043 0.00014712 0.00018289 0.00011213
8.001e-05 0.00026652 0.00010984 0.00018149 0.0001191 0.00010862 6.7563e-05 0.00013353 8.8091e-05 0.00011873 8.6177e-05 0.00014165
9.7691e-05 0.00014071 8.5166e-05 0.00010965 8.1041e-05 0.00015832 6.1496e-05 0.00011323 0.00033347 6.7184e-05 0.00010274 0.00016438
0.00011759 8.0701e-05 0.00015113 0.00013784 9.8161e-05 0.00013 8.1223e-05 0.00014997 0.00013153 0.00020701 0.00010946 8.5023e-05
7.8721e-05 0.00017262 0.00012575 0.00010082 0.00015856 0.00012022 0.00019603 9.1747e-05 0.0001146 0.00025862 0.00012908 0.00025963
0.00010958 6.7231e-05 0.00012907 0.00015315 9.4275e-05 0.00019972 7.1623e-05 0.00014025 0.00011092 0.00015744 5.0362e-05 0.00022549
0.00010874 9.8868e-05 0.00010793 0.00015434 8.0185e-05 6.8354e-05 0.00014551 0.00014525 0.00017218 0.00012177 9.0671e-05 0.00015715
0.00012027 0.00015784 0.00022979 0.00010367 5.3866e-05 0.00024136
</DataArray>
<DataArray type="Float32" Name="coordinationNumber" NumberOfComponents="1" format="ascii">
6 1 4 1 3 1 5 2 6 3 4 1
2 4 2 3 2 6 1 1 1 1 4 2
3 2 1 2 3 7 2 4 3 4 6 2
3 4 3 5 5 3 4 6 4 5 4 4
4 4 5 7 3 3 2 3 2 5 2 4
3 2 2 2 3 4 5 4 4 6 2 5
4 3 2 2 2 1 1 4 2 1 1 3
1 2 1 1 1 1
</DataArray>
<DataArray type="Float32" Name="poreLabel" NumberOfComponents="1" format="ascii">
-1 2 -1 2 -1 2 2 -1 -1 -1 -1 2
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 3
3 -1 -1 -1 3 -1
</DataArray>
</PointData>
<CellData Scalars="process rank" Vectors="velocity_liq (m/s)">
<DataArray type="Float32" Name="velocity_liq (m/s)" NumberOfComponents="3" format="ascii">
0.779514 0.779514 -0.779514 1.16371 -0 -1.16371 1.1927 -1.1927 -0 1.2201 -0 -0