templated reports, spatial / temporal convergence

parametrized jupytext templates can be converted to and run as jupyter
notebooks. These in turn can be used as is or rendered as HTML.
removed previous benchmark files, calculated on the fly with ogs6py
different markers for richardson extrapolation values

docs/conf.py

@@ -84,7 +84,11 @@ panels_add_bootstrap_css = False
 ### apidoc / autodoc settings ###
 apidoc_module_dir = "../ogstools"
 apidoc_output_dir = "reference"
-apidoc_excluded_paths = ["../**/examples/**", "./**/tests/**"]
+apidoc_excluded_paths = [
+    "../**/examples/**",
+    "./**/tests/**",
+    "../**/**/templates/**",
+]
 apidoc_separate_modules = True
 apidoc_module_first = True
 apidoc_extra_args = ["--force", "--implicit-namespaces"]

docs/examples/howto_studies/plot_convergence_study.py

-"""
-Convergence study
-=================
-
-This script performs a convergence study and generates plots to analyze the
-convergence of numerical simulations. It uses data from the following benchmark
-with multiple discretizations to evaluate the accuracy of the numerical
-solutions.
-https://www.opengeosys.org/docs/benchmarks/elliptic/elliptic-neumann/
-
-Here is some theoretical background for the Richardson extrapolation:
-https://www.grc.nasa.gov/www/wind/valid/tutorial/spatconv.html
-"""
-
-# %%
-# First import the dependencies and adjust some plot settings.
-
-import numpy as np
-
-import ogstools.meshplotlib as mpl
-from ogstools.propertylib import Scalar, Vector
-from ogstools.studies.convergence import (
-    convergence_metrics,
-    examples,
-    grid_convergence,
-    plot_convergence,
-    plot_convergence_errors,
-    richardson_extrapolation,
-)
-
-mpl.setup.reset()
-mpl.setup.show_element_edges = True
-mpl.setup.ax_aspect_ratio = 1
-
-# %%
-# Let's have a look at the different discretizations. The 3 finest will be used
-# for the Richadson extrapolation. The coarsest of those will be used for the
-# topology to evaluate the results. We inspect the primary variable,
-# in this case the hydraulic head.
-
-mesh_property = Scalar(
-    data_name="pressure",
-    data_unit="m",
-    output_unit="m",
-    output_name="hydraulic head",
-)
-fig = mpl.plot(np.reshape(examples.meshes, (2, 3)), mesh_property)
-
-# %%
-# Now we calculate the convergence ratio on the whole mesh. Values close to
-# 1 indicate that we are in asymptotic range of convergence. We see, that in the
-# bottom right corner, there is some small discrepancy, which is explained in
-# the benchmark documentation by "incompatible boundary conditions imposed on
-# the bottom right corner of the domain." In the end we will see the
-# implications for the secondary variable, the velocity.
-
-topology = examples.meshes[-3]
-conv_results = grid_convergence(examples.meshes, mesh_property, topology)
-fig = mpl.plot(conv_results, "grid_convergence")
-
-# %%
-# The Richardson extrapolation can be easily calculated. Again, it uses the 3
-# finest meshes from the given list of meshes.
-
-richardson = richardson_extrapolation(examples.meshes, mesh_property, topology)
-analytical = examples.analytical_solution(topology)
-fig = mpl.plot([richardson, analytical], mesh_property)
-fig.axes[0].set_title("Richardson extrapolation")
-fig.axes[1].set_title("Analytical Solution")
-fig.show()
-
-# %%
-# Now we can compute some convergence metrics and display them in a table, ...
-
-metrics = convergence_metrics(examples.meshes, richardson, mesh_property)
-metrics.style.format("{:,.5g}").hide()
-
-# %%
-# ... plot the converging values in absolute scale ...
-
-mpl.core.plt.rcdefaults()
-fig = plot_convergence(metrics, mesh_property)
-
-# %%
-# ... and the relative errors in loglog-scale. Note: since the Minimum doesn't
-# change in the different discretizations, the error is zero, thus there is no
-# curve for it in this plot.
-
-fig = plot_convergence_errors(metrics)
-
-
-# %%
-# Now let's inspect the velocity field. We see, that in the bottom right corner,
-# the velocity magnitude seems to be steadily increasing.
-
-mesh_property = Vector("v", "m/s", "m/s", "velocity")
-mpl.setup.num_streamline_interp_pts = None
-fig = mpl.plot(np.reshape(examples.meshes, (2, 3)), mesh_property)
-
-# %%
-# Looking at the grid convergence, we see the bottom left and right corners
-# deviating quite a bit from the desired value of 1. Thus we know, at these
-# points the mesh isn't properly converging (at least for the velocity field).
-
-conv_results = grid_convergence(examples.meshes, mesh_property, topology)
-fig = mpl.plot(conv_results, "grid_convergence")
-
-# %%
-# The Richardson extrapolation shows an anomalous value for the velocity in
-# the bottom right corner, hinting at a singularity there which is caused by
-# the incompatibility of the boundary conditions on the bottom and right in
-# this singular point. Regardsless of this, the benchmark gives a convergent
-# solution for the pressure field, but not for the velocity field.
-richardson = richardson_extrapolation(examples.meshes, mesh_property, topology)
-fig = mpl.plot(richardson, mesh_property)

docs/examples/howto_studies/plot_convergence_study_nuclear_decay.py

+"""
+Convergence study (temporal refinement)
+=======================================
+
+This example shows one possible implementation of how to do a convergence study
+with temporal refinement. For this, a simple model using a time dependent heat
+source on one side and constant temperature on the opposite side was set up.
+The heat source is generated with the `nuclearwasteheat` model.
+
+Here is some theoretical background for the topic of grid convergence:
+https://www.grc.nasa.gov/www/wind/valid/tutorial/spatconv.html
+
+At least three meshes from simulations of increasing temporal refinement are
+required for the convergence study. The topology has to stay the same.
+
+The below code cells will generate the simulation results and are evaluated for
+convergence at the end.
+
+First, the required packages are imported and an output directory is created:
+"""
+
+# %%
+from pathlib import Path
+from shutil import rmtree
+from tempfile import mkdtemp
+
+import matplotlib.pyplot as plt
+import numpy as np
+import vtuIO
+from IPython.display import HTML
+from ogs6py import ogs
+
+from ogstools import meshlib, msh2vtu, physics, propertylib, studies, workflow
+
+temp_dir = Path(mkdtemp(prefix="nuclear_decay"))
+
+# %% [markdown]
+# Let's Visualize the temporal evolution of the source term and it's
+# discretization in the simulations. We see, that with coarse time steps, the
+# applied heat will be overestimated at first and underestimated once the heat
+# curve has reached its maximum. The same is true for the resulting temperature.
+
+# %%
+n_refinements = 4
+t_end = 180
+sec_per_yr = 365.25 * 86400
+time = np.append(0, np.geomspace(1, t_end, num=100)) * sec_per_yr
+heat = physics.nuclearwasteheat.repo_2020_conservative.heat
+fig, ax = plt.subplots(figsize=(8, 4))
+# print(heat(time))
+ax.plot(time / sec_per_yr, heat(time) / 1e3, lw=2, label="reference", color="k")
+
+for r in range(n_refinements):
+    dt = 30.0 / (2.0**r)
+    time = np.linspace(0, t_end, int(t_end / dt) + 1) * sec_per_yr
+    edges = np.append(0, time) / sec_per_yr
+    ax.stairs(heat(time) / 1e3, edges, label=f"{dt=}", baseline=None, lw=1.5)
+ax.set_xlabel("time / yrs")
+ax.set_ylabel("heat / kW")
+ax.legend()
+fig.show()
+
+# %% [markdown]
+# The mesh and its boundaries are generated easily via gmsh and msh2vtu:
+
+# %%
+msh_path = temp_dir / "square.msh"
+meshlib.rect_mesh(lengths=100, n_edge_cells=[10, 1], out_name=msh_path)
+_ = msh2vtu.msh2vtu(msh_path, output_path=temp_dir, log_level="ERROR")
+
+# %% [markdown]
+# Let's run the different simulations with increasingly fine temporal
+# discretization via ogs6py, extract the temperature evolution via vtuIO and
+# plot it:
+
+# %%
+r_range = range(n_refinements)
+results_tmax = [temp_dir / f"nuclear_decay_tmax_{r}.vtu" for r in r_range]
+results_qmax = [temp_dir / f"nuclear_decay_qmax_{r}.vtu" for r in r_range]
+fig, ax = plt.subplots(figsize=(8, 4))
+
+for r in range(n_refinements):
+    model = ogs.OGS(
+        PROJECT_FILE=temp_dir / "default.prj",
+        INPUT_FILE=studies.convergence.examples.nuclear_decay_prj,
+    )
+    dt = 30.0 / (2.0**r)
+    model.replace_text(str(dt * sec_per_yr), ".//delta_t")
+    model.write_input()
+    script_path = Path(studies.convergence.examples.nuclear_decay_bc).parent
+    ogs_args = f"-m {temp_dir} -o {temp_dir} -s {script_path}"
+    model.run_model(write_logs=False, args=ogs_args)
+
+    pvd_path = str(temp_dir / "nuclear_decay.pvd")
+    pvd = meshlib.MeshSeries(pvd_path)
+    result_qmax = pvd.read_closest(30 * sec_per_yr)
+    result_tmax = pvd.read_closest(150 * sec_per_yr)
+    result_qmax.add_field_data(dt, "timestep_size")
+    result_tmax.add_field_data(dt, "timestep_size")
+    result_qmax.save(results_qmax[r])
+    result_tmax.save(results_tmax[r])
+
+    pvdio = vtuIO.PVDIO(pvd_path, dim=2, interpolation_backend="vtk")
+    max_temperature = propertylib.presets.temperature(
+        pvdio.read_time_series("temperature", {"pt0": [0, 0, 0]})["pt0"]
+    )
+    ts = pvdio.timesteps / sec_per_yr
+    ax.plot(ts, max_temperature, lw=1.5, label=f"{dt=}")
+ax.set_xlabel("time / yrs")
+ax.set_ylabel("max T / °C")
+ax.legend()
+fig.show()
+
+# %% [markdown]
+# Temperature convergence at maximum heat production (t=30 yrs)
+# -------------------------------------------------------------
+#
+# The grid convergence at this timepoint deviates significantly from 1,
+# meaning the convergence is suboptimal (at least on the left boundary where the
+# heating happens). The chosen timesteps are still to coarse to reach an
+# asymptotic rate of convergence. The model behavior at this early part of the
+# simulation is still very dynamic and needs finer timesteps to be captured with
+# great accuracy. Nevertheless, the maximum temperature converges (sublinearly)
+# from overestimated values, as expected.
+
+# %%
+report_name = str(temp_dir / "report.ipynb")
+studies.convergence.run_convergence_study(
+    output_name=report_name,
+    mesh_paths=results_qmax,
+    topology_path=results_qmax[-3],
+    property_name="temperature",
+    refinement_ratio=2.0,
+)
+HTML(workflow.jupyter_to_html(report_name, show_input=False))
+
+# %% [markdown]
+# Temperature convergence at maximum temperature (t=150 yrs)
+# ----------------------------------------------------------
+#
+# The temperature convergence at this timepoint is much closer to 1, signifying
+# a good convergence behaviour. The temperature gradient at this point in time
+# is already settled and can be the solution convergences to a good degree with
+# the chosen timesteps. Despite the good grid convergence, the maximum
+# temperature converges only sublinearly, but as expected, from underestimated
+# values.
+
+# %%
+studies.convergence.run_convergence_study(
+    output_name=report_name,
+    mesh_paths=results_tmax,
+    topology_path=results_tmax[-3],
+    property_name="temperature",
+    refinement_ratio=2.0,
+)
+HTML(workflow.jupyter_to_html(report_name, show_input=False))
+
+# %%
+
+# sphinx_gallery_start_ignore
+
+# Removing the created files to keep the code repository clean for developers.
+# If you want to use the created jupyter notebook further, skip this step.
+rmtree(temp_dir)
+
+# sphinx_gallery_end_ignore

