.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/howto_plot/plot_analytical_solution.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_plot_plot_analytical_solution.py: How to compare results with reference data or an analytical solution ==================================================================== .. sectionauthor:: Florian Zill (Helmholtz Centre for Environmental Research GmbH - UFZ) This example provides recipes for frequently needed plots in the OpenGeoSys benchmark section. It is a repeating pattern, to compare numerical results with analytical or reference data and evaluate the errors. The goal here is to have a standardized recipe for cleaner code and less repetition. .. GENERATED FROM PYTHON SOURCE LINES 14-24 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import ogstools as ot from ogstools import examples temp = ot.variables.temperature .. GENERATED FROM PYTHON SOURCE LINES 30-34 Compute the analytical solution ------------------------------- Note, that the relative errors are calculated with the reference values in data units, i.e. in Kelvin. .. GENERATED FROM PYTHON SOURCE LINES 36-47 .. code-block:: Python results = examples.load_meshseries_diffusion_3D() x = results[0].points[:, 0] analytical_func = examples.anasol.heat_conduction_temperature ref_values = np.asarray([analytical_func(x, tv) for tv in results.timevalues]) abs_error = ref_values - results["temperature"] rel_error = abs_error / ref_values np.testing.assert_array_less(np.abs(abs_error), 6) np.testing.assert_array_less(np.abs(rel_error), 0.02) results = results.scale(time=("s", "h")) .. GENERATED FROM PYTHON SOURCE LINES 48-49 Write the data into the results, to leverage plotting features. .. GENERATED FROM PYTHON SOURCE LINES 51-55 .. code-block:: Python results.point_data[temp.anasol.data_name] = ref_values results.point_data[temp.abs_error.data_name] = abs_error results.point_data[temp.rel_error.data_name] = rel_error .. GENERATED FROM PYTHON SOURCE LINES 56-63 Comparing data over a line for multiple timesteps ------------------------------------------------- If the resulting data is not already one-dimensional we need to extract the points which make up the line of interest. We can do so, by using ``extract`` and the point indices or by using ``extract_probe`` and a set of points. We need to further select the timesteps we want to plot. This is done by indexing the MeshSeries object. .. GENERATED FROM PYTHON SOURCE LINES 65-78 .. code-block:: Python fig, axs = plt.subplots(1, 3, figsize=[40, 10]) # extract every second timestep and only the points on the x-axis x_edge = results[::2].extract( (results[0].points[:, 1] == 0) & (results[0].points[:, 2] == 0) ) labels = [f"{tv:.1f} h" for tv in x_edge.timevalues] axs[0].plot([], [], "--k", label="analytical\nsolution") ot.plot.line(x_edge, temp, ax=axs[0], marker="o", labels=labels) ot.plot.line(x_edge, temp.anasol, ax=axs[0], ls="--") ot.plot.line(x_edge, temp.abs_error, ax=axs[1]) ot.plot.line(x_edge, temp.rel_error, ax=axs[2]) fig.tight_layout() .. image-sg:: /auto_examples/howto_plot/images/sphx_glr_plot_analytical_solution_001.png :alt: plot analytical solution :srcset: /auto_examples/howto_plot/images/sphx_glr_plot_analytical_solution_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 79-83 Comparing data of multiple points over time ------------------------------------------- Here, we use the same strategy as before, but in the ``plot.line`` function, we specify, that we want to plot over the time dimension. .. GENERATED FROM PYTHON SOURCE LINES 85-96 .. code-block:: Python fig, axs = plt.subplots(1, 3, figsize=(40, 10)) pts = np.asarray([[0.1, 0, 1], [0.3, 0, 1], [0.5, 0, 1]]) probe = ot.MeshSeries.extract_probe(results, pts) labels = ot.plot.utils.justified_labels(pts) axs[0].plot([], [], "--k", label="analytical\nsolution") ot.plot.line(probe, "time", temp, ax=axs[0], marker="o", labels=labels) ot.plot.line(probe, "time", temp.anasol, ax=axs[0], ls="--") ot.plot.line(probe, "time", temp.abs_error, ax=axs[1]) ot.plot.line(probe, "time", temp.rel_error, ax=axs[2]) fig.tight_layout() .. image-sg:: /auto_examples/howto_plot/images/sphx_glr_plot_analytical_solution_002.png :alt: plot analytical solution :srcset: /auto_examples/howto_plot/images/sphx_glr_plot_analytical_solution_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 97-102 Comparing data of a 2D slice in a contourplot --------------------------------------------- Again, we prepare the results, by extract the 2D data, we want to inspect. This can be done by slicing the original 3D results and selecting the timestep of interest. .. GENERATED FROM PYTHON SOURCE LINES 104-112 .. code-block:: Python fig, axs = plt.subplots(1, 3, figsize=[40, 10], sharey=True) results_y_slice = results.transform(lambda mesh: mesh.slice("y")) y_slice = results_y_slice[results.closest_timestep(20)] y_slice.plot_contourf(temp, fig=fig, ax=axs[0]) y_slice.plot_contourf(temp.abs_error, fig=fig, ax=axs[1]) y_slice.plot_contourf(temp.rel_error, fig=fig, ax=axs[2]) fig.tight_layout() .. image-sg:: /auto_examples/howto_plot/images/sphx_glr_plot_analytical_solution_003.png :alt: plot analytical solution :srcset: /auto_examples/howto_plot/images/sphx_glr_plot_analytical_solution_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 113-117 Comparing the transient data of a line -------------------------------------- Here, we extract a line and plot a timeslice, to evaluate the data spatially and temporally. .. GENERATED FROM PYTHON SOURCE LINES 119-129 .. code-block:: Python fig, axs = plt.subplots(1, 3, figsize=[40, 10], sharey=True) x_edge = results.extract( (results[0].points[:, 1] == 0) & (results[0].points[:, 2] == 0) ) x_edge.plot_time_slice("x", "time", temp, fig=fig, ax=axs[0]) x_edge.plot_time_slice("x", "time", temp.abs_error, fig=fig, ax=axs[1]) x_edge.plot_time_slice("x", "time", temp.rel_error, fig=fig, ax=axs[2]) fig.tight_layout() .. image-sg:: /auto_examples/howto_plot/images/sphx_glr_plot_analytical_solution_004.png :alt: plot analytical solution :srcset: /auto_examples/howto_plot/images/sphx_glr_plot_analytical_solution_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 130-134 We can also increase the resolution, by using a large number of probing points and by resampling with more timesteps. Using a logarithmic scaling for the time is beneficial here, as most of the changes in the results happen in the very beginning. .. GENERATED FROM PYTHON SOURCE LINES 136-145 .. code-block:: Python fig, axs = plt.subplots(1, 3, figsize=[40, 10], sharey=True) line_pts = np.linspace([0, 0, 0], [1, 0, 0], num=100) ms_line = ot.MeshSeries.extract_probe(results, line_pts) ms_line = ot.MeshSeries.resample(ms_line, np.geomspace(1, 3000, num=100)) for i, var in enumerate([temp, temp.abs_error, temp.rel_error]): ms_line.plot_time_slice( "time", "x", var, fig=fig, ax=axs[i], time_logscale=True ) fig.tight_layout() .. image-sg:: /auto_examples/howto_plot/images/sphx_glr_plot_analytical_solution_005.png :alt: plot analytical solution :srcset: /auto_examples/howto_plot/images/sphx_glr_plot_analytical_solution_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 5.017 seconds) .. _sphx_glr_download_auto_examples_howto_plot_plot_analytical_solution.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_analytical_solution.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_analytical_solution.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_analytical_solution.zip `