Commit 5cea017d authored by Timo Koch's avatar Timo Koch
Browse files

[test][md] Adjust test to reproduce Figure from reference paper

parent 41eeda7d
......@@ -79,4 +79,4 @@ dumux_add_test(NAME test_md_embedded_1d3d_1p1p_tpfatpfa_kernel
CMAKE_GUARD dune-foamgrid_FOUND
CMD_ARGS params.input -Vtk.OutputName test_md_embedded_1d3d_1p1p_tpfatpfa_kernel)
dune_symlink_to_source_files(FILES "params.input")
dune_symlink_to_source_files(FILES "params.input" "convergence.py")
#!/usr/bin/env python3
"""
# Script to reproduce Figure 5 of
# Koch et al (2020) JCP
# https://doi.org/10.1016/j.jcp.2020.109370
#
# The scipt also prints a table with errors and
# convergence rates
"""
import numpy as np
import subprocess
import sys
import multiprocessing as mp
import os
import matplotlib
import matplotlib.pyplot as plt
plt.style.use('ggplot')
font = {'family' : 'sans-serif',
'weight' : 'normal',
'size' : 8}
matplotlib.rc('font', **font)
def make_program_call(executable, args):
call = ['./' + executable]
for k, v in args.items():
call += ['-' + k, v]
return call
# remove the old log file and run
def run_simulation(params, filter=[]):
print("Running ", params["exec"], params["Vessel.Grid.Cells"], "...")
if params["exec"] not in filter:
with open(os.devnull, 'w') as devnull:
subprocess.run(['rm', '-f', params["Problem.OutputFilename"]], stdout=devnull)
subprocess.run(make_program_call(params["exec"], params))
executables = ["test_md_embedded_1d3d_1p1p_tpfatpfa_average", "test_md_embedded_1d3d_1p1p_tpfatpfa_surface", "test_md_embedded_1d3d_1p1p_tpfatpfa_kernel"]
method = {"test_md_embedded_1d3d_1p1p_tpfatpfa_average":1, "test_md_embedded_1d3d_1p1p_tpfatpfa_surface":3, "test_md_embedded_1d3d_1p1p_tpfatpfa_kernel":5}
label = {"test_md_embedded_1d3d_1p1p_tpfatpfa_average":"LS", "test_md_embedded_1d3d_1p1p_tpfatpfa_surface":"CSS", "test_md_embedded_1d3d_1p1p_tpfatpfa_kernel":"DS"}
marker = {"test_md_embedded_1d3d_1p1p_tpfatpfa_average":"o", "test_md_embedded_1d3d_1p1p_tpfatpfa_surface":"^", "test_md_embedded_1d3d_1p1p_tpfatpfa_kernel":"d"}
cells = [10, 20, 40, 80]
radius = 0.1
runs = []
res = {}
for e in executables:
for c in cells:
runs.append({"Vessel.Grid.Cells":str(c//2), "Tissue.Grid.Cells":str(c)+" "+str(c)+" "+str(c//2), "SpatialParams.Radius":str(radius), "Problem.OutputFilename":e+str(c)+".log", "exec":e});
res[e] = {}
for run in runs:
run_simulation(run)
name, numcells, result = (run["exec"], run["Vessel.Grid.Cells"], np.genfromtxt(run["Problem.OutputFilename"]))
print("Read file: ", run["Problem.OutputFilename"])
res[name][int(numcells)] = result
# produce output suitable for latex table
# Table 1
table1 = [["", "", "", "", "", "", ""] for i in range(len(cells)+1)]
table2 = [["", "", "", "", "", "", ""] for i in range(len(cells)+1)]
table3 = [["", "", "", "", "", "", ""] for i in range(len(cells)+1)]
dpi = 300.0
fig, axes = plt.subplots(1, 3, dpi=dpi, figsize=(8, 4))
for exec, result in res.items():
p3d = []
p1d = []
q = []
h = []
for cells, norms in sorted(result.items()):
h.append(norms[1])
p3d.append(norms[2])
p1d.append(norms[3])
q.append(norms[4])
rates3d = [(np.log10(p3d[i+1])-np.log10(p3d[i]))/(np.log10(h[i+1])-np.log10(h[i])) for i in range(len(p3d)-1)]
rates1d = [(np.log10(p1d[i+1])-np.log10(p1d[i]))/(np.log10(h[i+1])-np.log10(h[i])) for i in range(len(p1d)-1)]
ratesq = [(np.log10(q[i+1])-np.log10(q[i]))/(np.log10(h[i+1])-np.log10(h[i])) for i in range(len(q)-1)]
hR = np.array(h)/radius
axes[0].plot(hR, p3d, "--" + marker[exec], label=label[exec])
axes[0].set_title(r"$\vert\vert p^\mathbb{M}_{t,e} -p^\mathbb{M}_t \vert\vert_2$")
axes[1].plot(hR, p1d, "--" + marker[exec], label=label[exec])
axes[1].set_title(r"$\vert\vert p_{v,e} - p_v \vert\vert_2$")
axes[2].plot(hR, q, "--" + marker[exec], label=label[exec])
axes[2].set_title(r"$\vert\vert q_e - q \vert\vert_2$")
ofs = method[exec]
for i in range(len(h)):
table1[i][0] = "{:.4f}".format(h[i])
table2[i][0] = "{:.4f}".format(h[i])
table3[i][0] = "{:.4f}".format(h[i])
table1[i][ofs] = "{:.4e}".format(p3d[i])
table2[i][ofs] = "{:.4e}".format(p1d[i])
table3[i][ofs] = "{:.4e}".format(q[i])
if i > 0:
table1[i][ofs+1] = "{:.4e}".format(rates3d[i-1])
table2[i][ofs+1] = "{:.4e}".format(rates1d[i-1])
table3[i][ofs+1] = "{:.4e}".format(ratesq[i-1])
table1[len(h)][ofs] = "mean rate"
table2[len(h)][ofs] = "mean rate"
table3[len(h)][ofs] = "mean rate"
table1[len(h)][ofs+1] = "{:.4e}".format(np.mean(rates3d[1:]))
table2[len(h)][ofs+1] = "{:.4e}".format(np.mean(rates1d[1:]))
table3[len(h)][ofs+1] = "{:.4e}".format(np.mean(ratesq[1:]))
def print_table(table):
for row in table:
print(" & ".join(row) + "\\\\")
x = np.linspace(np.min(hR), np.max(hR), 10)
axes[0].plot(x, np.power(x*radius*0.4, 1.5), "--k", label=r"$\Delta$ 1.5")
axes[0].plot(x, np.power(x*radius*0.3, 2), "-.k", label=r"$\Delta$ 2")
axes[1].plot(x, np.power(x*radius*0.4, 1.5), "--k", label=r"$\Delta$ 1.5")
axes[1].plot(x, np.power(x*radius*0.3, 2), "-.k", label=r"$\Delta$ 2")
axes[2].plot(x, np.power(x*radius*2.5, 2), "--k", label=r"$\Delta$ 2")
axes[2].plot(x, np.power(x*radius*1.2, 2.5), "-.k", label=r"$\Delta$ 2.5")
for ax in axes:
ax.set_xscale("log")
ax.set_xlabel("$h/r_v$")
ax.set_xlim([np.max(hR), np.min(hR)])
ax.xaxis.set_minor_formatter(plt.NullFormatter())
ax.set_yscale("log")
ax.legend()
fig.tight_layout(rect=[0.03, 0.07, 1, 0.93], pad=0.4, w_pad=2.0, h_pad=1.0)
fig.savefig("rates_r01.pdf")
print("p3d")
print_table(table1)
print("p1d")
print_table(table2)
print("q")
print_table(table3)
......@@ -49,16 +49,9 @@ struct L2Norm
for (const auto& element : elements(gg.gridView()))
{
const auto geometry = element.geometry();
const auto center = geometry.center();
const auto elemSol = elementSolution(element, sol, gg);
// maybe exclude some elements from the norm
static const bool excludeInnerBulk = getParam<bool>("Problem.NormExcludeInnerBulk");
static const Scalar radius = getParam<Scalar>("SpatialParams.Radius") + 0.01;
using GridView = std::decay_t<decltype(gg.gridView())>;
if (int(GridView::dimension) == 3 && excludeInnerBulk && std::sqrt(center[0]*center[0] + center[1]*center[1]) < radius)
continue;
const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(geometry.type(), order);
for(auto&& qp : quad)
{
......@@ -83,25 +76,14 @@ struct L2Norm
for (const auto& element : elements(gg.gridView()))
{
const auto geometry = element.geometry();
const auto center = geometry.center();
// maybe exclude some elements from the norm
static const bool excludeInnerBulk = getParam<bool>("Problem.NormExcludeInnerBulk");
static const Scalar radius = getParam<Scalar>("SpatialParams.Radius") + 0.01;
using GridView = std::decay_t<decltype(gg.gridView())>;
if (int(GridView::dimension) == 3 && excludeInnerBulk && std::sqrt(center[0]*center[0] + center[1]*center[1]) < radius)
continue;
const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(geometry.type(), order);
for(auto&& qp : quad)
{
const auto globalPos = geometry.global(qp.position());
const Scalar pe = problem.exactSolution(globalPos);
norm += pe*pe*qp.weight()*geometry.integrationElement(qp.position());
}
for (auto&& qp : quad)
norm += 1.0*qp.weight()*geometry.integrationElement(qp.position());
}
return std::sqrt(norm);
return norm;
}
};
......
......@@ -25,13 +25,11 @@
#include <config.h>
#include <ctime>
#include <iostream>
#include <fstream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/istl/io.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
......@@ -39,9 +37,9 @@
#include <dumux/geometry/diameter.hh>
#include <dumux/linear/seqsolverbackend.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/grid/gridmanager_yasp.hh>
#include <dumux/io/grid/gridmanager_foam.hh>
#include <dumux/multidomain/traits.hh>
#include <dumux/multidomain/fvassembler.hh>
......@@ -178,9 +176,6 @@ int main(int argc, char** argv)
lowDimProblem->addVtkOutputFields(lowDimVtkWriter);
lowDimVtkWriter.write(0.0);
// an output file for the L2-norm
std::ofstream outFile_;
// compute hmax of both domains
double hMaxBulk = 0.0;
for (const auto& element : elements(bulkGridView))
......@@ -238,7 +233,7 @@ int main(int argc, char** argv)
// output the source terms
bulkProblem->computeSourceIntegral(sol[bulkIdx], *bulkGridVariables);
lowDimProblem->computeSourceIntegral(sol[lowDimIdx], *lowDimGridVariables);
auto sourceNorm = lowDimProblem->computeSourceIntegral(sol[lowDimIdx], *lowDimGridVariables);
// write vtk output
bulkVtkWriter.write(1.0);
......@@ -253,14 +248,16 @@ int main(int argc, char** argv)
// ouput result to terminal
std::cout << "-----------------------------------------------------\n";
std::cout << " L2_1d: " << norm1D << " hmax_1d: " << hMaxLowDim << '\n'
<< " L2_3d: " << norm3D << " hmax_3d: " << hMaxBulk << '\n';
std::cout << " L2_p1d: " << norm1D << " hmax_1d: " << hMaxLowDim << '\n'
<< " L2_p3d: " << norm3D << " hmax_3d: " << hMaxBulk << '\n'
<< " L2_q : " << sourceNorm << " hmax_3d: " << hMaxBulk << '\n';
std::cout << "-----------------------------------------------------" << std::endl;
// ... and file.
outFile_.open("1p_1p_norm.log", std::ios::app);
outFile_ << hMaxBulk << " " << hMaxLowDim << " " << norm3D << " " << norm1D << '\n';
outFile_.close();
{
std::ofstream outFile(getParam<std::string>("Problem.OutputFilename"));
outFile << hMaxBulk << " " << hMaxLowDim << " " << norm3D << " " << norm1D << " " << sourceNorm << '\n';
}
////////////////////////////////////////////////////////////
// finalize, print dumux message to say goodbye
......
[MixedDimension]
NumCircleSegments = 100
NumCircleSegments = 256
IntegrationOrder = 2
KernelWidthFactor = 1.0
WriteIntegrationPointsToFile = false
KernelIntegrationCRL = 0.01
[Vessel.Grid]
LowerLeft = 0 0 0
......@@ -20,9 +22,9 @@ Name = 1d
Name = 3d
[Problem]
NormIntegrationOrder = 2
NormExcludeInnerBulk = true
NormIntegrationOrder = 1
EnableGravity = false
OutputFilename = log.txt
[Vtk]
OutputName = test_md_embedded1d3d_1p1p
......
......@@ -294,9 +294,11 @@ public:
//! Called after every time step
//! Output the total global exchange term
void computeSourceIntegral(const SolutionVector& sol, const GridVariables& gridVars)
Scalar computeSourceIntegral(const SolutionVector& sol, const GridVariables& gridVars)
{
PrimaryVariables source(0.0);
SolutionVector source(sol.size());
SolutionVector sourceExact(sol.size());
Scalar volume = 0.0;
for (const auto& element : elements(this->gridGeometry().gridView()))
{
auto fvGeometry = localView(this->gridGeometry());
......@@ -309,11 +311,21 @@ public:
{
auto pointSources = this->scvPointSources(element, fvGeometry, elemVolVars, scv);
pointSources *= scv.volume()*elemVolVars[scv].extrusionFactor();
source += pointSources;
source[scv.dofIndex()] += pointSources;
auto a = scv.corner(0)[2];
auto b = scv.corner(1)[2];
if (a > b) std::swap(a, b);
sourceExact[scv.dofIndex()] += a*(1.0+0.5*a) - b*(1.0+0.5*b);
volume += scv.volume();
}
}
std::cout << "Global integrated source (1D): " << std::accumulate(source.begin(), source.end(), 0.0) << '\n';
std::cout << "Global integrated source exact (1D): " << std::accumulate(sourceExact.begin(), sourceExact.end(), 0.0) << '\n';
std::cout << "Global integrated source (1D): " << source << '\n';
source -= sourceExact;
const auto norm = source.two_norm()/sourceExact.two_norm();
std::cout << "relative L2 Norm source (1D): " << norm << '\n';
return source.two_norm()/volume;
}
/*!
......
......@@ -119,6 +119,7 @@ class TissueProblem : public PorousMediumFlowProblem<TypeTag>
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using GridView = typename GridGeometry::GridView;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
......@@ -184,7 +185,10 @@ public:
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
values.setAllDirichlet();
if (globalPos[2] > this->gridGeometry().bBoxMax()[2] - eps_ || globalPos[2] < this->gridGeometry().bBoxMin()[2] + eps_)
values.setAllNeumann();
else
values.setAllDirichlet();
return values;
}
......@@ -202,6 +206,49 @@ public:
return values;
}
/*!
* \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 scvf The sub control volume face
*
* Negative values mean influx.
* E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
*/
template<class ElementVolumeVariables, class ElemFluxVarsCache>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElemFluxVarsCache&,
const SubControlVolumeFace& scvf) const
{
NumEqVector flux(0.0);
// integrate over the scvf to compute the flux
const auto geometry = scvf.geometry();
Scalar derivative = 0.0;
const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension-1>::rule(geometry.type(), 4);
for(auto&& qp : quad)
{
auto globalPos = geometry.global(qp.position());
globalPos[2] = 0; // the derivative in z-direction is the exact solution evaluated with z=0
derivative += exactSolution(globalPos)*qp.weight()*geometry.integrationElement(qp.position());
}
const auto globalPos = scvf.ipGlobal();
if (globalPos[2] > this->gridGeometry().bBoxMax()[2] - eps_ )
flux[0] = -derivative/scvf.area();
else
flux[0] = derivative/scvf.area();
return flux;
}
// \}
/*!
......@@ -343,7 +390,7 @@ public:
//! The exact solution
Scalar exactSolution(const GlobalPosition &globalPos) const
{
const auto r = std::sqrt(globalPos[0]*globalPos[0] + globalPos[1]*globalPos[1]);
const auto r = std::hypot(globalPos[0], globalPos[1]);
static const auto R = getParam<Scalar>("SpatialParams.Radius");
if (CouplingManager::couplingMode == Embedded1d3dCouplingMode::kernel)
......
......@@ -4,8 +4,8 @@
<Piece NumberOfLines="20" NumberOfPoints="21">
<CellData Scalars="p">
<DataArray type="Float32" Name="p" NumberOfComponents="1" format="ascii">
1.02519 1.07551 1.12578 1.17601 1.22619 1.27633 1.32643 1.3765 1.42653 1.47653 1.52651 1.57645
1.62637 1.67626 1.72613 1.77598 1.8258 1.8756 1.92538 1.97513
1.02514 1.07537 1.12558 1.17574 1.22588 1.27599 1.32607 1.37612 1.42615 1.47615 1.52613 1.57609
1.62603 1.67595 1.72585 1.77574 1.8256 1.87545 1.92528 1.9751
</DataArray>
<DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
0 0 0 0 0 0 0 0 0 0 0 0
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment