.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/howto_postprocessing/plot_convergence_study_nuclear_decay.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_howto_postprocessing_plot_convergence_study_nuclear_decay.py: Spatial & temporal refinement - nuclear decay ============================================= This example shows one possible implementation of how to do a convergence study with spatial and 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 :py:mod:`ogstools.physics.nuclearwasteheat` model. Here is some theoretical background for the topic of grid convergence: `Nasa convergence reference `_ `More comprehensive reference `_ At least three meshes from simulations of increasing refinement are required for the convergence study. The third finest mesh is chosen per default as the topology to evaluate the results on. 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 `Temperature convergence at maximum heat production (t=30 yrs)`_. First, the required packages are imported and a temporary output directory is created: .. GENERATED FROM PYTHON SOURCE LINES 32-45 .. code-block:: Python from pathlib import Path from tempfile import mkdtemp import matplotlib.pyplot as plt import numpy as np from IPython.display import HTML from scipy.constants import Julian_year as sec_per_yr import ogstools as ogs from ogstools import examples, meshlib, msh2vtu, physics, studies, workflow temp_dir = Path(mkdtemp(prefix="nuclear_decay")) .. GENERATED FROM PYTHON SOURCE LINES 46-49 Let's run the different simulations with increasingly fine spatial and temporal discretization via ogs6py. The mesh and its boundaries are generated easily via gmsh and :py:mod:`ogstools.msh2vtu`. First some definitions: .. GENERATED FROM PYTHON SOURCE LINES 51-61 .. code-block:: Python n_refinements = 4 time_step_sizes = [30.0 / (2.0**r) for r in range(n_refinements)] prefix = "stepsize_{0}" sim_results = [] msh_path = temp_dir / "rect.msh" script_path = Path(examples.pybc_nuclear_decay).parent prj_path = examples.prj_nuclear_decay ogs_args = f"-m {temp_dir} -o {temp_dir} -s {script_path}" edge_cells = [5 * 2**i for i in range(n_refinements)] .. GENERATED FROM PYTHON SOURCE LINES 62-63 Now the actual simulations: .. GENERATED FROM PYTHON SOURCE LINES 65-78 .. code-block:: Python for dt, n_cells in zip(time_step_sizes, edge_cells, strict=False): ogs.meshlib.rect( lengths=100.0, n_edge_cells=[n_cells, 1], out_name=msh_path ) _ = msh2vtu.msh2vtu(msh_path, output_path=temp_dir, log_level="ERROR") prj = ogs.Project(output_file=temp_dir / "default.prj", input_file=prj_path) prj.replace_text(str(dt * sec_per_yr), ".//delta_t") prj.replace_text(prefix.format(dt), ".//prefix") prj.write_input() prj.run_model(write_logs=False, args=ogs_args) sim_results += [temp_dir / (prefix.format(dt) + "_rect_domain.xdmf")] .. rst-class:: sphx-glr-script-out .. code-block:: none OGS finished with project file /tmp/nuclear_decaycp0tgweo/default.prj. Execution took 1.7088615894317627 s Project file written to output. OGS finished with project file /tmp/nuclear_decaycp0tgweo/default.prj. Execution took 1.7861390113830566 s Project file written to output. OGS finished with project file /tmp/nuclear_decaycp0tgweo/default.prj. Execution took 1.9428648948669434 s Project file written to output. OGS finished with project file /tmp/nuclear_decaycp0tgweo/default.prj. Execution took 2.2398030757904053 s Project file written to output. .. GENERATED FROM PYTHON SOURCE LINES 79-81 Let's extract the temperature evolution and the applied heat via vtuIO and plot both: .. GENERATED FROM PYTHON SOURCE LINES 83-111 .. code-block:: Python time = np.append(0.0, np.geomspace(1.0, 180.0, num=100)) repo = physics.nuclearwasteheat.repo_2020_conservative heat = repo.heat(time, time_unit="yrs", power_unit="kW") fig, (ax1, ax2) = plt.subplots(figsize=(8, 8), nrows=2, sharex=True) ax2.plot(time, heat, lw=2, label="reference", color="k") for sim_result, dt in zip(sim_results, time_step_sizes, strict=False): mesh_series = ogs.MeshSeries(sim_result) results = {"heat_flux": [], "temperature": []} for ts in mesh_series.timesteps: mesh = mesh_series.mesh(ts) results["temperature"] += [np.max(mesh.point_data["temperature"])] max_T = ogs.variables.temperature.transform(results["temperature"]) # times 2 due to symmetry, area of repo, to kW results["heat_flux"] += [np.max(mesh.point_data["heat_flux"][:, 0])] tv = np.asarray(mesh_series.timevalues("a")) ax1.plot(tv, max_T, lw=1.5, label=f"{dt=}") edges = np.append(0, tv) mean_t = 0.5 * (edges[1:] + edges[:-1]) applied_heat = repo.heat(mean_t, time_unit="yrs", power_unit="kW") ax2.stairs(applied_heat, edges, lw=1.5, label=f"{dt=}", baseline=None) ax2.set_xlabel("time / yrs") ax1.set_ylabel("max T / °C") ax2.set_ylabel("heat / kW") ax1.legend() ax2.legend() fig.show() .. image-sg:: /auto_examples/howto_postprocessing/images/sphx_glr_plot_convergence_study_nuclear_decay_001.png :alt: plot convergence study nuclear decay :srcset: /auto_examples/howto_postprocessing/images/sphx_glr_plot_convergence_study_nuclear_decay_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 112-122 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 range 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 quadratically, as expected. .. GENERATED FROM PYTHON SOURCE LINES 124-134 .. code-block:: Python report_name = temp_dir / "report.ipynb" studies.convergence.run_convergence_study( output_name=report_name, mesh_paths=sim_results, timevalue=30 * sec_per_yr, variable_name="temperature", refinement_ratio=2.0, ) HTML(workflow.jupyter_to_html(report_name, show_input=False)) .. raw:: html
report


.. GENERATED FROM PYTHON SOURCE LINES 135-141 Temperature convergence at maximum temperature (t=150 yrs) ---------------------------------------------------------- The temperature convergence at this timevalue is much closer to 1, indicating a better convergence behaviour, which is due to the temperature gradient now changing only slowly. Convergence order is again quadratic. .. GENERATED FROM PYTHON SOURCE LINES 143-152 .. code-block:: Python studies.convergence.run_convergence_study( output_name=report_name, mesh_paths=sim_results, timevalue=150 * sec_per_yr, variable_name="temperature", refinement_ratio=2.0, ) HTML(workflow.jupyter_to_html(report_name, show_input=False)) .. raw:: html
report


.. GENERATED FROM PYTHON SOURCE LINES 153-160 Convergence evolution over all timesteps ---------------------------------------------------------- We can also run the convergence evaluation on all timesteps and look at the relative errors (between finest discretization and Richardson extrapolation) and the convergence order over time to get a better picture of the transient model behavior. .. GENERATED FROM PYTHON SOURCE LINES 162-167 .. code-block:: Python mesh_series = [meshlib.MeshSeries(sim_result) for sim_result in sim_results] evolution_metrics = studies.convergence.convergence_metrics_evolution( mesh_series, ogs.variables.temperature, units=["s", "yrs"] ) .. rst-class:: sphx-glr-script-out .. code-block:: none 0%| | 0/7 [00:00` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_convergence_study_nuclear_decay.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_convergence_study_nuclear_decay.zip `