Commit 9339340e authored by Dennis Gläser's avatar Dennis Gläser Committed by Timo Koch
Browse files

[test][multidomain][facet] Add fracture test with analytical solution

parent 32add64f
add_subdirectory(analytical)
add_subdirectory(threedomain)
dune_symlink_to_source_files(FILES "grids" "facetcoupling_1p1p.input" "convergencetest.py")
# executable for tpfa tests
add_executable(test_facetcoupling_tpfa_1p1p EXCLUDE_FROM_ALL test_facetcoupling_fv_1p1p.cc)
target_compile_definitions(test_facetcoupling_tpfa_1p1p PUBLIC BULKTYPETAG=OnePBulkTpfa LOWDIMTYPETAG=OnePLowDimTpfa)
# exectuable for box tests
add_executable(test_facetcoupling_box_1p1p EXCLUDE_FROM_ALL test_facetcoupling_fv_1p1p.cc)
target_compile_definitions(test_facetcoupling_box_1p1p PUBLIC BULKTYPETAG=OnePBulkBox LOWDIMTYPETAG=OnePLowDimBox)
dune_add_test(NAME test_facet_1p1p_tpfa_convergence
CMAKE_GUARD "( dune-foamgrid_FOUND AND dune-alugrid_FOUND AND gmsh_FOUND )"
TARGET test_facetcoupling_tpfa_1p1p
COMMAND ./convergencetest.py
CMD_ARGS test_facetcoupling_tpfa_1p1p 1e-4)
dune_add_test(NAME test_facet_1p1p_box_convergence
CMAKE_GUARD "( dune-foamgrid_FOUND AND dune-alugrid_FOUND AND gmsh_FOUND )"
TARGET test_facetcoupling_box_1p1p
COMMAND ./convergencetest.py
CMD_ARGS test_facetcoupling_box_1p1p 1e-4)
set(CMAKE_BUILD_TYPE Release)
#install sources
install(FILES
test_facetcoupling_fv_1p1p.cc
bulkProblem.hh
lowdimproblem.hh
spatialparams.hh
convergencetest.py
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/multidomain/facet/1p_1p/analytical)
// -*- 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 2 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 MixedDimension
* \ingroup MixedDimensionFacet
* \brief The problem for the bulk domain in the single-phase facet coupling test
*/
#ifndef DUMUX_TEST_TPFAFACETCOUPLING_ONEP_BULKPROBLEM_HH
#define DUMUX_TEST_TPFAFACETCOUPLING_ONEP_BULKPROBLEM_HH
#include <dune/alugrid/grid.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/multidomain/facet/box/properties.hh>
#include <dumux/multidomain/facet/cellcentered/tpfa/properties.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/porousmediumflow/1p/model.hh>
#include "spatialparams.hh"
namespace Dumux {
// forward declarations
template<class TypeTag> class OnePBulkProblem;
namespace Properties {
// create the type tag nodes
NEW_TYPE_TAG(OnePBulk, INHERITS_FROM(OneP));
NEW_TYPE_TAG(OnePBulkTpfa, INHERITS_FROM(OnePBulk, CCTpfaFacetCouplingModel));
NEW_TYPE_TAG(OnePBulkBox, INHERITS_FROM(OnePBulk, BoxFacetCouplingModel));
// Set the grid type
SET_TYPE_PROP(OnePBulk, Grid, Dune::ALUGrid<2, 2, Dune::cube, Dune::nonconforming>);
// Set the problem type
SET_TYPE_PROP(OnePBulk, Problem, OnePBulkProblem<TypeTag>);
// set the spatial params
SET_TYPE_PROP(OnePBulk, SpatialParams, OnePSpatialParams< typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
typename GET_PROP_TYPE(TypeTag, Scalar) >);
// the fluid system
SET_PROP(OnePBulk, FluidSystem)
{
private:
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
public:
using type = FluidSystems::OnePLiquid< Scalar, Components::Constant<1, Scalar> >;
};
} // end namespace Properties
/*!
* \ingroup OnePTests
* \brief Test problem for the incompressible one-phase model
* with coupling across the bulk grid facets
*/
template<class TypeTag>
class OnePBulkProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using PrimaryVariables = typename GridVariables::PrimaryVariables;
using Scalar = typename GridVariables::Scalar;
using FVGridGeometry = typename GridVariables::GridGeometry;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVGridGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVGridGeometry::SubControlVolumeFace;
using GridView = typename FVGridGeometry::GridView;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
public:
OnePBulkProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
std::shared_ptr<typename ParentType::SpatialParams> spatialParams,
const std::string& paramGroup = "")
: ParentType(fvGridGeometry, spatialParams, paramGroup)
, lowDimPermeability_(getParam<Scalar>("LowDim.SpatialParams.Permeability"))
, aperture_(getParam<Scalar>("Problem.FractureAperture"))
{}
//! Specifies the type of boundary condition at a given position
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
values.setAllDirichlet();
return values;
}
//! Specifies the type of interior boundary condition at a given position
BoundaryTypes interiorBoundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const
{
BoundaryTypes values;
values.setAllNeumann();
return values;
}
//! Evaluate the source term at a given position
NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
{
using std::cos;
using std::cosh;
Scalar u = (1.0 - lowDimPermeability_)*cos(globalPos[0])*cosh(aperture_/2);
return NumEqVector(u);
}
//! evaluates the exact solution at a given position.
Scalar exact(const GlobalPosition& globalPos) const
{
using std::cos;
using std::cosh;
const auto x = globalPos[0];
const auto y = globalPos[1];
return lowDimPermeability_*cos(x)*cosh(y) + (1.0 - lowDimPermeability_)*cos(x)*cosh(aperture_/2);
}
//! evaluates the exact gradient at a given position.
GlobalPosition exactGradient(const GlobalPosition& globalPos) const
{
using std::cos;
using std::sin;
using std::cosh;
using std::sinh;
const auto x = globalPos[0];
const auto y = globalPos[1];
GlobalPosition gradU;
gradU[0] = -lowDimPermeability_*sin(x)*cosh(y) + (lowDimPermeability_ - 1.0)*sin(x)*cosh(aperture_/2);
gradU[1] = lowDimPermeability_*cos(x)*sinh(y);
return gradU;
}
//! evaluates the Dirichlet boundary condition for a given position
PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{ return PrimaryVariables(exact(globalPos)); }
//! evaluates the Neumann boundary condition for a boundary segment
template<class ElementVolumeVariables>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf) const
{
auto pos = scvf.ipGlobal();
const Scalar k = this->spatialParams().permeabilityAtPos(pos);
const auto gradU = exactGradient(pos);
return NumEqVector( -1.0*k*(gradU*scvf.unitOuterNormal()) );
}
//! evaluate the initial conditions
PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{ return PrimaryVariables(1.0); }
//! returns the temperature in \f$\mathrm{[K]}\f$ in the domain
Scalar temperature() const
{ return 283.15; /*10°*/ }
//! sets the pointer to the coupling manager.
void setCouplingManager(std::shared_ptr<CouplingManager> cm)
{ couplingManagerPtr_ = cm; }
//! returns reference to the coupling manager.
const CouplingManager& couplingManager() const
{ return *couplingManagerPtr_; }
private:
std::shared_ptr<CouplingManager> couplingManagerPtr_;
Scalar lowDimPermeability_;
Scalar aperture_;
};
} // end namespace Dumux
#endif
#!/usr/bin/env python
from math import *
import subprocess
import sys
import numpy as np
if len(sys.argv) > 5 or len(sys.argv) < 3:
sys.stderr.write("Invalid number of arguments given. Please provide the following arguments (in this order):\n\
- the name of the executable (either test_facetcoupling_tpfa_1p1p or test_facetcoupling_box_1p1p)\n\
- between one and three fracture permeabilities you want to consider")
sys.exit(1)
# name of the executable
execName = sys.argv[1]
# array of permeabilities
k = []
for i in range(2, len(sys.argv)):
k.append(sys.argv[i])
# loop over the different permeability ratios
for permIndex in range(0, len(k)):
# remove the old log file
subprocess.call(['rm', execName + '.log'])
print("Removed old log file (" + execName + '.log' + ")!")
# for each given number of cells, create new .geo file, mesh & run simulation
numCells = np.linspace(10, 100, 10)
for cells in numCells:
geoFile = open("grids/hybridgrid.geo", 'r')
lines = geoFile.readlines()
firstLine = lines[0]
firstLine = firstLine.split(" ")
if firstLine[0] != "numElemsPerSide":
sys.stderr.write("Geo file is expected to have - numElemsPerSide = xx - in the first line!\n")
sys.exit(1)
# write a temporary geo-file using num elements
tmpGeoFile = open("grids/tmp.geo", 'w')
lines[0] = "numElemsPerSide = " + str(int(cells)) + ";\n"
for line in lines:
tmpGeoFile.write(line)
geoFile.close()
tmpGeoFile.close()
subprocess.call(['gmsh', '-2', 'grids/tmp.geo'])
subprocess.call(['./' + execName, 'facetcoupling_1p1p.input',
'-Grid.File', 'grids/tmp.msh',
'-Grid.NumElemsPerSide', str(int(cells)),
'-LowDim.SpatialParams.Permeability', str(k[permIndex]),
'-Problem.OutputFileName', execName + '.log'])
# remove geo and msh file
subprocess.call(['rm', 'grids/tmp.geo'])
subprocess.call(['rm', 'grids/tmp.msh'])
# check the rates and append them to the log file
logfile = open(execName + '.log', "r+")
eps = []
l2_matrix = []
l2_fracture = []
for line in logfile:
line = line.strip("\n")
line = line.split(",")
eps.append(float(line[0]))
l2_matrix.append(float(line[1]))
l2_fracture.append(float(line[3]))
matrix_rates = []
fracture_rates = []
logfile.truncate(0)
logfile.write("Matrix domain:\n")
logfile.write("-"*50 + "\n")
logfile.write("n\ta/L\t\terror_p\t\trate_p\n")
logfile.write("-"*50 + "\n")
for i in range(len(l2_matrix)-1):
if isnan(l2_matrix[i]) or isinf(l2_matrix[i]):
sys.stderr.write("l2 error norm is not a number!\n")
sys.exit(1)
if not (l2_matrix[i] < 1e-12 or l2_matrix[i] < 1e-12):
rateP = (log(l2_matrix[i])-log(l2_matrix[i+1]))/(log(eps[i])-log(eps[i+1]))
message = "{}\t{:0.4e}\t{:0.4e}\t{:0.4e}\n".format(i, eps[i], l2_matrix[i], rateP)
logfile.write(message)
matrix_rates.append(rateP)
else:
sys.stderr.write("error < 1e-12 does not seem reasonable!\n")
sys.exit(1)
# test is failed if the rate in the matrix is lower than 1.9
if abs(matrix_rates[len(matrix_rates)-1]) < 1.9:
sys.stderr.write("\n\nComputed a convergence rate of " + str(abs(matrix_rates[len(matrix_rates)-1])) + ", which is below the treshold of 1.9. Test fails...\n")
sys.exit(1)
# write last output for the matrix
i = len(l2_matrix)-1
message = "{}\t{:0.4e}\t{:0.4e}\n".format(i, eps[i], l2_matrix[i])
logfile.write(message)
logfile.write("\n\nFracture domain:\n")
logfile.write("-"*50 + "\n")
logfile.write("n\ta/L\t\terror_p\t\trate_p\n")
logfile.write("-"*50 + "\n")
for i in range(len(l2_fracture)-1):
if isnan(l2_fracture[i]) or isinf(l2_fracture[i]):
sys.stderr.write("l2 error norm is not a number!\n")
sys.exit(1)
if not (l2_fracture[i] < 1e-12 or l2_fracture[i] < 1e-12):
rateP = (log(l2_fracture[i])-log(l2_fracture[i+1]))/(log(eps[i])-log(eps[i+1]))
message = "{}\t{:0.4e}\t{:0.4e}\t{:0.4e}\n".format(i, eps[i], l2_fracture[i], rateP)
logfile.write(message)
fracture_rates.append(rateP)
else:
sys.stderr.write("error < 1e-12 does not seem reasonable!\n")
sys.exit(1)
# write last output for the matrix
i = len(l2_fracture)-1
message = "{}\t{:0.4e}\t{:0.4e}\n".format(i, eps[i], l2_fracture[i])
logfile.write(message)
logfile.close()
print("\nComputed the following convergence rates for:")
subprocess.call(['cat', execName + '.log'])
[Problem]
UseInteriorDirichletBCs = false
EnableGravity = false
FractureAperture = 1e-3
[Grid]
NumElemsPerSide = 25
File = ./grids/hybridgrid.msh
[Bulk]
Problem.Name = 1p_1p_bulk
SpatialParams.Permeability = 1
[LowDim]
SpatialParams.Permeability = 1e4
Problem.Name = 1p_1p_lowdim
[L2Error]
QuadratureOrder = 1
[Assembly]
NumericDifference.BaseEpsilon = 1e10
[Output]
EnableVTK = false # do not write .vtk files per default
numElemsPerSide = 25;
dx = 1/numElemsPerSide;
Point(1) = {0.5, -0.5, 0, dx};
Point(2) = {0.5, 0, 0, dx};
Point(3) = {0.5, 0.5, 0, dx};
Point(4) = {-0.5, 0.5, 0, dx};
Point(5) = {-0.5, 0, 0, dx};
Point(6) = {-0.5, -0.5, 0, dx};
Line(1) = {4, 5};
Transfinite Line{1} = Round(numElemsPerSide/2) + 1;
Line(2) = {5, 2};
Transfinite Line{2} = numElemsPerSide + 1;
Line(3) = {2, 3};
Transfinite Line{3} = Round(numElemsPerSide/2) + 1;
Line(4) = {3, 4};
Transfinite Line{4} = numElemsPerSide + 1;
Line(5) = {5, 6};
Transfinite Line{5} = Round(numElemsPerSide/2) + 1;
Line(6) = {6, 1};
Transfinite Line{6} = numElemsPerSide + 1;
Line(7) = {1, 2};
Transfinite Line{7} = Round(numElemsPerSide/2) + 1;
// we want line 2 to be explicitly meshed
Physical Line(1) = {2};
// the two surfaces
Line Loop(8) = {2, -7, -6, -5};
Plane Surface(9) = {8};
Line Loop(10) = {2, 3, 4, 1};
Plane Surface(11) = {10};
Transfinite Surface{9} = {1, 2, 5, 6};
Transfinite Surface{11} = {2, 3, 4, 5};
Recombine Surface{9};
Recombine Surface{11};
Physical Surface(1) = {9, 11};
This diff is collapsed.
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
121
1 0.5 -0.5 0
2 0.5 0 0
3 0.5 0.5 0
4 -0.5 0.5 0
5 -0.5 0 0
6 -0.5 -0.5 0
7 -0.5 0.4 0
8 -0.5 0.3000000000006937 0
9 -0.5 0.2000000000008325 0
10 -0.5 0.1000000000004163 0
11 -0.4000000000002772 0 0
12 -0.3000000000005547 0 0
13 -0.2000000000008318 0 0
14 -0.1000000000011089 0 0
15 -1.375122238300719e-12 0 0
16 0.09999999999889131 0 0
17 0.1999999999991684 0 0
18 0.2999999999994456 0 0
19 0.3999999999997228 0 0
20 0.5 0.09999999999977846 0
21 0.5 0.1999999999994866 0
22 0.5 0.2999999999994725 0
23 0.5 0.3999999999997362 0
24 0.4000000000002772 0.5 0
25 0.3000000000005547 0.5 0
26 0.2000000000008318 0.5 0
27 0.1000000000011089 0.5 0
28 1.375122238300719e-12 0.5 0
29 -0.09999999999889131 0.5 0
30 -0.1999999999991684 0.5 0
31 -0.2999999999994456 0.5 0
32 -0.3999999999997228 0.5 0
33 -0.5 -0.09999999999977846 0
34 -0.5 -0.1999999999994866 0
35 -0.5 -0.2999999999994725 0
36 -0.5 -0.3999999999997362 0
37 -0.4000000000002772 -0.5 0
38 -0.3000000000005547 -0.5 0
39 -0.2000000000008318 -0.5 0
40 -0.1000000000011089 -0.5 0
41 -1.375122238300719e-12 -0.5 0
42 0.09999999999889131 -0.5 0
43 0.1999999999991684 -0.5 0
44 0.2999999999994456 -0.5 0
45 0.3999999999997228 -0.5 0
46 0.5 -0.4 0
47 0.5 -0.3000000000006937 0
48 0.5 -0.2000000000008325 0
49 0.5 -0.1000000000004163 0
50 0.3999999999997228 -0.3999999999999737 0
51 0.2999999999994456 -0.3999999999999473 0
52 0.1999999999991684 -0.3999999999999208 0
53 0.09999999999889134 -0.3999999999998946 0
54 -1.375122238300719e-12 -0.3999999999998681 0
55 -0.1000000000011089 -0.3999999999998417 0
56 -0.2000000000008318 -0.3999999999998154 0
57 -0.3000000000005547 -0.3999999999997891 0
58 -0.4000000000002772 -0.3999999999997627 0
59 0.3999999999997228 -0.3000000000005716 0
60 0.2999999999994456 -0.3000000000004495 0
61 0.1999999999991684 -0.3000000000003273 0
62 0.09999999999889131 -0.3000000000002052 0
63 -1.375122238300719e-12 -0.3000000000000831 0
64 -0.1000000000011089 -0.299999999999961 0
65 -0.2000000000008318 -0.299999999999839 0
66 -0.3000000000005547 -0.2999999999997167 0
67 -0.4000000000002772 -0.2999999999995946 0
68 0.3999999999997228 -0.200000000000698 0
69 0.2999999999994456 -0.2000000000005633 0
70 0.1999999999991684 -0.2000000000004287 0
71 0.09999999999889131 -0.2000000000002942 0
72 -1.375122238300719e-12 -0.2000000000001595 0
73 -0.1000000000011089 -0.2000000000000249 0
74 -0.2000000000008318 -0.1999999999998904 0
75 -0.3000000000005547 -0.1999999999997557 0
76 -0.4000000000002772 -0.1999999999996211 0
77 0.3999999999997228 -0.1000000000003525 0
78 0.2999999999994456 -0.1000000000002887 0
79 0.1999999999991684 -0.1000000000002249 0
80 0.09999999999889131 -0.1000000000001611 0
81 -1.375108360512911e-12 -0.1000000000000973 0
82 -0.1000000000011089 -0.1000000000000336 0
83 -0.2000000000008319 -0.09999999999996981 0
84 -0.3000000000005547 -0.099999999999906 0
85 -0.4000000000002772 -0.09999999999984224 0
86 0.3999999999998338 0.09999999999984224 0
87 0.2999999999996674 0.09999999999990603 0
88 0.1999999999995011 0.09999999999996978 0
89 0.09999999999933482 0.1000000000000336 0
90 -8.250899963258007e-13 0.1000000000000974 0
91 -0.1000000000006654 0.1000000000001612 0
92 -0.2000000000004992 0.1000000000002249 0
93 -0.3000000000003329 0.1000000000002887 0
94 -0.4000000000001663 0.1000000000003525 0
95 0.3999999999999445 0.1999999999996212 0
96 0.2999999999998891 0.1999999999997557 0
97 0.1999999999998337 0.1999999999998903 0
98 0.09999999999977838 0.2000000000000249 0
99 -2.750299987752669e-13 0.2000000000001595 0
100 -0.1000000000002219 0.2000000000002941 0
101 -0.2000000000001665 0.2000000000004288 0
102 -0.3000000000001111 0.2000000000005633 0
103 -0.4000000000000554 0.2000000000006979 0
104 0.4000000000000555 0.2999999999995946 0
105 0.3000000000001111 0.2999999999997167 0
106 0.2000000000001664 0.2999999999998389 0
107 0.1000000000002219 0.299999999999961 0
108 2.750022431996513e-13 0.3000000000000831 0
109 -0.09999999999977835 0.3000000000002053 0
110 -0.1999999999998338 0.3000000000003273 0
111 -0.2999999999998892 0.3000000000004495 0
112 -0.3999999999999446 0.3000000000005715 0
113 0.4000000000001663 0.3999999999997627 0
114 0.3000000000003328 0.3999999999997889 0
115 0.2000000000004992 0.3999999999998153 0
116 0.1000000000006654 0.3999999999998417 0
117 8.250483629623773e-13 0.3999999999998681 0
118 -0.09999999999933486 0.3999999999998946 0
119 -0.1999999999995011 0.399999999999921 0
120 -0.2999999999996674 0.3999999999999472 0
121 -0.3999999999998337 0.3999999999999736 0
$EndNodes
$Elements
110
1 1 2 1 4 3 24
2 1 2 1 4 24 25
3 1 2 1 4 25 26
4 1 2 1 4 26 27
5 1 2 1 4 27 28
6 1 2 1 4 28 29
7 1 2 1 4 29 30
8 1 2 1 4 30 31
9 1 2 1 4 31 32
10 1 2 1 4 32 4
11 3 2 1 9 1 45 50 46
12 3 2 1 9 45 44 51 50
13 3 2 1 9 44 43 52 51
14 3 2 1 9 43 42 53 52
15 3 2 1 9 42 41 54 53
16 3 2 1 9 41 40 55 54
17 3 2 1 9 40 39 56 55
18 3 2 1 9 39 38 57 56
19 3 2 1 9 38 37 58 57
20 3 2 1 9 37 6 36 58
21 3 2 1 9 46 50 59 47
22 3 2 1 9 50 51 60 59
23 3 2 1 9 51 52 61 60
24 3 2 1 9 52 53 62 61
25 3 2 1 9 53 54 63 62
26 3 2 1 9 54 55 64 63
27 3 2 1 9 55 56 65 64
28 3 2 1 9 56 57 66 65