Convergence study (spatial refinement)

This example shows one possible implementation of how to do a convergence study. It uses the project file from the following benchmark with multiple discretizations to evaluate the accuracy of the numerical solutions. ogs: elliptic neumann benchmark

Here is some theoretical background for the topic of grid convergence:

Nasa convergence reference

More comprehensive reference

At least three meshes of increasing refinement are required for the convergence study. The three finest meshes are used to calculated the Richardson extrapolation. The third coarsest mesh will be used for the topology to evaluate the results. Its nodes should be shared by the finer meshes, otherwise interpolation will influence the results. With unstructured grids this can be achieved as well with refinement by splitting.

The results to analyze are generated on the fly with the following code. If you are only interested in the convergence study, please skip to Hydraulic pressure convergence.

First, the required packages are imported and an output directory is created:

from pathlib import Path
from tempfile import mkdtemp

from IPython.display import HTML
from ogs6py import ogs

from ogstools import meshlib, meshplotlib, msh2vtu, propertylib, workflow
from ogstools.studies import convergence
from ogstools.studies.convergence.examples import (
    steady_state_diffusion_analytical_solution,
)

meshplotlib.setup.reset()
temp_dir = Path(mkdtemp(suffix="steady_state_diffusion"))
report_name = str(temp_dir / "report.ipynb")
result_paths = []

The meshes and their boundaries are generated easily via gmsh and ogstools.msh2vtu. Then we run the different simulations with increasingly fine spatial discretization via ogs6py and store the results for the convergence study.

refinements = 6
edge_cells = [2**i for i in range(refinements)]
for n_edge_cells in edge_cells:
    msh_path = temp_dir / "square.msh"
    meshlib.gmsh_meshing.rect(
        n_edge_cells=n_edge_cells, structured_grid=True, out_name=msh_path
    )
    msh2vtu.msh2vtu(filename=msh_path, output_path=temp_dir, log_level="ERROR")

    model = ogs.OGS(
        PROJECT_FILE=temp_dir / "default.prj",
        INPUT_FILE=convergence.examples.steady_state_diffusion_prj,
    )
    prefix = "steady_state_diffusion_" + str(n_edge_cells)
    model.replace_text(prefix, ".//prefix")
    model.write_input()
    ogs_args = f"-m {temp_dir} -o {temp_dir}"
    model.run_model(write_logs=False, args=ogs_args)
    result_paths += [str(temp_dir / (prefix + ".pvd"))]

OGS finished with project file /tmp/tmphul2v_1usteady_state_diffusion/default.prj.
Execution took 0.1162104606628418 s

OGS finished with project file /tmp/tmphul2v_1usteady_state_diffusion/default.prj.
Execution took 0.1116323471069336 s

OGS finished with project file /tmp/tmphul2v_1usteady_state_diffusion/default.prj.
Execution took 0.11133050918579102 s

OGS finished with project file /tmp/tmphul2v_1usteady_state_diffusion/default.prj.
Execution took 0.11393499374389648 s

OGS finished with project file /tmp/tmphul2v_1usteady_state_diffusion/default.prj.
Execution took 0.11359739303588867 s

OGS finished with project file /tmp/tmphul2v_1usteady_state_diffusion/default.prj.
Execution took 0.12082672119140625 s

Here we calculate the analytical solution on one of the meshes:

analytical_solution_path = temp_dir / "analytical_solution.vtu"
solution = steady_state_diffusion_analytical_solution(
    meshlib.MeshSeries(result_paths[-1]).read(0)
)
meshplotlib.setup.show_element_edges = True
fig = meshplotlib.plot(solution, propertylib.presets.hydraulic_height)
solution.save(analytical_solution_path)

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.

convergence.run_convergence_study(
    output_name=report_name,
    mesh_paths=result_paths,
    property_name="hydraulic_height",
    timevalue=1,
    refinement_ratio=2.0,
    reference_solution_path=str(analytical_solution_path),
)
HTML(workflow.jupyter_to_html(report_name, show_input=False))

Grid convergence¶

If the shown values are approximately 1, this means that the results are in asymptotic range of convergence.

Grid comparison¶

Visualizing the requested mesh property on the 3 finest discretizations:

Richardson extrapolation.¶

Visualizing the Richardson extrapolation of the requested mesh property. If a reference solution is provided, the difference between the two is shown as well. Otherwise the difference between the finest discretization and the Richardson extrapolation is shown.

Convergence metrics¶

