Note

Go to the end to download the full example code.

Spatial refinement - steady state diffusion#

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

import ogstools as ogs
from ogstools import examples, msh2vtu, variables, workflow
from ogstools.studies import convergence

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"
    ogs.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.Project(
        output_file=temp_dir / "default.prj",
        input_file=examples.prj_steady_state_diffusion,
    )
    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/tmppnsa7lz4steady_state_diffusion/default.prj.
Execution took 0.13386845588684082 s
Project file written to output.

OGS finished with project file /tmp/tmppnsa7lz4steady_state_diffusion/default.prj.
Execution took 0.14445233345031738 s
Project file written to output.

OGS finished with project file /tmp/tmppnsa7lz4steady_state_diffusion/default.prj.
Execution took 0.132127046585083 s
Project file written to output.

OGS finished with project file /tmp/tmppnsa7lz4steady_state_diffusion/default.prj.
Execution took 0.14403414726257324 s
Project file written to output.

OGS finished with project file /tmp/tmppnsa7lz4steady_state_diffusion/default.prj.
Execution took 0.14754223823547363 s
Project file written to output.

OGS finished with project file /tmp/tmppnsa7lz4steady_state_diffusion/default.prj.
Execution took 0.14854931831359863 s
Project file written to output.

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

analytical_solution_path = temp_dir / "analytical_solution.vtu"
solution = examples.analytical_diffusion(
    ogs.MeshSeries(result_paths[-1]).mesh(0)
)
ogs.plot.setup.show_element_edges = True
fig = ogs.plot.contourf(solution, variables.hydraulic_head)
solution.save(analytical_solution_path)
plot convergence study steady state diffusion

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,
    variable_name="hydraulic_head",
    timevalue=1,
    refinement_ratio=2.0,
    reference_solution_path=str(analytical_solution_path),
)
HTML(workflow.jupyter_to_html(report_name, show_input=False))
report