Commit 7dea9881 authored by Timo Koch's avatar Timo Koch
Browse files

[richards][root] Add more output for comparison in benchmark study Schnepf et al

parent 5d1d423b
dune_symlink_to_source_files(FILES params.input convergencetest.py)
dune_symlink_to_source_files(FILES params.input convergencetest.py run_plot_m31.py)
dumux_add_test(NAME test_1p_rootbenchmark_tpfa
LABELS porousmediumflow 1p
......
......@@ -39,6 +39,7 @@
#if HAVE_GNUPLOT
#include <dumux/io/gnuplotinterface.hh>
#endif
#include <dumux/io/container.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
......@@ -106,6 +107,22 @@ int main(int argc, char** argv)
solver.solve(x);
vtkWriter.write(1.0);
std::vector<double> pos(leafGridView.size(0)), pressureHead(leafGridView.size(0));
for (const auto& element : elements(leafGridView))
{
const auto eIdx = leafGridView.indexSet().index(element);
pos[eIdx] = element.geometry().center()[GridGeometry::GridView::dimensionworld-1]*100;
pressureHead[eIdx] = 100*(x[eIdx][0] - 1e5)/(1000*9.81);
}
// Add Dirichlet BC
const auto pCollar = getParam<double>("Problem.RootCollarPressure");
const auto psiCollar = 100*(pCollar - 1e5)/(1000*9.81);
pos.push_back(0.0);
pressureHead.push_back(psiCollar);
writeContainerToFile(pos, "depth.txt", 15);
writeContainerToFile(pressureHead, "pressure_head.txt", 15);
// write more output for benchmark plot
problem->outputL2Norm(x);
......
#!/usr/bin/env python3
import sys
import numpy as np
import subprocess
# Run benchmark M3.1 root scenario
subprocess.run(["./test_1p_rootbenchmark_tpfa", "-Grid.Refinement", "5"])
try:
import matplotlib.pyplot as plt
except ImportError:
print("Matplotlib not found so no result plots will be created.")
sys.exit(0)
# Create Fig. of Schnepf et at 2020 benchmark
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(5, 5), dpi=72)
ax = axes
# set labels
ax.set_xlabel(r"root water pressure head (cm)")
ax.set_ylabel(r"depth (cm)")
depth = np.genfromtxt("depth.txt")
head = np.genfromtxt("pressure_head.txt")
ax.plot(head, depth)
textOutput = ""
textOutput += ",".join(f"{i:.10e}" for i in depth) + "\n"
textOutput += ",".join(f"{i:.10e}" for i in head) + "\n"
with open("root.txt", "w") as out:
out.write(textOutput)
fig.tight_layout()
fig.savefig("root.png")
......@@ -155,7 +155,7 @@ int main(int argc, char** argv)
}
};
const auto updateDeltaEtaOutput = [&](auto& deltaEtaNum, auto& thetaNum)
const auto updateDeltaEtaOutput = [&](auto& deltaEtaNum, auto& depthNum, auto& thetaNum)
{
const auto t = timeLoop->time()/86400.0;
auto fvGeometry = localView(*gridGeometry);
......@@ -164,6 +164,7 @@ int main(int argc, char** argv)
{
const auto eIdx = leafGridView.indexSet().index(element);
const auto pos = element.geometry().center()[GridGeometry::GridView::dimensionworld -1];
depthNum[eIdx] = pos*100;
deltaEtaNum[eIdx] = solutionInfiltration->deltaEta(pos*100, t);
fvGeometry.bindElement(element);
......@@ -180,7 +181,9 @@ int main(int argc, char** argv)
GnuplotInterface<double> gnuplotInfilt(openGnuplotWindow);
gnuplotInfilt.setOpenPlotWindow(openGnuplotWindow);
std::vector<double> deltaEtaNum(leafGridView.size(0)), thetaNum(leafGridView.size(0));
GnuplotInterface<double> gnuplotInfilt2(openGnuplotWindow);
gnuplotInfilt2.setOpenPlotWindow(openGnuplotWindow);
std::vector<double> deltaEtaNum(leafGridView.size(0)), thetaNum(leafGridView.size(0)), depthNum(leafGridView.size(0));
std::vector<double> analyticalWaterContent;
if (scenario == BenchmarkScenario::infiltration)
......@@ -230,7 +233,7 @@ int main(int argc, char** argv)
updateAnalyticalWaterContent(analyticalWaterContent);
vtkWriter->write(timeLoop->time());
updateDeltaEtaOutput(deltaEtaNum, thetaNum);
updateDeltaEtaOutput(deltaEtaNum, depthNum, thetaNum);
const auto posRef = interpolate<InterpolationPolicy::LinearTable>(solutionInfiltration->refWaterContent(), thetaNum, pos);
std::cout << Fmt::format("Simulated x_a: {}cm\n", posRef);
......@@ -241,6 +244,12 @@ int main(int argc, char** argv)
);
gnuplotInfilt.plot(Fmt::format("infiltration_{}_{:d}", soilType, checkPointCounter));
gnuplotInfilt2.addDataSetToPlot(thetaNum, depthNum,
Fmt::format("theta_depth_num_{}_{:d}.dat", soilType, checkPointCounter),
Fmt::format("w p t 'DuMux {:.1f} days'", timeLoop->time()/86400.0)
);
gnuplotInfilt2.plot(Fmt::format("infiltration_depth_{}_{:d}", soilType, checkPointCounter));
}
///////////////////////////////////////
......
......@@ -5,9 +5,9 @@ import numpy as np
import subprocess
# Run benchmark M2.1 infiltration scenarios
subprocess.run(["./test_richards_benchmark_tpfa", "params_infiltration_sand.input"])
subprocess.run(["./test_richards_benchmark_tpfa", "params_infiltration_loam.input"])
subprocess.run(["./test_richards_benchmark_tpfa", "params_infiltration_clay.input"])
#subprocess.run(["./test_richards_benchmark_tpfa", "params_infiltration_sand.input", "-Grid.Refinement", "2"])
#subprocess.run(["./test_richards_benchmark_tpfa", "params_infiltration_loam.input", "-Grid.Refinement", "2"])
#subprocess.run(["./test_richards_benchmark_tpfa", "params_infiltration_clay.input", "-Grid.Refinement", "2"])
try:
import matplotlib.pyplot as plt
......@@ -57,3 +57,42 @@ axes[2].legend()
fig.tight_layout()
plt.savefig("benchmark_infiltration.png")
# Create Fig. of Schnepf et al 2020
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 5), dpi=72, sharey=True)
textOutput = ""
# set labels
axes[0].set_title(r"$\theta$")
axes[0].set_ylabel(r"depth (cm)")
axes[0].set_ylim([-200.0, 0.0])
# first plot
for i, (t, marker) in enumerate(zip([0.1, 0.2, 0.3], ["x", "D", "o"])):
theta, depth = np.genfromtxt("theta_depth_num_sand_{}.dat".format(i+1)).T
textOutput += ",".join(f"{i:.10e}" for i in depth) + "\n"
textOutput += ",".join(f"{i:.10e}" for i in theta) + "\n"
axes[0].plot(theta, depth, linestyle='none', fillstyle='none', marker=marker, label="{:.1f} days".format(t))
axes[0].legend()
# second plot
for i, (t, marker) in enumerate(zip([0.2, 0.5, 1.0], ["x", "D", "o"])):
theta, depth = np.genfromtxt("theta_depth_num_loam_{}.dat".format(i+1)).T
textOutput += ",".join(f"{i:.10e}" for i in depth) + "\n"
textOutput += ",".join(f"{i:.10e}" for i in theta) + "\n"
axes[1].plot(theta, depth, linestyle='none', fillstyle='none', marker=marker, label="{:.1f} days".format(t))
axes[1].legend()
# third plot
for i, (t, marker) in enumerate(zip([0.1, 0.2, 0.5], ["x", "D", "o"])):
theta, depth = np.genfromtxt("theta_depth_num_clay_{}.dat".format(i+1)).T
textOutput += ",".join(f"{i:.10e}" for i in depth) + "\n"
textOutput += ",".join(f"{i:.10e}" for i in theta) + "\n"
axes[2].plot(theta, depth, linestyle='none', fillstyle='none', marker=marker, label="{:.1f} days".format(t))
axes[2].legend()
with open("benchmark_infiltration_depth.txt", "w") as out:
out.write(textOutput)
fig.tight_layout()
plt.savefig("benchmark_infiltration_depth.png")
......@@ -5,10 +5,10 @@ import numpy as np
import subprocess
# Run benchmark M2.2 evaporation scenarios
subprocess.run(["./test_richards_benchmark_tpfa", "params_evaporation_sand.input"])
subprocess.run(["./test_richards_benchmark_tpfa", "params_evaporation_loam1.input"])
subprocess.run(["./test_richards_benchmark_tpfa", "params_evaporation_loam2.input"])
subprocess.run(["./test_richards_benchmark_tpfa", "params_evaporation_clay.input"])
subprocess.run(["./test_richards_benchmark_tpfa", "params_evaporation_sand.input", "-Grid.Refinement", "2"])
subprocess.run(["./test_richards_benchmark_tpfa", "params_evaporation_loam1.input", "-Grid.Refinement", "2"])
subprocess.run(["./test_richards_benchmark_tpfa", "params_evaporation_loam2.input", "-Grid.Refinement", "2"])
subprocess.run(["./test_richards_benchmark_tpfa", "params_evaporation_clay.input", "-Grid.Refinement", "2"])
# Create equally spaced plots to improve test robustness.
# In the case of slightly different time steps this makes sure
......@@ -46,13 +46,19 @@ axes[0,0].set_ylim([0.0, 0.1])
axes[1,0].set_ylim([0.0, 0.3])
plots = {"sand": (0,0), "loam1": (0,1), "loam2": (1,0), "clay": (1,1)}
textOutput = ""
for soil, index in plots.items():
t, e = np.genfromtxt("rate_actual_{}.dat".format(soil)).T
t, e = t, e*0.1
textOutput += ",".join(f"{i:.10e}" for i in t) + "\n"
textOutput += ",".join(f"{i:.10e}" for i in e) + "\n"
axes[index].plot(t, e, "k", label="DuMux")
t, e = np.genfromtxt("rate_analytical_{}.dat".format(soil)).T
t, e = t, e*0.1
axes[index].plot(t, e, marker="x", linestyle='none', label="Benchmark")
axes[index].legend()
with open("benchmark_evaporation.txt", "w") as out:
out.write(textOutput)
plt.savefig("benchmark_evaporation.png")
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