mean element length	maximum	minimum	abs. error (max)	abs. error (min)	abs. error (L2 norm)	rel. error (max)	rel. error (min)	rel. error (L2 norm)
1	1.75	1	0.074685	0	0.58656	0.04458	0	0.053572
0.5	1.6857	1	0.0104	0	0.19746	0.0062075	0	0.018034
0.25	1.678	1	0.0026531	0	0.06575	0.0015836	0	0.0060051
0.125	1.676	1	0.00065028	0	0.0068033	0.00038815	0	0.00062136
0.0625	1.6755	1	0.00016145	0	0.0024512	9.6371e-05	0	0.00022387
0.03125	1.6754	1	4.0085e-05	0	0.002059	2.3927e-05	0	0.00018805
0	1.6753	1	0	0	9.1551e-16	0	0	8.3616e-17

Relative errors¶

Absolute values¶

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. The code cells from the templated notebook are show here for transparency.

convergence.run_convergence_study(
    output_name=report_name,
    mesh_paths=result_paths,
    property_name="velocity",
    timevalue=1,
    refinement_ratio=2.0,
)
HTML(workflow.jupyter_to_html(report_name, show_input=True))

# Import required modules and customize plot settings.

import numpy as np  # noqa: E402
import pyvista as pv  # noqa: E402

from ogstools import meshlib, meshplotlib, propertylib, studies  # noqa: E402

meshplotlib.setup.reset()
meshplotlib.setup.show_element_edges = True
meshplotlib.setup.fig_scale = 0.5
meshplotlib.setup.combined_colorbar = False

# Here, the meshes are read, a Property object is created from the property
# name and the Richardson extrapolation calculated.
# The 3 finest meshes of those provided will be used for the Richardson
# extrapolation.

mesh_series = [meshlib.MeshSeries(mesh_path) for mesh_path in mesh_paths]
timestep_sizes = [np.mean(np.diff(ms.timevalues)) for ms in mesh_series]
meshes = [ms.read_closest(timevalue) for ms in mesh_series]
topology: pv.UnstructuredGrid = meshes[-3]
mesh_property = propertylib.presets.get_preset(property_name, meshes[0])
richardson = studies.convergence.richardson_extrapolation(
    meshes, mesh_property, topology, refinement_ratio
)

Grid convergence¶

If the shown values are approximately 1, this means that the results are in asymptotic range of convergence.

fig = meshplotlib.plot(richardson, "grid_convergence")

Grid comparison¶

Visualizing the requested mesh property on the 3 finest discretizations:

fig = meshplotlib.plot(meshes[-3:], mesh_property)

Richardson extrapolation.¶

Visualizing the Richardson extrapolation of the requested mesh property. If a reference solution is provided, the difference between the two is shown as well. Otherwise the difference between the finest discretization and the Richardson extrapolation is shown.

fig = meshplotlib.plot(richardson, mesh_property)

data_key = mesh_property.data_name
if reference_solution_path is None:
    diff_mesh = meshlib.difference(
        richardson, topology.sample(meshes[-1]), mesh_property
    )
    fig = meshplotlib.plot(diff_mesh, mesh_property)
else:
    reference_solution = topology.sample(
        meshlib.MeshSeries(reference_solution_path).read_closest(timevalue)
    )
    diff_mesh = meshlib.difference(
        reference_solution, richardson, mesh_property
    )
    fig = meshplotlib.plot(diff_mesh, mesh_property)

Convergence metrics¶

metrics = studies.convergence.convergence_metrics(
    meshes, richardson, mesh_property, timestep_sizes
)
metrics.style.format("{:,.5g}").hide()

mean element length	maximum	minimum	abs. error (max)	abs. error (min)	abs. error (L2 norm)	rel. error (max)	rel. error (min)	rel. error (L2 norm)
1	1.0607	1.1938e-16	-333.69	-5.4421e-15	334.01	0.99683	0.97853	0.99759
0.5	1.1571	3.0629e-13	-333.59	3.0073e-13	333.6	0.99654	54.074	0.99637
0.25	1.5878	7.8016e-14	-333.16	7.2455e-14	333.16	0.99526	13.028	0.99508
0.125	2.0261	1.7428e-14	-332.73	1.1866e-14	332.73	0.99395	2.1337	0.99377
0.0625	2.4666	6.1454e-15	-332.29	5.8397e-16	332.29	0.99263	0.105	0.99245
0.03125	2.9076	1.4762e-15	-331.84	-4.0853e-15	331.84	0.99131	0.73457	0.99113
0	334.75	5.5615e-15	0	0	5.4177e-16	0	0	1.6181e-18

Relative errors¶

meshplotlib.core.plt.rcdefaults()
fig = studies.convergence.plot_convergence_errors(metrics)

Absolute values¶

fig = studies.convergence.plot_convergence(metrics, mesh_property)