docs/examples/howto_studies/plot_convergence_study_steady_state_diffusion.py

+"""
+Convergence study (spatial refinement)
+======================================
+
+This example shows one possible implementation of how to do a convergence study.
+It uses data from the following benchmark with multiple discretizations to
+evaluate the accuracy of the numerical solutions.
+https://www.opengeosys.org/docs/benchmarks/elliptic/elliptic-neumann/
+
+Here is some theoretical background for the topic of grid convergence:
+https://www.grc.nasa.gov/www/wind/valid/tutorial/spatconv.html
+
+At least three meshes of increasing refinement are required for the convergence
+study. The three finest meshes are used to calculated the Richardson
+extrapolation. It is recommended to use the third coarsest mesh for the topology
+to evaluate the results. The finer meshes should share the same nodes of the
+chosen topology, otherwise interpolation will influence the results. With
+unstructured grids this can be achieved as well with refinement by splitting.
+
+First, the required packages are imported and an output directory is created:
+"""
+
+# %%
+from pathlib import Path
+from shutil import rmtree
+from tempfile import mkdtemp
+
+from IPython.display import HTML
+from ogs6py import ogs
+
+from ogstools import meshlib, msh2vtu, studies, workflow
+from ogstools.studies.convergence.examples import (
+    steady_state_diffusion_analytical_solution,
+)
+
+temp_dir = Path(mkdtemp())
+report_name = str(temp_dir / "report.ipynb")
+result_paths = []
+
+# %% [markdown]
+# The meshes and their boundaries are generated easily via gmsh and msh2vtu.
+# Then we run the different simulations with increasingly fine spatial
+# discretization via ogs6py and store the results for the convergence study.
+
+# %%
+for i in range(6):
+    msh_path = temp_dir / "square.msh"
+    meshlib.gmsh_meshing.rect_mesh(
+        n_edge_cells=2**i, structured_grid=True, out_name=msh_path
+    )
+    msh2vtu.msh2vtu(
+        input_filename=msh_path, output_path=temp_dir, log_level="ERROR"
+    )
+
+    model = ogs.OGS(
+        PROJECT_FILE=temp_dir / "default.prj",
+        INPUT_FILE=studies.convergence.examples.steady_state_diffusion_prj,
+    )
+    model.write_input()
+    ogs_args = f"-m {temp_dir} -o {temp_dir}"
+    model.run_model(write_logs=False, args=ogs_args)
+
+    result = meshlib.MeshSeries(
+        str(temp_dir / "steady_state_diffusion.pvd")
+    ).read(-1)
+    result_path = temp_dir / f"steady_state_diffusion_{i}.vtu"
+    result.save(result_path)
+    result_paths += [str(result_path)]
+
+# %% [markdown]
+# Here we choose the topology to evaluate the reults on and calculate the
+# analytical solution on it:
+
+# %%
+topology_path = result_paths[-3]
+analytical_solution_path = temp_dir / "analytical_solution.vtu"
+steady_state_diffusion_analytical_solution(topology_path).save(
+    analytical_solution_path
+)
+
+# %% [markdown]
+# Hydraulic pressure convergence
+# ------------------------------
+#
+# The pressure field of this model is converging well. The convergence ratio
+# is approximately 1 on the whole mesh and looking at the relative errors we
+# see a quadratic convergence behavior.
+
+# %%
+studies.convergence.run_convergence_study(
+    output_name=report_name,
+    mesh_paths=result_paths,
+    topology_path=topology_path,
+    property_name="hydraulic_height",
+    refinement_ratio=2.0,
+    reference_solution_path=str(analytical_solution_path),
+)
+HTML(workflow.jupyter_to_html(report_name, show_input=False))
+
+# %% [markdown]
+# Darcy velocity convergence
+# --------------------------
+#
+# For the velocity we some discrepancy of the convergence ratio in the bottom
+# right corner. Thus we know, at these points the mesh isn't properly
+# converging (at least for the velocity field).
+
+# We see, that in the bottom right corner, the velocity magnitude seems to be
+# steadily increasing, which is also reflected in the Richardson extrapolation,
+# which shows an anomalous high value in this spot, hinting at a singularity
+# there. This is explained by the notion in the benchmark's documentation of
+# "incompatible boundary conditions imposed on the bottom right corner of the
+# domain." Regardless of this, the benchmark gives a convergent solution for
+# the pressure field.
+
+# %%
+
+studies.convergence.run_convergence_study(
+    output_name=report_name,
+    mesh_paths=result_paths,
+    topology_path=topology_path,
+    refinement_ratio=2.0,
+    property_name="velocity",
+)
+HTML(workflow.jupyter_to_html(report_name, show_input=True))
+
+# %%
+
+# sphinx_gallery_start_ignore
+
+# Removing the created files to keep the code repository clean for developers.
+# If you want to use the created jupyter notebook further, skip this step.
+rmtree(temp_dir)
+
+# sphinx_gallery_end_ignore

ogstools/studies/convergence/__init__.py

 # Author: Florian Zill (Helmholtz Centre for Environmental Research GmbH - UFZ)
 """functions to generate a convergence study."""
 
+from . import examples
 from .convergence import (
+    add_grid_spacing,
     convergence_metrics,
     grid_convergence,
     log_fit,
@@ -9,14 +11,16 @@ from .convergence import (
     plot_convergence_errors,
     richardson_extrapolation,
 )
-from .generate_report import execute_convergence_study
+from .study import run_convergence_study
 
 __all__ = [
+    "add_grid_spacing",
     "convergence_metrics",
-    "execute_convergence_study",
+    "examples",
     "grid_convergence",
     "log_fit",
     "plot_convergence",
     "plot_convergence_errors",
     "richardson_extrapolation",
+    "run_convergence_study",
 ]

ogstools/studies/convergence/convergence.py

 from copy import deepcopy
-from typing import Optional
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -26,29 +25,15 @@ def add_grid_spacing(mesh: pv.DataSet) -> pv.DataSet:
     dim = mesh.get_cell(0).dimension
     key = ["Length", "Area", "Volume"][dim - 1]
     _mesh = mesh.compute_cell_sizes()
-    _mesh["grid_spacing"] = _mesh.cell_data[key] ** (1.0 / dim)
+    _mesh.cell_data["grid_spacing"] = _mesh.cell_data[key] ** (1.0 / dim)
     return _mesh
 
 
-def get_refinement_ratio(
-    meshes: list[pv.DataSet], topology: pv.DataSet
-) -> pv.DataSet:
-    _meshes = resample(topology, [add_grid_spacing(mesh) for mesh in meshes])
-    spacings = np.stack([mesh["grid_spacing"] for mesh in _meshes], axis=0)
-    refinement_ratios = spacings[:-1] / spacings[1:]
-    mean_ratios = np.mean(refinement_ratios, axis=0)
-    result = deepcopy(topology)
-    result.clear_point_data()
-    result.clear_cell_data()
-    result["refinement_ratio"] = mean_ratios
-    return result
-
-
 def grid_convergence(
     meshes: list[pv.DataSet],
     property: Property,
     topology: pv.DataSet,
-    refinement_ratio: Optional[float] = None,
+    refinement_ratio: float,
 ) -> pv.DataSet:
     """
     Calculate the grid convergence field for the given meshes on the topology.
@@ -73,16 +58,14 @@ def grid_convergence(
     f3 = cast(_meshes[-3].point_data[property.data_name])
     f2 = cast(_meshes[-2].point_data[property.data_name])
     f1 = cast(_meshes[-1].point_data[property.data_name])
-    if refinement_ratio is None:
-        r = get_refinement_ratio(meshes, topology)["refinement_ratio"]
-    else:
-        r = np.ones(f1.shape) * refinement_ratio
+    r = np.ones(f1.shape) * refinement_ratio
     a = f3 - f2
     b = f2 - f1
     zeros = np.zeros_like
     ones = np.ones_like
     c = np.divide(a, b, out=ones(a), where=(b != 0))
-    p = np.log(np.abs(c)) / np.log(r)
+    with np.errstate(divide="ignore"):
+        p = np.log(np.abs(c)) / np.log(r)
     rpm1 = r**p - 1
     _gci23 = np.divide(np.abs(a), f3, out=zeros(a), where=(f3 != 0.0))
     _gci12 = np.divide(np.abs(b), f2, out=zeros(a), where=(f2 != 0.0))
@@ -105,7 +88,7 @@ def richardson_extrapolation(
     meshes: list[pv.DataSet],
     property: Property,
     topology: pv.DataSet,
-    refinement_ratio: Optional[float] = None,
+    refinement_ratio: float,
 ) -> pv.DataSet:
     """
     Estimate a better approximation of a property on a mesh.
@@ -118,7 +101,7 @@ def richardson_extrapolation(
     :param meshes:              At least three meshes with constant refinement.
     :param property:            The property to be extrapolated.
     :param topology:            The topology on which the extrapolation is done.
-    :param refinement_ratio:    If not given, it is calculated automatically
+    :param refinement_ratio:    Refinement ratio (spatial or temporal).
 
     :returns:                   Richardson extrapolation of the given property.
     """
@@ -155,9 +138,13 @@ def convergence_metrics(
     def _data(m: pv.DataSet):
         return property.magnitude.strip_units(m.point_data[property.data_name])
 
-    el_lens = [
-        np.mean(add_grid_spacing(mesh)["grid_spacing"]) for mesh in meshes
-    ] + [0.0]
+    if np.all(["timestep_size" in mesh.field_data for mesh in meshes]):
+        x = [mesh.field_data["timestep_size"][0] for mesh in meshes]
+        x_str = "time step size"
+    else:
+        x = [np.mean(add_grid_spacing(mesh)["grid_spacing"]) for mesh in meshes]
+        x_str = "mean element length"
+    x += [0.0]
     _meshes = meshes + [reference]
     maxs = [np.max(_data(m)) for m in _meshes]
     mins = [np.min(_data(m)) for m in _meshes]
@@ -170,9 +157,9 @@ def convergence_metrics(
             / np.linalg.norm(_data(reference), axis=0, ord=2)
         ]
     data = np.column_stack(
-        (el_lens, maxs, mins, rel_errs_max, rel_errs_min, rel_errs_l2)
+        (x, maxs, mins, rel_errs_max, rel_errs_min, rel_errs_l2)
     )
-    columns = ["mean element length", "maximum", "minimum"] + [
+    columns = [x_str, "maximum", "minimum"] + [
         f"rel. error ({x})" for x in ["max", "min", "L2 norm"]
     ]
 
@@ -182,8 +169,9 @@ def convergence_metrics(
 def log_fit(x: np.ndarray, y: np.ndarray) -> tuple[float, np.ndarray]:
     if np.all(np.isnan(y)):
         return 0.0, y
-    _x = x[np.invert(np.isnan(x))]
-    _y = y[np.invert(np.isnan(y))]
+    indices = np.invert(np.isnan(x))
+    _x = x[indices]
+    _y = y[indices]
     params = np.polyfit(np.log10(_x), np.log10(_y), 1)
     fit_vals = 10 ** (params[0] * np.log10(x) + params[1])
     return params[0], fit_vals
@@ -192,8 +180,12 @@ def log_fit(x: np.ndarray, y: np.ndarray) -> tuple[float, np.ndarray]:
 def plot_convergence(metrics: pd.DataFrame, property: Property) -> plt.Figure:
     "Plot the absolute values of the convergence metrics."
     fig, axes = plt.subplots(2, 1, sharex=True)
-    metrics.plot(ax=axes[0], x=0, y=1, c="r", style="-o", grid=True)
-    metrics.plot(ax=axes[1], x=0, y=2, c="b", style="-o", grid=True)
+    metrics.iloc[:-1].plot(ax=axes[0], x=0, y=1, c="r", style="-o", grid=True)
+    axes[0].plot(metrics.iloc[-1, 0], metrics.iloc[-1, 1], "r^")
+    axes[0].legend(["maximum", "Richardson\nextrapolation"])
+    metrics.iloc[:-1].plot(ax=axes[1], x=0, y=2, c="b", style="-o", grid=True)
+    axes[1].plot(metrics.iloc[-1, 0], metrics.iloc[-1, 2], "b^")
+    axes[1].legend(["minimum", "Richardson\nextrapolation"])
     y_label = property.output_name + " / " + property.output_unit
     fig.supylabel(y_label, fontsize="medium")
     fig.tight_layout()
@@ -203,6 +195,7 @@ def plot_convergence(metrics: pd.DataFrame, property: Property) -> plt.Figure:
 def plot_convergence_errors(metrics: pd.DataFrame) -> plt.Figure:
     "Plot the relative errors of the convergence metrics in loglog scale."
     plot_df = metrics.replace(0.0, np.nan)
+    plot_df[plot_df < 1e-12] = np.nan
     x_vals = plot_df.iloc[:, 0].to_numpy()
     fig, ax = plt.subplots()
     for i, c in enumerate("rbg"):

ogstools/studies/convergence/convergence_study_template.md

----
-jupyter:
-  jupytext:
-    text_representation:
-      extension: .md
-      format_name: markdown
-      format_version: '1.3'
-      jupytext_version: 1.15.2
-  kernelspec:
-    display_name: .venv
-    language: python
-    name: python3
----
-
-# Convergence study
-
-This script performs a convergence study and generates plots to analyze the
-convergence of numerical simulations.
-Below are the custom parameters for this report.
-
-```python tags=["parameters"]
-mesh_paths: str = ""
-timestep: int = 0
-property_name: str = ""
-```
-
-Import required modules and setup plots.
-
-```python
-import pyvista as pv  # noqa: E402
-
-import ogstools.meshplotlib as mpl  # noqa: E402
-from ogstools.meshlib import MeshSeries # noqa: E402
-from ogstools.propertylib import THM  # noqa: E402
-from ogstools.studies.convergence import (  # noqa: E402
-    convergence_metrics,
-    plot_convergence,
-    plot_convergence_errors,
-    richardson_extrapolation,
-)
-
-mpl.setup.reset()
-mpl.setup.show_element_edges = True
-mpl.setup.ax_aspect_ratio = 1
-mpl.setup.fig_scale = 0.5
-```
-
-Read the meshes, get Property object from property name, define topology and calculate Richardson extrapolation.
-
-```python
-meshes = []
-for mesh_path in mesh_paths:
-  if mesh_path.split(".")[-1] in ["xdmf", "pvd"]:
-    meshes += [MeshSeries(mesh_path).read(timestep)]
-  else:
-    meshes += [pv.read(mesh_path)]
-     
-mesh_property = THM.find_property(property_name)
-mesh_property.output_unit = ""
-mesh_property.data_unit = ""
-topology = meshes[-3]
-richardson = richardson_extrapolation(meshes, mesh_property, topology)
-```
-
-Plotting the grid convergence.
-
-```python
-fig = mpl.plot(richardson, "grid_convergence")
-```
-
-Plotting the 3 finest discretizations.
-
-```python
-fig = mpl.plot(meshes[-3:], mesh_property)
-```
-
-Plotting the Richardson extrapolation.
-
-```python
-fig = mpl.plot(richardson, mesh_property)
-```
-
-Table of convergence metrics.
-
-```python
-mpl.core.plt.rcdefaults()
-metrics = convergence_metrics(meshes, richardson, mesh_property)
-metrics.style.format("{:,.5g}").hide()
-```
-
-Absolute convergence metrics.
-
-```python
-fig = plot_convergence(metrics, mesh_property)
-```
-
-Relative errors in loglog-scale.
-
-```python
-fig = plot_convergence_errors(metrics)
-```

ogstools/studies/convergence/examples/__init__.py

-import importlib.resources as pkg_resources
-from copy import deepcopy
-
-import numpy as np
-import pyvista as pv
-
-mesh_paths = [
-    str(pkg_resources.files(__name__) / f"square_1e0_neumann_{i}.vtu")
-    for i in range(6)
+from importlib import resources
+
+from .steady_state_diffusion import (
+    analytical_solution as steady_state_diffusion_analytical_solution,
+)
+
+_prefix = resources.files(__name__)
+steady_state_diffusion_prj = _prefix / "steady_state_diffusion.prj"
+nuclear_decay_prj = _prefix / "nuclear_decay.prj"
+nuclear_decay_bc = _prefix / "decay_boundary_conditions.py"
+
+__all__ = [
+    "steady_state_diffusion_analytical_solution",
+    "steady_state_diffusion_prj",
+    "nuclear_decay_prj",
+    "nuclear_decay_bc",
 ]
-meshes = [pv.read(mesh_path) for mesh_path in mesh_paths]
-
-
-def _c_k(k):
-    return 0.5 * (2 * k - 1) * np.pi
-
-
-def _a_k(k):
-    return 2 / (_c_k(k) ** 2 * np.cosh(_c_k(k)))
-
-
-def _h(points):
-    result = np.ones(len(points))
-    for k in np.arange(1, 100):
-        c_k_val = _c_k(k)
-        sin_c_k = np.sin(c_k_val * points[:, 1])
-        sinh_c_k = np.sinh(c_k_val * points[:, 0])
-        result += _a_k(k) * sin_c_k * sinh_c_k
-    return result
-
-
-def analytical_solution(target_mesh: pv.DataSet) -> pv.DataSet:
-    new_mesh = deepcopy(target_mesh)
-    new_mesh.clear_point_data()
-    points = new_mesh.points
-    new_mesh.point_data["pressure"] = _h(points)
-    return new_mesh

ogstools/studies/convergence/examples/decay_boundary_conditions.py

+from typing import TYPE_CHECKING
+
+try:
+    import ogs.callbacks as OpenGeoSys
+except ModuleNotFoundError:
+    if not TYPE_CHECKING:
+        import OpenGeoSys
+
+import ogstools.physics.nuclearwasteheat as nuclear
+
+boundary_len = 1500  # m
+repo_edge_len = 1500  # m
+
+
+class T_RepositorySourceTerm(OpenGeoSys.SourceTerm):
+    def getFlux(self, t, coords, primary_vars):  # noqa: ARG002
+        value = (
+            nuclear.repo_2020_conservative.heat(t)
+            / repo_edge_len
+            / boundary_len
+            / 2.0  # due to symmetry
+        )
+        derivative = [0.0] * len(primary_vars)
+        return (value, derivative)
+
+
+T_source_term = T_RepositorySourceTerm()

ogstools/studies/convergence/examples/nuclear_decay.prj

+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <meshes>
+        <mesh>square_domain.vtu</mesh>
+        <mesh>square_physical_group_left.vtu</mesh>
+        <mesh>square_physical_group_bottom.vtu</mesh>
+        <mesh>square_physical_group_right.vtu</mesh>
+    </meshes>
+	<python_script>decay_boundary_conditions.py</python_script>
+    <processes>
+        <process>
+            <name>HeatConduction</name>
+            <type>HEAT_CONDUCTION</type>
+            <integration_order>2</integration_order>
+            <process_variables>
+                <process_variable>temperature</process_variable>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable internal_name="heat_flux" output_name="heat_flux"/>
+            </secondary_variables>
+        </process>
+    </processes>
+    <time_loop>
+        <processes>
+            <process ref="HeatConduction">
+                <nonlinear_solver>basic_picard</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <reltol>1.e-6</reltol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>FixedTimeStepping</type>
+                    <t_initial> 0.0 </t_initial>
+                    <!-- 180 yrs -->
+                    <t_end> 5680368000 </t_end>
+                    <timesteps>
+                        <pair>
+                            <repeat>1</repeat>
+                            <!-- 30 yrs -->
+                            <delta_t>946728000</delta_t>
+                        </pair>
+                    </timesteps>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>nuclear_decay</prefix>
+            <timesteps>
+                <pair>
+                    <repeat> 1 </repeat>
+                    <each_steps> 1 </each_steps>
+                </pair>
+            </timesteps>
+            <variables>
+                <variable> temperature </variable>
+                <variable> heat_flux </variable>
+            </variables>
+            <suffix>_ts_{:timestep}_t_{:time}</suffix>
+        </output>
+    </time_loop>
+    <media>
+        <medium id="0">
+            <phases/>
+            <properties>
+                <property>
+                    <name>thermal_conductivity</name>
+                    <type>Constant</type>
+                    <value>2</value>
+                </property>
+                <property>
+                    <name>specific_heat_capacity</name>
+                    <type>Constant</type>
+                    <value>1000</value>
+                </property>
+                <property>
+                    <name>density</name>
+                    <type>Constant</type>
+                    <value>2000</value>
+                </property>
+            </properties>
+        </medium>
+    </media>
+    <parameters>
+        <parameter>
+            <name>T0</name>
+            <type>Constant</type>
+            <value>300</value>
+        </parameter>
+        <parameter>
+            <name>HeatSource</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>temperature</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>T0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>square_physical_group_right</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>T0</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+            <source_terms>
+                <source_term>
+                    <mesh>square_physical_group_left</mesh>
+					<type>Python</type>
+                    <source_term_object>T_source_term</source_term_object>
+                </source_term>
+            </source_terms>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_picard</name>
+            <type>Picard</type>
+            <max_iter>10</max_iter>
+            <linear_solver>direct_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>direct_linear_solver</name>
+            <eigen>
+                <solver_type>SparseLU</solver_type>
+                <scaling>true</scaling>
+            </eigen>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>

ogstools/studies/convergence/examples/readme.md

-# how to generated the meshes
-
-```shell
-max=7
-for (( i=0; i <= $max; ++i ))
-do
-    n=$((2**$i))
-    generateStructuredMesh -e quad -o square_$i.vtu --lx 1 --ly 1 --lz 0 --nx $n --ny $n --nz 1
-done
-```
-
-# how to run the simulations
-
-```python
-from ogs6py import ogs
-
-from ogstools.meshlib import MeshSeries
-
-for i in range(8):
-    model = ogs.OGS(
-        INPUT_FILE="./square_neumann.prj",
-        PROJECT_FILE=f"./square_neumann_{i}.prj",
-    )
-    model.replace_mesh("square.vtu", f"square_{i}.vtu")
-    model.replace_text(f"square_1e0_neumann_{i}", "./time_loop/output/prefix")
-    model.write_input()
-    model.run_model(logfile=f"out_{i}.log")
-
-base_id = 3
-base_mesh = MeshSeries("./square_1e0_neumann_0.pvd").read(-1)
-base_mesh.save("./square_neumann_convergence_study_res_0.vtu")
-for i in range(base_id + 1, 8):
-    mesh = MeshSeries(f"./square_1e0_neumann_{i-base_id}.pvd").read(-1)
-    base_mesh = base_mesh.sample(mesh)
-    base_mesh.save(f"./square_neumann_convergence_study_res_{i-base_id}.vtu")
-```

ogstools/studies/convergence/examples/square_1e0_neumann_0.vtu

-<?xml version="1.0"?>
-<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
-  <UnstructuredGrid>
-    <FieldData>
-      <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="20" format="appended" RangeMin="45"                   RangeMax="103"                  offset="0"                   />
-    </FieldData>
-    <Piece NumberOfPoints="4"                    NumberOfCells="1"                   >
-      <PointData>
-        <DataArray type="Float64" Name="pressure" format="appended" RangeMin="1"                    RangeMax="1.75"                 offset="84"                  />
-        <DataArray type="Float64" Name="v" NumberOfComponents="2" format="appended" RangeMin="1.1938124471e-16"     RangeMax="1.0606601718"         offset="156"                 />
-      </PointData>
-      <CellData>
-      </CellData>
-      <Points>
-        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="1.4142135624"         offset="252"                 />
-      </Points>
-      <Cells>
-        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="324"                 />
-        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="396"                 />
-        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="456"                 />
-      </Cells>
-    </Piece>
-  </UnstructuredGrid>
-  <AppendedData encoding="base64">
-   _AQAAAAAAAAAAgAAAAAAAABQAAAAAAAAAHAAAAAAAAAA=eF4z0zPRM9E1NDQ31E1PTjEyNjdPSQMALjgE6g==AQAAAAAAAAAAgAAAAAAAACAAAAAAAAAAFQAAAAAAAAA=eF5jYACBD/YMaPR/MPhtDwBiYQrCAQAAAAAAAAAAgAAAAAAAAEAAAAAAAAAAJQAAAAAAAAA=eF5jYACBKTZgimEOlN5i8+M/CDzf/wdKw8T/Qfm/oDQABQ0h+Q==AQAAAAAAAAAAgAAAAAAAAGAAAAAAAAAAFQAAAAAAAAA=eF5jYMAHPtjjlcaQh/ER4gCW5AS9AQAAAAAAAAAAgAAAAAAAACAAAAAAAAAAEwAAAAAAAAA=eF5jYIAARijNDKWZoDQAAHgABw==AQAAAAAAAAAAgAAAAAAAAAgAAAAAAAAACwAAAAAAAAA=eF5jYYAAAAAoAAU=AQAAAAAAAAAAgAAAAAAAAAEAAAAAAAAACQAAAAAAAAA=eF7jBAAACgAK
-  </AppendedData>
-</VTKFile>

ogstools/studies/convergence/examples/square_1e0_neumann_1.vtu

-<?xml version="1.0"?>
-<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
-  <UnstructuredGrid>
-    <FieldData>
-      <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="20" format="appended" RangeMin="45"                   RangeMax="103"                  offset="0"                   />
-    </FieldData>
-    <Piece NumberOfPoints="9"                    NumberOfCells="4"                   >
-      <PointData>
-        <DataArray type="Float64" Name="pressure" format="appended" RangeMin="1"                    RangeMax="1.6857142857"         offset="84"                  />
-        <DataArray type="Float64" Name="v" NumberOfComponents="2" format="appended" RangeMin="1.9053550706e-16"     RangeMax="1.1571428571"         offset="192"                 />
-      </PointData>
-      <CellData>
-      </CellData>
-      <Points>
-        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="1.4142135624"         offset="380"                 />
-      </Points>
-      <Cells>
-        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="480"                 />
-        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="576"                 />
-        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="648"                 />
-      </Cells>
-    </Piece>
-  </UnstructuredGrid>
-  <AppendedData encoding="base64">
-   _AQAAAAAAAAAAgAAAAAAAABQAAAAAAAAAHAAAAAAAAAA=eF4z0zPRM9E1NDQ31E1PTjEyNjdPSQMALjgE6g==AQAAAAAAAAAAgAAAAAAAAEgAAAAAAAAALwAAAAAAAAA=eF5jYACBD/YMOOj4j6JA9Nle5oqj7BXHn3DxAPZSIPpi37n+R9f6H7/sAUWjFeg=AQAAAAAAAAAAgAAAAAAAAJAAAAAAAAAAagAAAAAAAAA=eF5jYACBFhswxbB6D4ReYPOha/2PrvU39kP4W2xMVjRbrWj+tP8LinjJHjOw+KP9VrJXHGWvXNz/GSz/Yr9o/EcgerrfDiz+EKp+iQ0PVJwLTB/df28B3ysg2n9tW+7tbbmn9wMAeYNAHA==AQAAAAAAAAAAgAAAAAAAANgAAAAAAAAAKAAAAAAAAAA=eF5jYMAHHthjF/+AQxwG0PXB+OjiH3CIo8vDAEwduvgHDHEAmlgN1Q==AQAAAAAAAAAAgAAAAAAAAIAAAAAAAAAAJgAAAAAAAAA=eF5jYIAARijNAqWZ0cSZoDQrDnUwPjuUZkMTh+njQFMHAA1AAEE=AQAAAAAAAAAAgAAAAAAAACAAAAAAAAAAEwAAAAAAAAA=eF5jYYAADijNA6UFoDQAAqAAKQ==AQAAAAAAAAAAgAAAAAAAAAQAAAAAAAAADAAAAAAAAAA=eF7j5OTkBAAAXgAl
-  </AppendedData>
-</VTKFile>

ogstools/studies/convergence/examples/square_1e0_neumann_2.vtu

-<?xml version="1.0"?>
-<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
-  <UnstructuredGrid>
-    <FieldData>
-      <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="20" format="appended" RangeMin="45"                   RangeMax="103"                  offset="0"                   />
-    </FieldData>
-    <Piece NumberOfPoints="25"                   NumberOfCells="16"                  >
-      <PointData>
-        <DataArray type="Float64" Name="pressure" format="appended" RangeMin="1"                    RangeMax="1.6779678311"         offset="84"                  />
-        <DataArray type="Float64" Name="v" NumberOfComponents="2" format="appended" RangeMin="2.1138016631e-16"     RangeMax="1.5877798805"         offset="336"                 />
-      </PointData>
-      <CellData>
-      </CellData>
-      <Points>
-        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="1.4142135624"         offset="876"                 />
-      </Points>
-      <Cells>
-        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="1028"                />
-        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="1208"                />
-        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="1320"                />
-      </Cells>
-    </Piece>
-  </UnstructuredGrid>
-  <AppendedData encoding="base64">
-   _AQAAAAAAAAAAgAAAAAAAABQAAAAAAAAAHAAAAAAAAAA=eF4z0zPRM9E1NDQ31E1PTjEyNjdPSQMALjgE6g==AQAAAAAAAAAAgAAAAAAAAMgAAAAAAAAAmgAAAAAAAAA=eF5jYACBD/YMRNJ1Mz+yF134YL9BV8/W4exHe07N7X0THT/bi24S+PEo8htc3cT0iSK61R/tN5qdfatk/dl+J09WwLqFX+0VVZ2ZZnD8RJhnpmjV+fyjvSpHrTmX5Bf7cKGfh1PufLPf9Ozy5cnZv+DqjqS8T7Lm+WR/foWdyZaEL/b77a9fLrX/bn8y6PnpLzd+2QMAFAlK6w==AQAAAAAAAAAAgAAAAAAAAJABAAAAAAAAcwEAAAAAAAA=eF5jYACBGXvAFMMaGzDVMMdmib7xPwe+U/sh/I02P249usR85Q6Uv83GzmNnaQ/PKwif4ZBNwAmnB93pP/dvQdbHMN3mP1jfmf0slWKeFzcd3x9do2P2p+Py/klbM/PEb97cP0Ww63WX0aP9D4xWOO5se7Y/0eOMV/ehF/tdvNpYDAQ/7ecF67sONa9sTxtUH8OFn0rFhQf3L5qS9J5l1oP9sy6XXfKPu7R//87rgVdMn++v8uJfNfPanf02h0/0+Zq/3p91iSf/sNuD/bfShRdNrrgHMS9g2p6JCl/uaqQ82B+4bP/pk4qb9ifkx6vnXX68/4KieLpI5KH9OlXRhWpeL/fnGlw2q7hzcv+FRdfqKirf7N9zy2nWGdaz++XUqsVvJz6AmOewao/V4l8XLjQ93M+lDNK/aP+b7PmHX555sn+XVoZB/auN+18sXxTH9PDl/l2R3d+dj+zYHzKj4cD302/2X2K9zWIXsXs/AGtuwjY=AQAAAAAAAAAAgAAAAAAAAFgCAAAAAAAAUAAAAAAAAAA=eF590LERwDAMQlFvlv27jEBJyQipaEh+6PzkO0465y/39e0CN3jAm+3pe13gBg94s3v037rADR7wZu/U3nWBGzzgO2/auy5wg+flD1SCLSk=AQAAAAAAAAAAgAAAAAAAAAACAAAAAAAAZAAAAAAAAAA=eF5djTcOwFAMhX567/X+J80Qs+DlSQjhlP7LYuvYSjyPbeTBi9hWHryM7eTxh94Q24vTG+XB6U3y4Pyd5fGH3hq7iNPb5MHp7fLg/D3k8YfeFXuK07vlwek98uD8feV9Nx4DAQ==AQAAAAAAAAAAgAAAAAAAAIAAAAAAAAAAMwAAAAAAAAA=eF5jYYAADijNA6UFoLQIlJaA0jJQWgFKq0BpDSitA6UNoLQJlLaA0jZQ2gFKAwBmgAIhAQAAAAAAAAAAgAAAAAAAABAAAAAAAAAACwAAAAAAAAA=eF7j5EQFAATYAJE=
-  </AppendedData>
-</VTKFile>

ogstools/studies/convergence/examples/square_1e0_neumann_3.vtu

-<?xml version="1.0"?>
-<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
-  <UnstructuredGrid>
-    <FieldData>
-      <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="20" format="appended" RangeMin="45"                   RangeMax="103"                  offset="0"                   />
-    </FieldData>
-    <Piece NumberOfPoints="81"                   NumberOfCells="64"                  >
-      <PointData>
-        <DataArray type="Float64" Name="pressure" format="appended" RangeMin="1"                    RangeMax="1.6759650555"         offset="84"                  />
-        <DataArray type="Float64" Name="v" NumberOfComponents="2" format="appended" RangeMin="1.2807441219e-15"     RangeMax="2.0260704397"         offset="872"                 />
-      </PointData>
-      <CellData>
-      </CellData>
-      <Points>
-        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="1.4142135624"         offset="2536"                />
-      </Points>
-      <Cells>
-        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="2824"                />
-        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="3308"                />
-        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="3556"                />
-      </Cells>
-    </Piece>
-  </UnstructuredGrid>
-  <AppendedData encoding="base64">
-   _AQAAAAAAAAAAgAAAAAAAABQAAAAAAAAAHAAAAAAAAAA=eF4z0zPRM9E1NDQ31E1PTjEyNjdPSQMALjgE6g==AQAAAAAAAAAAgAAAAAAAAIgCAAAAAAAALgIAAAAAAAA=eF6tzN1LU2EAx3F1XkxkYgpluygjwUQwo+ZyZL9EbaamEjaVrDFErDXSi0gE8WWlNfVCY4UlmQsv7IWychbTrWXlUEMlWpJGmi+4UtlznnO2s8bE6mb+A31vvnefgIB/EQT8p7+rt2eIJAQxqlHGqSGQFcnvc4/+/vCS/hZLsN8juERPMqiLEI69tDGw/JJScQNFRdm4uFjE+R3WXlVXXk6wMaE2R3wmiG/J1ghyGBRMfDkRMs1gdUf/4ut6iua23sXrJSw6hL378ts5TKZ3dtMMt9+xpUUuGvUEKybD3iuJDIK/83tEswyqDrQ6A7soNFJBdlATi1f9Boerh0OidvNFTbgbC/a8XM033u+M3ZmLlAwSdMUk9ehqGThGLAPPT1EYnWEJt1NZIHe4jSo46KJXnsm7XZB8LNleFs8jb2hmPVTw2+9onprOv/9KcNenaEkYYtB0pEAm66NQ7p531L1hYfwge5jKcQhRTlFvkRtmhfr0IZ5HyrBnYcS+5cxcXjONOwmuptXI1Q4GsUMXj09SisL1WnlUFIczWbxSpXIhbidfmDntRmnKzxtRDR4QffMxW6nX7zS6lypEGwSf3pJ4NogiKTu3UidmsZzVUXwzh0Pjk6Ob1QYXzOFa7Y9oHrPCB6W+KQ/OPa68N2/ccizXRJnVgQzSjVbhaCiFdFffsimWRUvrtmRylsOctT2/c8CFg+oxfVgyj/ZgX3PcqgcXEtcGrVNe/AFgoyggAQAAAAAAAAAAgAAAAAAAABAFAAAAAAAAvwQAAAAAAAA=eF5Nk2tUEmYch3UezVmz1kXLaccuU48brax5QY66aKXZqNWaucrSLmbNDcvEkkrRSFtWHsWZl0pUMMsaTZsi+cMSM9SNQgUzLiISItpOSRtLayZ+6P3yP895PzzP+z/ntbJ6d04LJscsmWVaUUnR0UIDm8mHhUME88O2HR4LbJviVIHPtCF9blSXhUPYpNBpZ4pvH1ZO3bNIKkaw2kOlnWKGYMwz7SopxTjFN0n6rciR2Y9OsZBUxyIc2vullfDs+173cpL9pLcBHLpJoeqpR53i6I7uY3dxThVZHcEQo1Hzw56aTiGu3a5r9hV2gv0mV1a/qxk+Az+W5FEVyA9NKW9JFyMcV90e2Gph9AomZ/pIsTZjjVuw1gD6ry5Mb4IG0YuvskO+fgH7635Zmvt6+M+JZ4idzWie9D609EikpJIp75om6kZPwx8wc5IKT+S2Y5v/9MAWait+MqHlwBcSXN5h/1mmtRRv2jKTIsKlEOWeqhJxniB2iWvJYIQcOrP52UlhH6inr3jyM9QQftwxOsuggxub9Wxr+DNkNZDsHCf2FXd/5SnP9UO4lL6k+FLyCDqTD8QI8h9P7U0giL35CTFmmhQs098L9r2sQXtZ2MiRfZ1YWnOw5ylFBKs4bqNA3g1pw0qu9qAEyfm1FzffegInJwb5z8puWLvcaL3ip8aqff2Swj0KOH2vCZZu0cJ8ajT6jZsGmz5MIPzeMIg+V048V6KFe6BtrbnRCDvB6zixfgCrA16msPfLLD2/lQmUYTHkQr0M4deexJWSeeh69IIgpvUgwbfy1QxaE2J550palz8FSfLL9LIKMe4QXSV2i9WgK4NYKufHyGc2OZ4zaZCkTe439Xej6rn087yyib1UGNibQnrhrSx+vshpCPrQao/jh5TIs1bR4o8Po010fv01lgr0jq3O2f69lp6FnaSChz/X6u/14mlmPHNLeBUK+woW7IxVIKj+zCOvZD4+dQoznnVRQ8a9bOLZNCMxbBUt0H3i3TVvCWSJGBvG0zbK27TgN218aeRIQMykR7Oc9TBlvU11oErBXkBO26Ubgnz2Ay+/7V1wqNPJvQaH4X123Lt2ejdWFJ3fSTMrLD3fiQTe/s20cqYSzjOWstLLS7Gf+N921/lqNH6T4OHF46GIMs+KeqIPNjfST3bk8lGQsKP9YmI/Kqi2zhR2E7KimtJTzQMY7aELHbgiEAsurKsf16PS1/fVo7FWLFMWrflnjxGV5Fct1gFtuJc87/iQ9wiWLx2n+Nm0o0Ryeay6VGXpCb8jyHZk7rrkrEZv/8xF2dU5OFOx3kY9pEa22p+9+1bpxP8hxZs3aGBjv4L919EqJDqk0zL8tCjq4O4dO8qDcBll3TG2Dg6lJR50Si1uaVyWzKEPgri6c9WFhDrYdSXK/uUa0dnGywjJ4MPHy1BMpowgQLlfHhXUAGW4PDVvwj/ZE6oguVfZrg2NUmPOwDt/GkRtc4dn+PSBodVdqb+eA18+kxGQo0GLNee+IaYQxHkfRdjQtIjjfHVMtbMUQcMtihdyHXZHzt3y7aIK5LttliaUD4KRFBk6m1CJ6yscA3LuGrGy9nWkJ6EKKXc/mHkkcgRFDqdTFqqr8D+cdapRAQAAAAAAAAAAgAAAAAAAAJgHAAAAAAAAtgAAAAAAAAA=eF591KEVQzEMQ9G/WVfvCIaGBgYGBh6hSETtq1huQnQS53n+5f367QGe4AXe4AO+4AeueA+t3QM8wQu8wQd8wQ9c8XvQOfcAT/ACb/ABX/ADV/wdqZd7gCd4gTf4gC/4gSs+B+rlHuAJXuANPuALfuCKz7F6uQd4ghd4gw/4gh+44v+QerkHeIIXeIMP+IIfuOL/qHq5B3iCF3iDD/iCH7jvK+rlHuAJXuANPuALfl/+Aa9FoFE=AQAAAAAAAAAAgAAAAAAAAAAIAAAAAAAASQEAAAAAAAA=eF5dz8VSHQAUBFEsaHB3d5cQnAcEkuD//zUsuGczs5mqrt50U9P3muu767uCt9T3hIe31v8MD2+r7w0P/1HfFx7eXt8fHt5RPxAe3lk/GJ5OPcP1Q8H1jISH6xkND9czFh6uZzw8XM9EeLieyfBw3VPh6dQzUz8dXM9seLieufBwPfPh4XoWwsP1LIaH61kKD9e9HJ5OPav1K8H1rIWH61kPD9ezER6+UL8ZHq5nKzxcz3Z4uO6d8HTq2avfDa5nPzxcz0F4uJ7D8HA9R+Hheo7Dw/WchIfr/hWeTj2/60+D6zkLD9dzHh6u5yI8XM9leLieq/BwPdfh4bpvwtOpp1F/G1zPXXi4nvvwcD0P4eF6/oSH63kMD9fzFB6u+294Ohv1/+v/BdfzHB6u5yU8XM9reLiet/BwPe/h4Xo+wsN1f4b3BaKQKAE=AQAAAAAAAAAAgAAAAAAAAAACAAAAAAAAmQAAAAAAAAA=eF4txRF0QgEAAMB6LwiCIAiCIAiCQRAEwSAIgiAYBEEwGARBEARBEAyCIAiCIAiCIBgEQRAEQRAEQRAMguDfyUVCgahjjjvhpFNOO+Osc/5w3gUXXfKny6646prr/nLDTbf87R+33XHXPfc98NAj/3rsiaeeee6Fl1557Y23/vPOex989MlnX3z1zXc//PS/Xw6Fg95TUR+CAQAAAAAAAAAAgAAAAAAAAEAAAAAAAAAADAAAAAAAAAA=eF7j5KQMAABJYAJB
-  </AppendedData>
-</VTKFile>

ogstools/studies/convergence/examples/square_1e0_neumann_4.vtu

-<?xml version="1.0"?>
-<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
-  <UnstructuredGrid>
-    <FieldData>
-      <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="20" format="appended" RangeMin="45"                   RangeMax="103"                  offset="0"                   />
-    </FieldData>
-    <Piece NumberOfPoints="289"                  NumberOfCells="256"                 >
-      <PointData>
-        <DataArray type="Float64" Name="pressure" format="appended" RangeMin="1"                    RangeMax="1.675476229"          offset="84"                  />
-        <DataArray type="Float64" Name="v" NumberOfComponents="2" format="appended" RangeMin="1.1281095627e-15"     RangeMax="2.4665711682"         offset="2780"                />
-      </PointData>
-      <CellData>
-      </CellData>
-      <Points>
-        <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0"                    RangeMax="1.4142135624"         offset="8556"                />
-      </Points>
-      <Cells>
-        <DataArray type="Int64" Name="connectivity" format="appended" RangeMin=""                     RangeMax=""                     offset="9332"                />
-        <DataArray type="Int64" Name="offsets" format="appended" RangeMin=""                     RangeMax=""                     offset="11096"               />
-        <DataArray type="UInt8" Name="types" format="appended" RangeMin=""                     RangeMax=""                     offset="11676"               />
-      </Cells>
-    </Piece>
-  </UnstructuredGrid>
-  <AppendedData encoding="base64">
-   _AQAAAAAAAAAAgAAAAAAAABQAAAAAAAAAHAAAAAAAAAA=eF4z0zPRM9E1NDQ31E1PTjEyNjdPSQMALjgE6g==AQAAAAAAAAAAgAAAAAAAAAgJAAAAAAAAwwcAAAAAAAA=eF7Nlfk/lukexz0NNRExLSNabB2lkjWp9CkHJ1IeemLspVGWlKXlHJSorFOhxSMRUsfSIIVKUpkQRQuiSaRUtLju5VloMPf5hX/h3L/cr/un+/35XN/v+5KR+d9DIPN/8jZVHjgvUCRoNheK2tUI3HnbCioWEXz0UrfQNSM4YLGmyseGIN890fB3Z4KUM9vvu/oRJNG/7XgeQXDCsiEqL4XAnzHMOFdIELs5a8fTeoLUDbLhSz4TBL56faP/ZwpLg/jz8gQUOt3EexyvUQjJtTINOUiPc2wr71aM5P4fltrxwtGUQHV/TewbPkFCwI3q2CCOr6MRVUkEHV/vBI0VEWRaLfOhmwmm6Chc8RAR7BbIvbDWoDCjIjRW1ZFC1VjphbcJFNxjhO/rnlA4mu65YkiTRhfvAlt9ikaysCJDX4/BraVNZxJk2XGOTcP+H+10CaxaV2qrOhBoO2Q//HaQYO236G7vXAKTJeiZ1kJQ7GHZXzJKcCnDWqBrREFkluMfG0Chvz9eqF9A4d3h7wPiQQr6ex/cD1xPo8l+QLAyh4apfUL+85kM8ubPSD6Sx+BUi4LrWTcWNhaXCkrXicY55ncGM0JjAnR9tezYRRDQ1/fEP5OgLVqXJ9fG5e54N2nfdArpR78+urCJQkvlFPMTydx3hVd7UieFDKW0/cv1aFzLM9SJPkYjqicraOlnGpEk2+TedgZOdxe8E31hkLV42Cw/lUVotXmmqqsIu0Kdz1hai8c51k26EmgJrn81fq9tFEEGeWP29h6BUcmpDd2yFNo25Brp21PItm/VbU3jenhhmTH6iQLf+Fq1Mpe7cjEjqOByF6t7lX2fxmDqnBrK/ziDWWNeezSVWFT+wZp0FbKY6xkfccdbBPE+p5FSAzEkjj0qVzQk4xzecjW3FtkRhK9+kLj8LIHilBX/kL4leGCs8yzEkIJsXaly1jEKJ8N69glfU7D67tR3cRWNzrtT72tm0/jDd+hZnRIDx2Bvk+h4BrkGh9cKp3PnnvbFuOgKi9e1/xFk8UUwO5fZt0NZjLKS5MoLn8TQcg111vxzgiNO4PPQSEDw56Gu5S4FBDGDU77vkKEwK2fhBjc3Lr/LX+4bb1KYu9FV6qXO5W8mO19x/Q+ZZ8xcJKYx5n5StzuYAZ+Ypj8SM3BeX9LRl8DCT9VqivdSEejfdqbt7xEhcEn14y35Ykwa/Kb09LgEjx2U1BaHS8c5CjRvvbT3JKirHX6TXkUQf2h25md1Cr9eVvJxiKJwLFPa2/2Bwhzt//7QuIVG8xrXYsd6Gpo7BJ7tlgxW9TmnttczONhx+iVxYVGvUarWw7JofkolBWeL8ObgMY8ADzHe999MDtGToDesuKdLQYrhsUtLUnlDE/5Qmq/n50vw3i+45lUTga6NQYGvCYWoh6NN87IpfFVXjNRSpqE2Gvck+ziNjDKFuHYegyr3rbYzYxl0Hj244u4MFkHxuT/aXWVxWVvezNxRhKadC28pyopRGvn0yMo6MX7kq2tNSpfAZvCsUfAhKUoUikaX7JvgUOPNWx2+m9sbrZKLcq8Iju3OGiq3pRB54rJN3h0K0dNdVCljGi9UzS7WldI4+9fFuW4m3P75CJy/1DD4qVY/skHAwpDlNazjepgzuV4p/6II8f6JOf/6RQy56g5PxwUSnGwrLWsWSTBYNHr45GsprKF/eWbrBEfdsq1ykaEEH9ZvtVP8QLD3gOBDozuForit0S2tFEwa+Y2rnGhY5eVEGLXRCIzSL+71YjDfvEr3EGGgIv+ThVkSixuz37VZG4hQGiez9CU3D9WaYhX5LDHiOu6QXj+uh4i61bCUotK4sKJBbwhZfFu3Aa3hcQ4lNwcDf85bouKVbf2DBHbPFU+c5/wUY1+jc47bT9fWLIfhAG4uN4/0fqRprEuu9YyJZtBSpxfx7WcW1nVxiQ6VLO673rar3SZCkM3jeR6zxFjOc0yb0y5GtPnVku5LEuwiqgiKkiLJN+qmk/8QdKqa+APbJzhWS99uAndPhN0qn54mIci6/uwXjQMUHvQPbhFJKDx2SQkwPEwjsyV4kExlMDzy61z+BQat608b+a9gkdI1+5RjJwtj4WjKrFjOj5fCNl23EEOvS37OTJ4EPnFdLnLPJLCXOScxKJZictuw5+z0IQwUNaXIpk5w1HjQBiOHCC7ndi+aMUIgb5064HaYwjyr4Iavk2h41v2u25fEeSKx9JGrGoMQreehPiWcp5LKG1TsWazx07ALo1hsXdnLj+P247u568K1bmJ8PnnELpSbi659Ppt7iQRl+k6n9Fqk4IW9ia65PYRZtj3T+q9PcNzt8PHN5Tz66ebNDE/OGzGRTJjdUQpO5fPToqdyfv72Wkn9DA2VEm+XZToMbAq3tTfdZtB3hieUdWURE7JhTdcYC+G9f093KBUhpm2Tp2Ug5yk3bb0EQwnCUyLOd/4gRbtQe3dDrxRXJv9TXNs8hL1jI4uUGyY4Ck9vfqJ/hIAXbuF3gMd5a0FU4iDnT5XShVllCty9rPHU6OU5GuaFeSOhugyUc9XWC+8yyE+4GrXFk0WPLbszX1aEbWbB5ES5CK+2723sDxbji3iFaZKZBCEv9mQslJeiPuO21cZ+KXrzc5L92obgG5Sw+F7zMP4Gp/dXtQ==AQAAAAAAAAAAgAAAAAAAABASAAAAAAAAyxAAAAAAAAA=eF5NV3k01d37FZKEZKgMRW+DVEqJSrcoRRlKMiVjGimRvFRCmRUZIiSRZMpciJttnsV1zeOdZCZDIdLvfl/++H3+OWuvddZZn73Pc579bBYW5ucKwv8WFmLj0nrmEwGWEk1yiUn4D/+IJcp/OxL2syVvCbN8IUaksKTrXC5dwq7uxGt2aX5RwTVLWOkDIcp4bX1qGGl5/xvC9GDoFY3k5iUsEUPk3XbNwGZ9xxKOTiFWviPQw6Z7lrBZMtHrAufci4O0ZRxL3HTn+HWydt/yea+JiQ1Vdty/B5axH/Ekv4/S3WujS7jwE+HjhzrN184Ty/+XRchu5Gg+ofZzeX8SwSPl2vh6ofllnEvwFNie5JPPVmjz/3nb5BHU/uOdBAmzz708TklYsWPCnvdKMnYdOTsymJCHNxu3u6lqfcTE6e60mEOlGBYhdK37mAKldyJ30sxrcKBowXraIw3lwpfnbjiQcN0vmd9tPh2zeW9yRl2accdWfUV4UCYKghNOa/S0g/pHT9Y07hOaXK8WvUzuAdXatVH0Sg5yVHRljNlpsF5YOTwhlA8p8903vNf1QUYt5NOh7YVInV9bIxw1ANm5d+d1OEoRud37rHXXCDx3rCbvNKmGYO1IZtXaCZTyyCjP7CGjl5fyfvzxNLR0lF1Ie2hwLP3Uff/CHFZMVVbNlgyAkrnn4g5+lkL2/3gv3zt7PdF9mXf9hiEXjV+JEOnqvUvYnw/Fo5oNoRJ5WLlmZPS+CBEdtgJ8q76VYLIjw8PI+iuqmmPDd36rhnf8CXYDAlCbYxx9SoiElKmIsMq6QiRLsMzosTQj/pmfwDn9YojTXI0yj7QjeYpNvkunFAop3sWuLd2YfcoQnh0rR7E91/MYRypKyUOWk9erEZQTtsCux8B6fjufV7vr4c11kvY5vB+/DgfLTbGSEaFsllpnNoxdP+6y0I63IV2sQuHFt3HU/Hihq+dGQfEzoVC68yS+pEyy5mn0Iz0uWUucWS/iHkH/mKoPI6pl+t382V+I2PHnk11cyZIeMiTiCJ/cN87yEtC/bg3sE0uEo1yI38zOUnBdj/M9pPwF0c4q/OOkUtgIRlaR/imBNHtOvGdFGXoby7oCeqtAfi1qwypYgX/6xS02OjTgvWjkijfllVDI5fdMlm3CB1L2huif1dhz/VynoWwbvAPJfKSiOlzSvsf13aMLjd1ndvOlNoD6TlQqIoqCuAQZybhVZHQOjnX6yNMxQ2w/Jr2zBQe4e3hTBL+jJ7qzret2B35TPVlCFAcRsskv0+0gBT81/qmtThvBO7mdnLuNGGA/PiGiveEHeIdWhK8oGQS/pDC3lsMEiNmj/JYlI7CTNlyoWDmJb/63cx/UVi3pIUsj1J13yOf7VQWNyVUE95l4LPArv5G/X42IU7HPejfmwmz7Lf996jWgGCe8iCQUwzTeSdHvaS1GCDbKHW8rscrWdf/k9m8QvlkRpnOmHr9mH5Xu2NoA6cVjjV7byOAw2+FPfUyCGld870BHC5zbhwJVD5LhZPX3iRqpA9PFxnmJJs3Im1P1m+HrxUPCb5a34m1IqpX6V6OBCoMHgzy0+E6Uvan+S1NkQKabx+CVGAXvtqm/Cvz2HZuCguVceekQ984pTzAbBMvk25jdxd/h3n8xpyBrGKHPS8vipIZhrWd5Ks98FO+jQunHPUYhrE5y0/s7it/iK8hyMfVLekzXEmyyOF46Ntbj7Yk958pffMDVmk1e7bIN4Nm8byg2ORuOzycrPmY2ADc5BUO3FuHIKiuDt7dIcHq056+saAW+pQQPvDZpxJjzUTs3tTrMWKwT5Ikg4+uisYl9LwmjzZH5QZLN8Le+6aN7shnHc7nkZ9a3gpXTXPdrYxuIa96EmDxuh1yf2vPhoC4cHHocbvdvF9jyulhXKFFgJc9beJmfgn2sP3m+qtDQkmcoe1CfBqWbedNp5xi4SglpTBfug7tBtJCXwnf8YOv0IBgPwKlKviOvk/lONF2uhfOOYHHr3bgTswMQic4PfLc4ikSrgi0jxoMo0S+6u9e3cUkPKo2oclfVYCKuEYdVtImXfrxHrhzFf3SuEYKeGeUfOD5jhlOx9ow3GUF2Li/Ni4COHv3zX3Wa4O7o1AeTMtjNeTk0XW3Gnwo3PbVVNfgbr/JSobAFBvVv+QxtG2B4S8pDyrYNJgawf/CIjPHzdS+73TvA3azUPfSuBX/2ZluOrenG6rPXL4r0taM/RaNyPQ8Fqi9/dxCPdGPBsnH2LY2KXHfWuwWqFLyU/WJqzsfA3pAL/u10Km5OrRV99vo7vrmJSXpeosPoq/X+7y6DoDE0BOx9GIjYq8HTmjyCuugK4YkLfZiVlusK1hqDWnfvU5mwPqRsCTHI4F72QZZHxMv5afN9u5rx+4zY703bYnHNNOvTDs9mHN5xSmyXXBbkb7CvtpBoAdfxozb5c0To8/ptHllogfbNSIWrDiWYULLnCJVoQ0Pvz8D7AZU41+ZvQg1sRwctOcOitg5hnce8t17qhNj+SokjoSSwJWYIhLl2Y9FbxydYrwm215LPbuSioHYsbc3sllZc1+Xltwqgop7d7/r0aDucn9AoekyeoezHPecTulD45f7w/M0++B6JZ0k+3IvtMoH6GujHPqVk10IiBaxbPB3jdYfQovA5oyOaivUq/M06m5i+rNwVyM/088PCIxZbno3B+Y7IRZ4bNEhanhmTu9i6pEdtNjHV0aD3qF0rLNO+5F+Vi0HHo+TkNnIr3N3sqypOZKBA1ulK+r9tENHeLdjKkofHf5WJ3ibtMBE5+4tiU4So7LWyOi86sHpRhr3yZxkODwrcPsLbhVB8VLcPqMYCS4LcM2o3OOZyOTgV61Ho2aF4jpMCs8OnXnn+IcFCdyY2/gIVh+e3dVMvN4GUG+wR1kCDF1vwKE96C2579J6e92Kg+qxRjdVQG6xtT2967/gddrVShJafHTAItNXgix5A3swxaePcLsgXvj15eGYIL4JLk6uFeqBvnpG6+t4oPDkp99mZ887vpxkmPtlj0HB6FWyk1AslPzcvLsP2JT0KOggOYjxj+o/aQWAUVJQbRuH4ogotrbcd4pdMPhVapcIwuCbQ3LsD5Mczynv7spGbG3E5wqETHdz31e8lFcBSQHK/WXoX6p0j9TVPlaC2udH0iXwPnK3/0HcoVmBtRzPx0CoK1CfSAkVP1qBd4Ouvh5upzLr5q+h8uh5iNY3Tly1oWMibVzgWRULVgtTq9FY6JBxKeBZ9yBjuk+Y8ynyae5O6ni4aNONDwLXI6fP9uL1BIIv2uwXEYJMXKucHkf3eexfFqA1XpfTFHtkN46r9n+Pvrdphkzdf+SZzFB4rVlULMOdBLr/Ypw9ax+A4qLr/vFYHnFWUKvw3d/6nh6sbg1CizPiZpdSJ8K2K1tmSEZC/I93gwOznBqztDN8DSSBa1fU0XewCcevGoKjMLESZCrZaqHZjRMg8ucMgD82W+/nV3XvQGjgcIOBZCAWS/6kADgoifW95WjFKUPS6yfbyJAXTslmVq5Ur0HRjZLaPgwZfA4fwQxHV2CzSLqJ6mg4O/5gPJ3rqYCXjLjWXwgBr1bqEBy4NuKrrFXZH6TsMn9RGFUyRkGctOvHydz+ETDaIPFMkw6TwOr9t2yCmHz77o6nbhMAxE/5G8jC6b1yTi2X28duvWyRNGaPo8vN+ysecl79ZTpaV/hgDX1SlXU1rM67aquuueNu1VB++dGIqIUEtsagLm/szrg51h2Cj0KB4l1Q3Trc895FjxCE6l5Oi0taNKHpmztl76Qh39HVZqOvBn6e93OOzn+FXVfR+H5O/0m1a7Gv3fNhc7DzmXE2Bn4/fhJZ3ISI1w4uSPlHRPsQ4w6dWAr2kw2Nnymjg5AuX82Eph3b7Y5rvAh3FtUqhJTGV+NHeP7Z4qQ+TFsKp0RI1CEuXKVLu+g66oFYR25063OYpUye6DmB0MM728f16+IatfKOtPASvF3cVqgMbMJHZbGK0eQQKl9n8SLIkOOhKD2RwjcFJdoOMFzNPjEbS7Q6yj2O1U5rd5zwS3NxaBDQ6upf0UC0ixvQIxSv8YfIfeikc6BmIh6vqX+2+3oNZPwezMp8YBGSQf/gI9SJS+tjje7pJGJLgk1pYSYEE937NwaAM2EhS6KJxFKjVqEhozX+GXYE6zfUmFW0fQ/qPmObB8HHr01taNISrKormZBdA1Uy5ttGIzjyPNZv1URG275ub1fBjYB3l4JOXniV4yHrmYBWlD4amFI5xuzJ8+9ey6qpuP0bI5UachyqgsP/vpiMDA7Do3Fv2pLgSnlPOTgohQxjI3l0ZzlmNmMIDfjL6Iwg5QZ2yXqhGcfIeUZLkGMov7my6xcxbXFm+q6f4x7GTHErxK69BKa+zI7l+OUepNhIM9ub8ok72ILWU6zmrli9O3Z4x3mLUC1uzfs147QgMsMcqM5j8zxXeczcujEXXyd+i655QUPB6Jj1WPwnpxry5X8Sp8FvxYNVt+3TQryhsGqZS4dn35tJUbhbubZU84F5Bg8DGs15cP5jzzNym3Rsa6Bj650rZh3V5mLz2+VfZbwaK5rVqe9d9xZxebZ+oyndwqY0LsNGAh+dLbQbS+qFHJ9CfHC7CPypRHyTkBhFqTa/Zo1eMLQcyB+pIQ6C1+KfEK5VApaPJnO/JCFiDz/PqjJeA9sUw9sixMRCMS/QPMPPowy6Rh+WbxjHxaOvKNJdSpO97RdBw6V3OkdmEWJ1p0WvvesHF/ZTN5PkTZLkvhncy+bcPL67c6BcI5w1cc4nGFLRWBCrL60Rgj1AuKbaNgn8DavYqP4/BLUt2wcJ/qehPC1V/PxiH40WnM5sO0fD0fu/YDUIS1q/U+uEiSoeT+sf5L49SkWTEuM+zhYEJL6mcTM0MzJGrnxqoMnPsXFVYy+ksjIn6ZaU8Z/rHAcGvBWKfYXqrLzxtrB9OSsKjXhXZEEoNDjhpOYhmp+PPvRVyYbQgrOk+P4SpNSvPb7H8AmZ0O0R+M4KEU7Uqeep5kPxAbetXHwPqN26rZeb1gGIFs/3bxyFXLBVeszIfgfuyDcIHlvVIsCA89H+W387kf8Iwh0p7botjBK5BVwUKEm8ZO/xh6pNcbpqtGEQBW6pVlou2L6LFX6Wu5qHiRl+rZbJPIKZer2UnJlHx8XTWqfWMELwlHnXktqKBOnR84MYBZn823sx1X5OORoW/AZFWUaiwczJZf4GB0aDFX9knYrAv/55D/P0+WL+nksXkYjGpqMennc2sj9V//d5yxCEhLDEpRmAAohTJ9XHJH/AAnY12PoNgZ72ZPbQxAfLDbBLTAsMQDTHh1lROxFdtB2Jc6giUBe/c/CyRBN9bH+oe6TD7B1WuhJeZ63cLm+p7So1jSmtbolhrEuTtGxP/x/8/PbR7iALRrW7XN1BAsf8ff3Okl3KS12hQMBy0TqOSqc+UmZtJSSwFY7sO5d288BDRLHwe//PLczPNhw57P4Hzu/GZlUQq7Lc9DxeneoB7SwKDz4mGA8fbJQ9I+2L812MJThM63C4W1Ota+GPRKcKY1ZwBfj5yctTRQLh2KnDHu/WBy8E8LWlPMKJcMg6tL/0OX4nObrf5l8iYtdJoEh/AjOl+z09vQnFHeGbt1+BBuO6wO7eHMwxZmgqzLGLDGOpYVccmEw52femL+z6NwKJ8jdGVNRHoER7VOKc/hoAf3g9vBUXgZsCs2NSucSx25kizFkbg/wCk39hsAQAAAAAAAAAAgAAAAAAAABgbAAAAAAAAJQIAAAAAAAA=eF592DuuHCEARcG3M2/ZSyAk7KADAgLUQgghhFiCIxLs8smmZrKr6Q8/P//r969/e4BH+AN/4Qme4QX+wSu8wTt8wCd8wTf8dO9wPt8e4BH+wF94gmd4gX/wCm/wDh/wCV/wDT/d/4Pzu9sDPMIf+AtP8Awv8A9e4Q3e4QM+4Qu+4af7OnR2uT3AI/yBv/AEz/AC/+AV3uAdPuATvuAbfrrvA2eX2wM8wh/4C0/wDC/wD17hDd7hAz7hC77hp/s+fHa5PcAj/IG/8ATP8AL/4BXe4B0+4BO+4Bt+up+Dzi63B3iEP/AXnuAZXuAfvMIbvMMHfMIXfMNP93Po2eX2AI/wB/7CEzzDC/yDV3iDd/iAT/iCb/jpfg84u9we4BH+wF94gmd4gX/wCm/wDh/wCV/wDT/d72Fnl9sDPMIf+AtP8Awv8A9e4Q3e4QM+4Qu+4af7PfjscnuAR/gDf+EJnuEF/sErvME7fMAnfME3/HSfQ5xdbg/wCH/gLzzBM7zAP3iFN3iHD/iEL/iGn/rlZ5fbAzzCH/gLT/AML/APXuEN3uEDPuELvuGn+xzu7HJ7gEf4A3/hCZ7hBf7BK7zBO3zAJ3zBN/x0n4OeXW4P8Ah/4C88wTO8wD94hTd4hw/4hC/4hp/uc+izy+0BHuEP/IUneIYX+Aev8Abv8AGf8AXf8Pv709nl9gCP8Af+whM8wwv8g1d4g3f4gE/4gu+//A9Uo1hvAQAAAAAAAAAAgAAAAAAAAAAgAAAAAAAACQUAAAAAAAA=eF5d2VcY/3MZh3Ffe+89/ypEk7IrikKEQjtUSqRCaEfTqow0kNlCSbSXaNOwtSM7FCWrIgc+r4Pf/Zy8r+u+7rP78JlrrkdvGrv02KXC5x67TDx8nrHLxsPnHbtcPHy+scvHw+cfu0I8fIGxK8bDFxy7Ujx8obErx8MXHrtKPHyRsavGwxcdu1o8fLGxq8fDFx+7Rjx8ibFz4uFLjl0zns56PnbsY8L1fFw8XM+14uF6rh0P13OdeLiej4+H67luPFzP9eLhej4hHq7nE+Phej4pHq7nk+Phej4lHq7nU+Phc8auHw/XfYN4Ouv59LFPC9dzw3i4nhvFw/XcOB6u5ybxcD03jYfruVk8XM/N4+F6PiMerucz4+F6PiserucW8XA9t4yH6/nseLiez4mH675VPJ31fO7YrcP1fF48XM9t4uF6bhsP13O7eLiez4+H67l9PFzPHeLher4gHq7njvFwPXeKh+u5czx8y7EvjIfr+aJ4uJ67xMN13zWeznq+eOxu4Xq+JB6u50vj4Xq+LB6u58vj4Xq+Ih6u5yvj4Xq+Kh6u5+7xcD33iIfruWc8XM9Xx8P1fE08XM/XxsP13Cservvr4ums595jXx+u5xvi4XruEw/Xc994uJ5vjIfruV88XM83xcP1fHM8XM+3xMP13D8evufYA+Lheh4YD9fzrfFwPQ+Kh+t5cDxc90Pi6azn28e+LVzPd8TD9XxnPFzPd8XD9Xx3PFzP98TD9XxvPFzPQ+Pheh4WD9fzffFwPd8fD9fzA/FwPT8YD9fzQ/FwPT8cD9f98Hg663nk2CPC9TwqHq7n0fFwPT8SD9fzo/FwPT8WD9fzmHi4nsfGww8be1w8XM/j4+F6fjwerucJ8XA9PxEP1/OT8XA9PxUP1/3T8XTW86SxJ4breXI8XM/PxMP1PCUeruep8XA9T4uH63l6PFzPM+Lhep4ZD9fzs/FwPT8XD9fz8/FwPb8QD9fzi/FwPc+Kh+t+djyd9fzS2HPC9fxyPFzPc+Phen4lHq7nefFwPb8aDz997PnxcD0viIfr+bV4uJ5fj4fr+Y14uJ7fjIfr+a14uJ7fjofr+Z14uO7fjaeznt8f+71wPX8QD9fzwni4nj+Mh+t5UTxcz4vj4Xr+KB6u54/j4Xr+JB6u50/j4Xr+LB6u58/j4Xr+Ih6u5yXxcD0vjYfr/st4Ouv567G/CtfzN/FwPS+Lh+t5eTz8orFXxMP1vDIerudV8XA9r46H63lNPFzPa+Phev42Hq7n7+Lhev4+Hq7nH+Lhev4xHq77n+LprOdfxv45XM/r4uF6Xh8P1/Ov8XA9b4iH63ljPFzPm+Lhet4cD9fzlni4nrfGw/W8LR6u59/i4XreHg/X8454uJ53xsN1/3s8nfW8a+w/wvW8Ox5+/dh/xsP1/Fc8XM974uF6/jserue98XA974uH63l/PFzPB+Lhej4YD9fzP/FwPf8bD9fzf/FwPR+Kh+v+cDyd9fTw+3+4npOH4DTL9Zx7mvVwPeeZZj1cz3mnWQ/Xc75p1sP1nH+a9XA9F5hmPVzPBadZD9dzoWnWw/VceJr1cD0XmWY9XM9Fp1kP13OxadbD9Vx8mvVw3ZeYZj2d9Vxq7JLprOfS8XA9l4mH67lsPFzP5eLhei4fD9dzhXi4nivGw/VcKR6u58rxcD1XiYfruWo8XM/V4uF6rh4P13ONeLjuc+I9AtvA4m4=AQAAAAAAAAAAgAAAAAAAAAAIAAAAAAAAkAEAAAAAAAA=eF4txRGAGgAAAMC2BUEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQBEEQBA8PDw9BEARBEARBEARBEARBEARBMOhOLhj4CDnsiKOOOe6Ek0457YyzzjnvgosuueyKq6657oabbrntjrvuue+Bhx557Imn/vK3f/zrmedeeOmV19546533Pvjok8+++Oqb73746ZffDvz5FHTIYUccdcxxJ5x0ymlnnHXOeRdcdMllV1x1zXU33HTLbXfcdc99Dzz0yGNPPPWXv/3jX88898JLr7z2xlvvvPfBR5989sVX33z3w0+//Hbg76egQw474qhjjjvhpFNOO+Osc8674KJLLrviqmuuu+GmW26746577nvgoUcee+Kpv/ztH/965rkXXnrltTfeeue9Dz765LMvvvrmux9++uW3A/8+BR1y2BFHHXPcCSedctoZZ51z3gUXXXLZFVddc90NN91y2x133XPfAw898tgTT/3lb//41zPPvfDSK6+98dY7733w0SefffHVN9/98NMvvx0IfvoP9CZ/hQ==AQAAAAAAAAAAgAAAAAAAAAABAAAAAAAADAAAAAAAAAA=eF7j5BzZAACFvAkB
-  </AppendedData>
-</VTKFile>

ogstools/studies/convergence/examples/steady_state_diffusion.prj

+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+    <meshes>
+        <mesh>square_domain.vtu</mesh>
+        <mesh>square_physical_group_left.vtu</mesh>
+        <mesh>square_physical_group_bottom.vtu</mesh>
+        <mesh>square_physical_group_right.vtu</mesh>
+    </meshes>
+    <processes>
+        <process>
+            <name>SteadyStateDiffusion</name>
+            <type>STEADY_STATE_DIFFUSION</type>
+            <integration_order>2</integration_order>
+            <process_variables>
+                <process_variable>pressure</process_variable>
+            </process_variables>
+            <secondary_variables>
+                <secondary_variable internal_name="darcy_velocity" output_name="velocity"/>
+            </secondary_variables>
+        </process>
+    </processes>
+    <media>
+        <medium id="0">
+            <phases/>
+            <properties>
+                <property>
+                    <name>diffusion</name>
+                    <type>Constant</type>
+                    <value>1</value>
+                </property>
+                <property>
+                    <name>reference_temperature</name>
+                    <type>Constant</type>
+                    <value>293.15</value>
+                </property>
+            </properties>
+        </medium>
+    </media>
+    <time_loop>
+        <processes>
+            <process ref="SteadyStateDiffusion">
+                <nonlinear_solver>basic_picard</nonlinear_solver>
+                <convergence_criterion>
+                    <type>DeltaX</type>
+                    <norm_type>NORM2</norm_type>
+                    <abstol>1.e-6</abstol>
+                </convergence_criterion>
+                <time_discretization>
+                    <type>BackwardEuler</type>
+                </time_discretization>
+                <time_stepping>
+                    <type>SingleStep</type>
+                </time_stepping>
+            </process>
+        </processes>
+        <output>
+            <type>VTK</type>
+            <prefix>steady_state_diffusion</prefix>
+            <variables>
+                <variable> pressure </variable>
+                <variable> velocity      </variable>
+            </variables>
+            <suffix>_ts_{:timestep}_t_{:time}</suffix>
+        </output>
+    </time_loop>
+    <parameters>
+        <parameter>
+            <name>p0</name>
+            <type>Constant</type>
+            <value>0</value>
+        </parameter>
+        <parameter>
+            <name>p_neumann</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+        <parameter>
+            <name>p_Dirichlet</name>
+            <type>Constant</type>
+            <value>1</value>
+        </parameter>
+    </parameters>
+    <process_variables>
+        <process_variable>
+            <name>pressure</name>
+            <components>1</components>
+            <order>1</order>
+            <initial_condition>p0</initial_condition>
+            <boundary_conditions>
+                <boundary_condition>
+                    <mesh>square_physical_group_left</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>p_Dirichlet</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <mesh>square_physical_group_bottom</mesh>
+                    <type>Dirichlet</type>
+                    <parameter>p_Dirichlet</parameter>
+                </boundary_condition>
+                <boundary_condition>
+                    <mesh>square_physical_group_right</mesh>
+                    <type>Neumann</type>
+                    <parameter>p_neumann</parameter>
+                </boundary_condition>
+            </boundary_conditions>
+        </process_variable>
+    </process_variables>
+    <nonlinear_solvers>
+        <nonlinear_solver>
+            <name>basic_picard</name>
+            <type>Picard</type>
+            <max_iter>10</max_iter>
+            <linear_solver>general_linear_solver</linear_solver>
+        </nonlinear_solver>
+    </nonlinear_solvers>
+    <linear_solvers>
+        <linear_solver>
+            <name>general_linear_solver</name>
+            <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+            <eigen>
+                <solver_type>CG</solver_type>
+                <precon_type>DIAGONAL</precon_type>
+                <max_iteration_step>10000</max_iteration_step>
+                <error_tolerance>1e-16</error_tolerance>
+            </eigen>
+            <petsc>
+                <prefix>gw</prefix>
+                <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+            </petsc>
+        </linear_solver>
+    </linear_solvers>
+</OpenGeoSysProject>

ogstools/studies/convergence/examples/steady_state_diffusion.py

+from copy import deepcopy
+from pathlib import Path
+from typing import Union
+
+import numpy as np
+import pyvista as pv
+
+
+def _c_k(k):
+    return 0.5 * (2 * k - 1) * np.pi
+
+
+def _a_k(k):
+    return 2 / (_c_k(k) ** 2 * np.cosh(_c_k(k)))
+
+
+def _h(points):
+    result = np.ones(len(points))
+    for k in np.arange(1, 100):
+        c_k_val = _c_k(k)
+        sin_c_k = np.sin(c_k_val * points[:, 1])
+        sinh_c_k = np.sinh(c_k_val * points[:, 0])
+        result += _a_k(k) * sin_c_k * sinh_c_k
+    return result
+
+
+def analytical_solution(topology: Union[Path, pv.DataSet]) -> pv.DataSet:
+    mesh = topology if isinstance(topology, pv.DataSet) else pv.read(topology)
+    new_mesh = deepcopy(mesh)
+    new_mesh.clear_point_data()
+    points = new_mesh.points
+    new_mesh.point_data["pressure"] = _h(points)
+    return new_mesh

ogstools/studies/convergence/generate_report.py → ogstools/studies/convergence/study.py

 import argparse
-import tempfile
 from pathlib import Path
+from typing import Optional
 
-import jupytext
-import papermill as pm
+from ogstools.definitions import ROOT_DIR
+from ogstools.workflow import jupytext_to_jupyter
 
 
-def execute_convergence_study(
-    output_path: str, mesh_paths: list[str], property_name: str, timestep: int
+def run_convergence_study(
+    output_name: Path,
+    mesh_paths: list[Path],
+    topology_path: Path,
+    property_name: str,
+    timestep: int = 0,
+    refinement_ratio: Optional[float] = None,
+    reference_solution_path: Optional[Path] = None,
+    prepare_only: bool = False,
+    progress_bar: bool = False,
 ) -> None:
     params = {
-        "mesh_paths": mesh_paths,
+        "mesh_paths": [str(mesh_path) for mesh_path in mesh_paths],
+        "topology_path": str(topology_path),
         "property_name": property_name,
         "timestep": timestep,
+        "refinement_ratio": refinement_ratio,
+        "reference_solution_path": reference_solution_path,
     }
-    parent = Path(__file__).resolve().parent
-    template_path = str(parent) + "/convergence_study_template.md"
-    nb = jupytext.read(template_path)
-    with tempfile.NamedTemporaryFile(delete=False, suffix=".ipynb") as temp:
-        jupytext.write(nb, temp.name, fmt="py:percent")
-    pm.execute_notebook(
-        input_path=temp.name,
-        output_path=output_path,
-        parameters=params,
+    template_path = Path(
+        str(ROOT_DIR) + "/studies/templates/convergence_study.py"
+    )
+    jupytext_to_jupyter(
+        template_path,
+        output_name,
+        params=params,
+        prepare_only=prepare_only,
+        progress_bar=progress_bar,
     )
-    Path(temp.name).unlink()
 
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser(description="Convergence Study")
     parser.add_argument("output", help="Path to the output notebook")
     parser.add_argument("--mesh_paths", nargs="+", help="List of mesh paths")
+    parser.add_argument("--topology_path", help="Path to topology mesh")
     parser.add_argument("--property_name", help="Name of the property")
     parser.add_argument("--timestep", help="Timestep to read")
 
     args = parser.parse_args()
 
-    execute_convergence_study(
-        args.output, args.mesh_paths, args.property_name, args.timestep
+    run_convergence_study(
+        args.output,
+        args.mesh_paths,
+        args.topology_path,
+        args.property_name,
+        args.timestep,
     )