diff --git a/docs/conf.py b/docs/conf.py index 11f1ba0a858d6d1807d47c25e7d1aa5f16320a62..152e7beccb9f6683dd21449480de1f415faf4859 100644 --- a/docs/conf.py +++ b/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"] diff --git a/docs/examples/howto_studies/plot_convergence_study.py b/docs/examples/howto_studies/plot_convergence_study.py deleted file mode 100644 index 13abc890513c68fb25dba2e6fa9de76e5155dd1f..0000000000000000000000000000000000000000 --- a/docs/examples/howto_studies/plot_convergence_study.py +++ /dev/null @@ -1,115 +0,0 @@ -""" -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) diff --git a/docs/examples/howto_studies/plot_convergence_study_nuclear_decay.py b/docs/examples/howto_studies/plot_convergence_study_nuclear_decay.py new file mode 100644 index 0000000000000000000000000000000000000000..c0df5d524ef95f44d04e311c8124f5414407527e --- /dev/null +++ b/docs/examples/howto_studies/plot_convergence_study_nuclear_decay.py @@ -0,0 +1,166 @@ +""" +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 diff --git a/docs/examples/howto_studies/plot_convergence_study_steady_state_diffusion.py b/docs/examples/howto_studies/plot_convergence_study_steady_state_diffusion.py new file mode 100644 index 0000000000000000000000000000000000000000..0882ad1c0813654f58ace2c99658b4e0961e8489 --- /dev/null +++ b/docs/examples/howto_studies/plot_convergence_study_steady_state_diffusion.py @@ -0,0 +1,135 @@ +""" +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 diff --git a/ogstools/studies/convergence/__init__.py b/ogstools/studies/convergence/__init__.py index ea96ca64ab771fe53b1b34b58697acfbce62e6dc..df2e65a5f8223d8e26f3325ee0169f24b0a45665 100644 --- a/ogstools/studies/convergence/__init__.py +++ b/ogstools/studies/convergence/__init__.py @@ -1,7 +1,9 @@ # 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", ] diff --git a/ogstools/studies/convergence/convergence.py b/ogstools/studies/convergence/convergence.py index ed15122f708ffd977ecf37d9108a5d758e0a2109..bcf963cac508e0ea718402eeb7421d9f8e266a34 100644 --- a/ogstools/studies/convergence/convergence.py +++ b/ogstools/studies/convergence/convergence.py @@ -1,5 +1,4 @@ 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"): diff --git a/ogstools/studies/convergence/convergence_study_template.md b/ogstools/studies/convergence/convergence_study_template.md deleted file mode 100644 index b30f5b47c392e584adb26338276956c9b2ad7bf4..0000000000000000000000000000000000000000 --- a/ogstools/studies/convergence/convergence_study_template.md +++ /dev/null @@ -1,101 +0,0 @@ ---- -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) -``` diff --git a/ogstools/studies/convergence/examples/__init__.py b/ogstools/studies/convergence/examples/__init__.py index 11a2d7de1365371004f0fcaba4ee054f8420c447..cce3be86779ad5f3ff7dde84efdc85317bc42eb6 100644 --- a/ogstools/studies/convergence/examples/__init__.py +++ b/ogstools/studies/convergence/examples/__init__.py @@ -1,37 +1,17 @@ -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 diff --git a/ogstools/studies/convergence/examples/decay_boundary_conditions.py b/ogstools/studies/convergence/examples/decay_boundary_conditions.py new file mode 100644 index 0000000000000000000000000000000000000000..03603dc6dffe6ec078e35d4969228b38417d9ed1 --- /dev/null +++ b/ogstools/studies/convergence/examples/decay_boundary_conditions.py @@ -0,0 +1,27 @@ +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() diff --git a/ogstools/studies/convergence/examples/nuclear_decay.prj b/ogstools/studies/convergence/examples/nuclear_decay.prj new file mode 100644 index 0000000000000000000000000000000000000000..d655b561e77df0758d6323cec396f8927453f9de --- /dev/null +++ b/ogstools/studies/convergence/examples/nuclear_decay.prj @@ -0,0 +1,139 @@ +<?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> diff --git a/ogstools/studies/convergence/examples/readme.md b/ogstools/studies/convergence/examples/readme.md deleted file mode 100644 index c61379f35cd4b7f24e12a3fecdda0d60615dc9ac..0000000000000000000000000000000000000000 --- a/ogstools/studies/convergence/examples/readme.md +++ /dev/null @@ -1,36 +0,0 @@ -# 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") -``` diff --git a/ogstools/studies/convergence/examples/square_1e0_neumann_0.vtu b/ogstools/studies/convergence/examples/square_1e0_neumann_0.vtu deleted file mode 100644 index 83210032940bf387e491306b710e1df49ff16516..0000000000000000000000000000000000000000 --- a/ogstools/studies/convergence/examples/square_1e0_neumann_0.vtu +++ /dev/null @@ -1,27 +0,0 @@ -<?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> diff --git a/ogstools/studies/convergence/examples/square_1e0_neumann_1.vtu b/ogstools/studies/convergence/examples/square_1e0_neumann_1.vtu deleted file mode 100644 index c80ca0760500c6f21aea728565d29c9744df97b9..0000000000000000000000000000000000000000 --- a/ogstools/studies/convergence/examples/square_1e0_neumann_1.vtu +++ /dev/null @@ -1,27 +0,0 @@ -<?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> diff --git a/ogstools/studies/convergence/examples/square_1e0_neumann_2.vtu b/ogstools/studies/convergence/examples/square_1e0_neumann_2.vtu deleted file mode 100644 index 35730e027e11d9a43c32b25a6a193da5c5723c58..0000000000000000000000000000000000000000 --- a/ogstools/studies/convergence/examples/square_1e0_neumann_2.vtu +++ /dev/null @@ -1,27 +0,0 @@ -<?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> diff --git a/ogstools/studies/convergence/examples/square_1e0_neumann_3.vtu b/ogstools/studies/convergence/examples/square_1e0_neumann_3.vtu deleted file mode 100644 index 08b85f5a69c134744f86a1d00051307d6f3b786b..0000000000000000000000000000000000000000 --- a/ogstools/studies/convergence/examples/square_1e0_neumann_3.vtu +++ /dev/null @@ -1,27 +0,0 @@ -<?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> diff --git a/ogstools/studies/convergence/examples/square_1e0_neumann_4.vtu b/ogstools/studies/convergence/examples/square_1e0_neumann_4.vtu deleted file mode 100644 index 9f671bf595828f31bec6eb44ec9a3b9d79a98950..0000000000000000000000000000000000000000 --- a/ogstools/studies/convergence/examples/square_1e0_neumann_4.vtu +++ /dev/null @@ -1,27 +0,0 @@ -<?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> diff --git a/ogstools/studies/convergence/examples/square_1e0_neumann_5.vtu b/ogstools/studies/convergence/examples/square_1e0_neumann_5.vtu deleted file mode 100644 index 0e33fa7fc29a8f5b67e70a274d7b8d9e4abc0df8..0000000000000000000000000000000000000000 --- a/ogstools/studies/convergence/examples/square_1e0_neumann_5.vtu +++ /dev/null @@ -1,27 +0,0 @@ -<?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="1089" NumberOfCells="1024" > - <PointData> - <DataArray type="Float64" Name="pressure" format="appended" RangeMin="1" RangeMax="1.675354863" offset="84" /> - <DataArray type="Float64" Name="v" NumberOfComponents="2" format="appended" RangeMin="1.9495981051e-15" RangeMax="2.9076480651" offset="10376" /> - </PointData> - <CellData> - </CellData> - <Points> - <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0" RangeMax="1.4142135624" offset="32156" /> - </Points> - <Cells> - <DataArray type="Int64" Name="connectivity" format="appended" RangeMin="" RangeMax="" offset="34828" /> - <DataArray type="Int64" Name="offsets" format="appended" RangeMin="" RangeMax="" offset="41584" /> - <DataArray type="UInt8" Name="types" format="appended" RangeMin="" RangeMax="" offset="43492" /> - </Cells> - </Piece> - </UnstructuredGrid> - <AppendedData encoding="base64"> - _AQAAAAAAAAAAgAAAAAAAABQAAAAAAAAAHAAAAAAAAAA=eF4z0zPRM9E1NDQ31E1PTjEyNjdPSQMALjgE6g==AQAAAAAAAAAAgAAAAAAAAAgiAAAAAAAABB4AAAAAAAA=eF7teOc71v8btxFlJw2bjERZGVl1GqUkWyQJGSlKSLJCKKFECUkys4lUkoyszCgVUZRKRT7jmij3+3c/+Lr/iLsn19Gho+M6X+/XxMLyvz8YsPz/z//7+TVKZp8DOwZPJSX9qjkx6MhQt+fnxoDL+T4ZzIeBcn/XoZ9rMUhi/VPkux4DzsVNuizCGBTZ7M0qE8Og60Z82hkpDMyEwo2tZTEQ5ufZZq+AwbSQjlnUNgxu3wpXf6+Kgcg7f43jmhh419pcUdDFIOj0nwTV3RgszihwxJlgIJtbnqdmhoFpkcMpEysMPKa4bvcfwmCXk0Hu9FEMmF/nPhV4YuCz8IapchoDulZHVf15DO5+ub4+LAaD550hJz5ewyBu2mnNqbsYnCvLnhCowuC3fH+Abzv6fgZDZ7u+YNB/jeeHhCAO9I9Dg8rn8P9w4Bscry1F9ysJiog6obs7/xYzhNC9Z7X4VSdFMdhkaN/7cjMGrhkWoh1bMWjw+/PsF7rrfbT1I/2dGAwIWaY0o3uuJfefu2j6P/y228ZYYjDrXs4x7IBwgOPp4a4YaIoonEv0wWBi7dZR3kAMouQn1YQiMBjKwmtqLmPAE/FGfCkV/dxe5PXaexi8+LBOZW05+n6fU0pkG9D/c3B/VVQ3Bgn3L6/2GcXAsE8iSXsWA35vmREvVhxGvzGVC0RwGJFpEPujhYNu2ZpHbxxx2BpByf8Rg4OWzPJM91Mcdt5oymPlJCCpYMlPNYj4D4dPPKv96Vzo+zyPP0IXwuBK0177NZIYtIgdEtFAd1u+NSqN3oHBDlEn1wUDDL4dauUv24eBVYq/zB1bdM9LjbYRF8SP/Sf/+ZzEQG6n2Sr3YAz+aYa69KF3wm3Mdzdex0A8wvEHoHcSvz4g7V+Gwbb9Hh/c0H3i+RZ/1F5hoGvhbcIyhoHASK0M7TcGq7mblY2XMdj85tNqqfU47HlvI/lOCYe17H9lO0xwaFti6Hi74pCyX+H40Ys4FDgJKF3Jw+Fz8NEzMj04uLvPbVhk4lCzJYmjQ4MA0eif93LCCbh4Zmt62TsCohTUjjaYkf/hoHvhZb0B4kGqGfNOHHp/s8grC18Qr6PWrHL00EK8tU1OEES8VV8lPcG0xoCV96qBInrfkLO8NpWInxu/ej6+hd71wOAPtfkkxNO7OvP92RiQ8XmjJhUYdC9Lk15NGFTkpA2ZD2JQyjTQXIv4WdrxrG2IikFODC/bI24cmAztm5PS6N341r5P0MEhLk9JasYGh5JZ/2LT0zjYyruV/E3EYYojTiSyHIePrn0ddoM4DLgnL+jR0c/X57oqyBEwpEjjdHck4APb7YGCNAISGYak8QcCZgSqg+W2kVCnSe/wuE5CvnpPWPgqyn84XFKXfHce6V/+YkjQojTS9+ZNLdnqGHw5nvLdyxgDGeOrIT52iDcftW7UIn2udxtcZY90uaEm+YZLAtLHx0e8U3cQbx9O1+CVyBcUN7sWtiIdLXMtC7zDQKIidNYCva/Cds2oMMTf9OZNdfXCOHRu61iQVUd3W6odmD+Aw7u7Xr77vXEIuHrvsU8sDtkRxta38nHAnArEJdpxMBSSKveeweGwhETAbwECjPL9O030CTBEslH1JeD9q+rl3PsEPClKnHg1jv5+eyZSZTMJCw1q6bn+JOQd3pcg20MCh9K2nIodFNjXm3S1omwFB5e76/UKkB4Gq0tqHRAPjIXi+5X0MdBIX9ZVQ36Vv2O32AUPdF/ZbY31FzBgWzJfvwn5kaSFtN2tfMTfu2Fd959ioBdgNGSL3ntGib+i7TsG3wW4hrkRr9935QfuR3fv1bvOdncHDluO/vmtaIl4/rDphJgvDhvY1JNzr+IQX38qd7YU8X3wTJtRHw7qHKz8PzEcdFKfGOkJE1BI5dh+z5iA2dgt1LSzBBSoRZ55mE+A9k9t97pRAjZKhp4I20hCXGDo8IYjJEQs3lt+WUQC7WXB1MwCCec2XEskjlAgq6fznGIXBWZGWjcHG1JXfJI+4du8EQPr5ASfVGUMxsIyqxP2oDvMjSyeOiO/e6fapxGEgdvDbE8BxHtumSFVtwIM7mwtczd4jvC6U321cQSDy5vJCdo88m0r+TQuHhy4E//tElTAYffEgIPKHhzK6pktcceRbmNtowwu4SBZ9OrKpQLk268KO6924kCVve2c/BuHTyGmkv3rCKiL+xGaZUDAd1/VGvOTBGQ+ykq1yCLg12xJlVg/AX1xbb/3c5Bw3PG3wiVjEkDWXDs1noQrYsslLoMknNcyJWakKfDLBX9uHE6B4SEnnaDP6LNGwPKoJRVWZUW1Z3Sv4GD6zFH0lQjitbztfLEGBnYPpzyrLDCwGcx8SjmB/PLhqqqUS8hHOlVb0nOQLxhecZRE/vZE/yqHMbrfUmG0S5RAePbct28XwOEeIyrwqAoO9R7Nx9jRux9q59w7fgaHzWcv9vGl4jDdsLq57REO+ZZrXyuN4XD7pdvVcFYCQniWPHElAtKuHBMbciDgIyN6+8l4AvScm4uFnhBQvXdxi84sAc73Hf5ayJMwnjNv8uk4CQacrjHx6N2faw9JU+ZIOIqtL5w1oID363np3TcR/yG6rgGnQDB7eYSMExW+Z5YVi/RQIemVjir3Xtp/OEhdvyT+QhwDyL5EpKFcj8qy80lHeTfB9+vFOOLBnuHM1miUZ457+hlZ1RhEXpyfNBpAPmFSueXmHOoDg5t/lfLhMLwrlrsQ3f9jk2FEnjUORqYd5e1BiO/mEfM7MnGA+/M8G17g0MRyrD/jGw6KqebPZvkJsA+q/nJAD/F628LA/Amkd17aDelMAs5OtlGVegn4dP3NsDcrCd2Xdx7bqU/CU/59gbyhJIR6mMZ6NJKwke/rTBcbBbjn1mNDVhTYpDrDIldAgRu/D4t7/aXAhsQAPR1XKlSb87SovUI8uHbIokOfBr0cCuomT1Zw+PmT89xd1IOmeaaCbADl/LbepztRHtTtfeLoE4XBI9fj+FIuwsNu6zcB5H8Z3tlXniO/1zv710WKA/mX9TtXh604TPKuNow5iINwy+SL1gAcyvdhnAboftVih42bWnAEQ5525E/kC+oS8XfXEyCg73rtgRF653En+gjSO1epgJkX0nvjEuezLJRrM3KsQs18JPw5VXxdcT8J7N9fFGtfJsGTc1uFRTcJJmIp1M/8FOjIKj3ig3R/NMEz52Q5BVrGEw77slFhr8laKSq63yDW6smVdio8gcrBW2o0qH7QInq+kAYRryINvKTp/+Fw1b7N8rgM6k0d1wc+IV/QO508kYL0MJH2TfNGIgYLW/cIL6IcMI9ruDI7hP59Kn4uhoZy9tD9pfdiKLdsFYsEjHEw5k0rMDuJQ+TwpEIp4n+EulCtXSMOFxRHzgZ8x0H7NFe5shABJ2VIzSfo/tocW9wwkICtHNMubEUEECqaijpjBMjf8OyXW0eCotFHBUELEnoDU1gskkgIydfzU+0nwTye6/dGIQqICLKKnnShgI2cYltqBQX8MjAuVlYqHP9BT93qTIUPN28eC2mkgsod3/v2MjS49eBnvcgNGliq7ODj5aCDSLR2t2PMCg5lm06s5pdH7z2k+XTnAeSTrh9+ZPpjkOVxsd0iHYNWP2n2oEYMMk/kCWxBPOAMq4J4LhxKizdz1aDcC48uMeg6gsM/CNo7F4eDiGvFGdMalOf5J1WZ40gHk0WRKjxI9wdvbV2Hci5SUUPvjR8BjtqOroUo51gkXOyfo/e/vV/czkGABOuFMxG1B0iIYg/U3HiVhJlwB88+lHdBTdop8msp8D5RMi/EiQINW4ydJB5QQDVQcHswkwLtBd5hn62p8PuO/arnlUgHWy/d9F5LA/PHoYs7QmnwpuZdg88MDXTz88p2H6NDlymf2/uxFRwmKkH+HspLObkgQyeUkwkN7wIsQzDI3tKhko30QJqWnHZBfVZpLu1NHY567FfB6SbEA5vLLRp3TXHwydh672wgDi8OnIo6lotDC+9GIrEfh7qqKEGJvzgU/cw23K1CwOHHssGi7gR4supcnkwngPKVLbKtj4CK1vZpKvL9DNx44yMjEprrUy/oRZMw/Ol+VU8rCSL93UVlnBT4JrYULIf07+pewnI5mwIBz1N1pGcRH9yClvYaUcF4oi3gQDYVbNSuvLu4QIW8X9QhVVcays2B4/mvaFBs31W4Vo8Ox4simopr6GDek2YRrcz4D4cLf7c7CyphcNGzdcDQHuHCL51dcRHtnNC1NkGlGIhGsC3WvUV9iHZ733UWHPCLb+pllHG488D1ZRrigTsmIcGK8l/p97HhVNTfHz8eSDyFfEDtqkBzrRgBNi4PyZtWBMyrv6fqIv9vmrITfvMc5WKNlf8NGgEyFimri3eQsOue2gbzABJkeoTXP6slYUrFY1iNTkIbC/MItpsC8ayWjwwSKbDzV4TX0VEKzL63rcvYTgU2A7WEnXFUEBfm+hc7SYUz9hsT24xoUPJd9p1GCQ1GygeOC66ng+QLUiIpng6nu/XkO5boUP77HT4SuoIDW+hCvPd21PcMmJKzhzEYfWfIXhOPfDLmW/1QLcqHFppSwCQGp96m/Knix+ENF79FxS4cKHL5EUmo7/7Mcev2v4fD6/ax55GvUZ932+LwkZ2A8g2cqjU6BHCI3quR9ydgrqbF6MQD5I/TLufypwiIWVhQWS2Jcn/IaLDrKAmST2qNRO6R4P7TkkvjCwmnXBrVrBQp8JydTaIsiAKMClLnUisFBLdk+G8SokKQb7xfmQ8Vbl2Vyg5vo8InloapV9I00NDxEx+KpcGRIGnP7t80OF3mrjN5mI76glWFYS8dDoS/blptwoDZlvUs+1tWcKi7Jl92VwX1yLudBk5oL9XRbrYdSsZgX+jHxFrkCzOstgdvoF781/Tk/rXiOHDm2vjbW6Dea9sQFxWFQ0KI6d76Whw0npYni/7AYdxLMGJEnICXU473+exRv7UI2jSbTMC91s667C6ExxpSzGwVCfzf/03LmJAgay1U6BRHwhZOoEihHEgjRn3LUA6EFZtEqB6mgKIvXetHIQUqXof/4qZS4M6G/qK3+6mgeGY09OJ9KnApDWQaLVFBQjxx96GjNLB2uHVmroUGf65sSdNVokMvEbDGI5MOOeGSrUk8DLjN2rTUF8sAodTdVtqszP9w8IzTxmvVMHiWoX/xsjsGsfVv3bPTMNj+ov/WOrTbR/Zon+ZD+XBa62NTJcrHb5EstA1HcVhXYml08gYO/knDD3o7cDClDSofX8JhVS4raatJwLn2KuuO0wQUeTdTp0oIyGouWR7+hvrhZBl/uxwJHf5edl+8EP/f7b15tJSE672t0h7zJJiNNb9U0KGA4R7+6xOxFJA6HZDxdJgC52V6787JUcHknvXFh2FUmDzzI0vnLRWsuZ6zPFanwS9H3sQjN2mQaTDVbcekQV3/iz89HnQQvJ1b+nWYDnsviAo0mjIgROxqiXczA9gzautIgxUcRo53F5egXV3x72VkoBcG5awdry9mon5pNiEz34v4oLts9g3tBP+59MRYTRzOj7LUzaJ8lF03qL3nPuLFl2izqveoR5uGMazWEuD/qsLE6gDqAyq2mm+QHyinjmLsLwm4yfK9j5uNhNU883brEA9uijUeNEV9INctqe5NL9LFS/+SGZSHD3ZPWT46RoGbm2OyvVEe7jn/ilcN9SG+hBpzIxsqZD0df9hTQoWoXplHi+w0eJDzTJ3NgwZq3hNlXJ004Bo19NRTpkNt/Meop4gHO323byxYw4CdsZ4uQpEM2FDbNytGY8D7BDWl7qAVHBxrnbUiNTGIN3soJH4Sg3Mvdx/9h/pzQ8VTFtM3qF/OHNLfgHLyek9m8XVAPllRZT9yHgcLpup33mocTkTScp3QDlQWL46alCFgTY+HZdcxAo5nCDzUvEuA4i8FLeuPBFz4s9fJWJyE0oTvKZquJDx+hOvZFJLwl0ft3svfJETa5XU+0aLA1ZwLnkcuUcBFyz/9wxAFhJnaOa6IB9KcioniiAdZ01i7/hsqtDhQE2dUabBYV2h/EPWC4slEpesUGvgnbPwzdJQOawVZcna+ooNRtXDmmC4Djj9Kpb+uYoCA/y660lYmXNEujyOKVnCQKOaTVdPGoMgozeqSHwbHkv4VCKMdZcmf7rNhDPWF7EXZbEHUk2uPVPeZIT88O2LWj3ZS8oON5S3PcbjVkLO+F+3/d7qu7oJID9cfqQs1BKD9e/aYy0wNAb8PJ2zrxAlYb4zLhmmSoPCKUbAd9eIG/0v1nC0kfLHuXFTnosCA45/0D/YU8Pjxmvq/XvyRRbbUiEKBnENlQk7ID+o2Ooqk51JBTK37iSTKReq5rBSJwzRId2Aee9JAg3v9jmcFpOjwuCuA1TSBDjF0r5IQKh0edNU0PfdmgEtMzfDWcQbwZF8w6DnEBPY62Welb1ZwGFA7wN+6E4OwttiAQtSflBT9bvQ9QH3RS5PqjHKi37PlWLAI6gNz1kP6djhUvjlUPHwNh8QFsWyXVziQUosTLJwEPN91S2TKhIDtuXciZWMJCM290j6H9MD6J1zLm5OE2Y2BeXWoHx1M6Dzy9wbKg75Fed9REqI7rsrIy1GgSiS0zC6AAn9LrytIoVywZ0yPV6NcSHfmLlQ+SQUar2hOcysVLnx9PZQkQQO/lw59jyNQPiY6ynp8pkGsDue/mr108DpTfbOpmg5RjqHl+RIMkJGs33w+hQHcY7cvmHEwwXFiS5hiFBNq055iEn9XcLgelFymivbVR5+8rzGBGLieyt4vUoFBbaDpLeHvGLxL5qjMl8aB5YJp7ZgzDoKBzDXTGWgnsjR6f3mLg/xd649MtJO7U7g+WNoS8O1IrCbnTQKMhw3U9o4Q6PtanFUXIWG3cCwDO0aCKiuPVFkxCZ09vaGxGOJFc3B7AdqHggVuL9STKMDma79sN06BZ/NRbWpqVEiRE38wdxnp4NBN72rUDxiCEdFZu2iAK1hO9OfQIDXE56QfKx3apV+UpfnQod/UyN75DR1K/vSNTRgxQDZoUEDlEQPUJLQ6nZWY4Om2FB9YyARW9ZO2oTILK3xo3tUUroeBT/RjDdFgDLbw1m/9h3blvSyP0AMoL9k4RLbIKOBgIj9xtMQDh/az7y+x5KE91WlrYPkZh8NDtMB6SQK03HudjrgS8NAx+Yx7HgED0BkyP412RLVclKISCZ/vWnRpop4ULZv1XPsZCcELc5/tOdAedOEuf2RLgankxNPh+RS4FrZeZgDp4TPrTv0hM7QTgjpOlORRwX1Nt+V5lI+61VEmR4/QIKE4mAhtpMEdFWGuOaQHOWHq6OsrdFAY3M6ym6SDp5EKx4HjDCjslOBdPcKAQC2fy7HmTPDbM0151cGEZSzm6aTJCg46zi2iRfrIJ03N+E6jPs1n9edpWB0Gzb5vBPF5DOi3ixTI7SgP58OVC0/h4PEkNVqhFIeXXLGV95E/ZktWBWkoEpDLb9vJ7ov0sH2hR6sK7Wa5xrdzBAFn6pucrXVJEBVJwmJiSIj3SVR6gPLhraRk6vQmCvyB9xePeVNg7q3xP50nFGitwHIzuKlg1RlqluNGBWWV2L++DVRgarHkyG2gQanqwfzvQcgXRja8HxyhQUduNOcS6suF0eaVcQV0GHupWhYtwAD++EWMPYoBDbGSJlIEA2K0/bDPPky4Sa3Lt/rKBEOZSO1z7v8PDl41BeUGGKTMrdcMC0W76tf2ouuP0a7q2sPkoaB9Ieq5b8MOHOIej7X0n8VhV1N2/WG0ozyItRmjGNqZLdzMoB0EqKn73jAPJuBIlaxqXAMBAYc3NRktExBPnlEuMSWhKH3G/2MKCZkPK12Wx0iQ4oa8nVtRX9okmPUghAIPOS5rBr5Cu0mXv7dBggq3Lx2byz1HhbDycbW9A1QgunOp3Uo0oD4sUTBLRHhk7l2amqWBPv/3ikpbOrwQoqq1N9Kh4eRtht5WBowHOpfJZTJAR1n2+w0eJtyPcxJIvsQEh/0l0lLIF8TFvydaha/gsKf7uHHKLgzWDbNGGYWjHpV6ul63AYPv+jfn79ARHzYpmyZo49AndHBaDuXlaGAp163HOIhz+/EwUU58fbrGLkQP7cjpqSM7LhLQLmtaYoH80TXGt22WiwSN8+xbNGxJSH65SmvPXRL2ZRyUMZkhoQtv+WCrjXrCzoqJ9MsUaLOUVFFBu+HnLf1Nu1SoEL7hMtfXeCq4TTudAOQLbLNcVmeRL/yM3+uehnxhKLoqroWFDtZOXQ+ETtBhGzXnUv5rOuiTPU4xuxggtiG2vbWSAZLbEkMCNzOhiXi9MTeTCa1T3gLWQguwK0o5PDF1BYcb7CUzVrsxeCTy4E1bBAarxpc6w1GfTpupGMxZwIBWELXqoB7alwdGr5SG4ZB32yFhtBHtDPVxWEZ78tPsS2yXEerLI68nn6Le9HA61eZmLwHZHWPB2DoStPoWPUhntKOYXquakT/W8PsvBZNoP5VTbukZUyBNvfKuyk3UF0O4TYO+U8CprJJF0YAKA67XLP1vob7cujHdcx71Bq6WX/IWNGhc9WhqoBLx4X6AdaggHdQUT4vsukAHvtd/CNUvdPj4PkXNy4oBdvTXKXOoN0Y/tOB6r8EE61jOku3lTPjkJOFDk1+AkZHtkrqFKzhIfHl3bQ7hYH7ZReZfJAbtTtc67ZvQ3jhWdGX9XwyEFRIvuKFdNbasV+t8Ee1sW147xRZ0/yGHiFk2AvDtppbdpgQQmlFTn5MI6H2mt9thmICotbetLUVJCElvG/3lQUKJ0ZCLZTUJjabjfNlLJFxddKkizCnQI5oeej6HAtjY77+mOAX6BgIYsfuokPAhY9tBtB++HY3Z0bRIhd1xXjakEw3srYpMuJA/StR7rRFG/rgcWFSlfRn1xqmfm8MwOiT9bGBQXBigrapq+mQA+cS3Gf8BYyaY9eRy7X3GBC23i09ktBbgMV1D2K9uBYfg95uVfQD15/SZ4yTa220p6uXtLzB4/m1JcPkf0svNxqIO1CN5k2rP68fg4PRFNj32JQ5RTX/4WlFvcB5z+LbRHOWkdZFAyQ2Ew9yhExnvCbg6fOX+amkS6Et+umKnSOBbDj9BqSfhq5s6az07yslv+p4X7VCP3s9yOKCIAkciz89WMikQ8Tb8kJ016g3zw21RpVTYwbVNy5qDBgonyq0nUX9e1ipUPtCBejRVcWuxIsrL2Zhq7lTkC/HHNyYv0iHlpJSXwUkGuFVOq2wbY4CHkJyvlyUT2IwP2FNRTnzb+uTgF8MFeLu5v0eteQUH3Uc67c0Ih8U1Dan5Uag3de8K7G3G4BAhddaFBfHg7Of0cCMcjsHL2f2xyBcbDMNm0K6yUBzXj+Qi4FKI1mYlSwLC/Kjs624RoNBgm+OIevSw/K8LEmhPtZwsn44/Q8LEvcKJWpSXO87FMDvWUGBVct7BSbQnizc3HpQppwD/nvU15csIh4Jgs2wHKqzv3zDFX0OFDUpj3LK8NOg18l0gTtHgmcuzQ5l9NDi07THnNnU6uF9ew9eKdoSM2qD1GXYGrPpqE7IngAHWvC3PbL8w4A/Lv44iByZMNoen7xtgQhw7zQXMFqDy70RoatcKDhwVO1F6o/4Y1VPJE40Bfz6aGS0Y8EhKF6uy4rDvRrvdJWMcPud3K8bF4XAkZ2zL4S4czGl/d2/mIWCBR/EQmzXa1/rdOyRvE/Bl+HNExgQB5ZuSim9sIaHXI2dyG+oNKqDBmdREwgERvsuDPKg3jGbvFz1KAUezhrSkKgpk+It9M2ZHOlhd9uHYESo8t2HQl+qo4HHO+LbuWuQHYn6Fqv40sLG+G8d4TQOlS6oT5Vp04E4e9nXNoUOle9XoFrSn/lrcSRU6zwDdFOsszR+oP4XknLntzASBHUE/jVGPXhyukzWyXIDOPrcT6X0rOPAqTb/jMMTAmin+RAXhIMXxoTQH4WBZU/wkGOFwfdzvxyjCYTzVy+B/v3/0EMRqhxEORskf9CsRDqcEh+pzEQ6Nl4am+xAOrQPCJxw/ERBXfP2BuwIJ9d3m/XyBJHQrCpdEviBBnO9Z1jAvBZY+Locru1BAyTk6tqqaAnajm6WCV1Ghce7G21xnKsiy/RI5WE+F7Veo228Lot4UdG4m8ywNdO2FM08P00DE9lPvtp10cIvK2fbzHh24XlWYN3Ax4LxNWmxlCMJjT13X2xkGdCnTL+u7MMH0qcm132+ZgGXn/P5ttQAlsVnfDQcW4P8Ae6DTpg==AQAAAAAAAAAAgAAAAAAAABBEAAAAAAAArD8AAAAAAAA=eF5Me3c81f/jPYm2kCgRaSiVVFJGEoUoJCsktBVFdiVC2SQjI5sQ2dmOvfee1x32TrLL736+vf/4+ec+zsN9vDy8znmecb3Q0FC/7JD/vxcaLpL4/70uxeR3Mq3vEkkLwv/hzRnibH02VikFCf/wnmzxOJdT9+I4M/7hF1H5Y8w/JEcmcv/hIgdx25P1lWpjRf+wVYr4BTfa6ZinZf/wWoo4J8MquwhN1T9MdBP/HX3zamBC7T8sVJTvUeHt5GPc+A8n5uV/jXa5eCCu+b+fV5vvuO2D0+6E1n9YMFpcZW5LSDTa/+HJHPFJ+ZqCzxOd/3BIVb7n9omsvad7/mGmsvySM91xpZ/6/mE733zdHGHLMvaBf3g0S/yHxMo1WhLxH5bLE0+sGjE8ep383/v9xINuxvv1NFH+4ZSv+SmaVu8evBz6h49W559TYmvkExn5hyWT8vn+WuUtHh37h72rxImBf80KFSb+u16auOubG6lPk6b+e39Gvuutvo9Sp2b/u35N/pPOpelneT//4a6mfDtv8W2XPX79wzQJ4laK5umrxb//4fAf4nIF0SLGHkv/3d/i/LMufBzTmmv/vb9UvLEv+wWBmbboH+4WtwDXnz6pTUUX2P4/3s/1iEcM/I/3IIxcEH9ZGxuEBpe74kLEIMyp8Z4RSkvAnqs/N+/nCIZzxZ/YpV0ZeB7WrsnwLBglz88xniLmQn1TkHB1RzACt/IkDJGK8PpelqSHdghqkmuPeBiUAY8l79SvhkBZR9jyzu9KjNNp/HiQ8wVeTyWP9nypxaQnQ2e0fyiWhrd8rbzbiNkrv3cb+IdhhPfe2S9hzeg4R2uzQzccn0vY3FTDW6G8tp20nSECLjZdhvYZ7cg1dHzeVh4B+SEvzSFCJxjCxBK/hUfiyI2KLV94e7AlXeBPQlAUdp3fKC/s0If0Y0aJq9nRIN/dcqFu4wDm796zOb0hFvtFz+5n6yBiyd/qTo3tV9y3fx7kcImM63aO+aMi8UgN5qfjKKFAc+pR7ROpb2hl1y6hNxiC+3H6qo+JSbAhrTKbHxrBNH9EgJZ8Cip8tlT/3jmGx39vdjYyp4HI1DZLc3ICLWdHczkUMsDknHOv+f0UXpZ5p2UfzkKv3hGVO2yz6B+hrf01kwvPO6EDZ4J+Yrn3J7MYRxFUT/rkkXR/4U16eE4GTRneR7lKir35jYTxqqBq7RrYl8iestm/BJ+lQKUvJ1px7uCcl6f7Kmp/7yLWniAj235yxw5BmqIA10+nV0pHcY1z6NGrQLoit3P/4/2/c7+hQ5xJ9R/voxOnfMIcg5BtcsfXrS0BZmepkvJKwG03ufFjjN8gWHKYlWYlHX/CHbVi7n6DYVZRxuGCXBwaKLXtqP4GkY90BfG1ReiSvi00cCMRNrfcM3Vky6A800IvNJWI8OdJbeNtleg5poi4eOr9a9vvV2ZZC4an7Ct2Tt/Rrb23M+BsI66Ze+Y7OSTjsobHqsD7ZlB02/u9b6UAUq5Kwg6t+Os+eE56JQVGMnJMVz+1Y7dc1ivj7FQUqs+cpcuh+oNsc7n0xzQ4240tP1vsxgmF7Tn5zunwmBK+O3yzD8bzztWdURkQ+i7cf6uOAAZ39V8To5lYPJh16EsmET7b8vi49bJAKFE7P7+XjJI0nUvMrDnYuOy0TcOfgkOCrhZdLHlwZY1Z8BEYwomul+fsnhfATFyJ9cbIMGIL9F6d4i2CwxPbSdmKUexsh6FxbzH+sh64Wt4yjuEScv3w31JIdFvJMLNNgVy33fjmjwrcLdyefi1kBp/u6SocvFIDUwHrY/miP0E4ZbBJek8jGnYf03vfO4cy1Y4G/b8tyDPe55ReOI8t618u3rjYBd3tRLHA3wvosOftvuRABMtFR5KBxzIq+FPeNl0fAdu2BlIN1S9o28saTak+pZQ9WJ6t+QeOeeyMhOn0/3wyJf9Uc76T1mI61FOig7ZdD8KD+FSJIaYMzFa/KKLVSEDXtG5joVgGbrx7S6RtTAepgm1N1SoDYGr9zeeQi2WTv3QylVS9m24JIcQWYcTf/u2DY5lgfSL8RJ+jDN84zdlKwzIRFK1ULxlWiZpAI9P04z9g17h4p+FCLRTVc5Z7G34gNuqR26fxBvDEvNwt5ZYFg2RvhlM6zUjQnNx59142OB4wSfPdakVfHccGUa0cZNC6/DHTaceV/bIEl5e5yNie8jzAphMl/OIO5PQ88O4naWmkdAOMddt+8hRg5YxoeQdNH7rWd72tLy5E7NY4i+jnBGz/8M4+U6kIJss6Bw3tiMi1DjwTZFCM+C/P6zz7SOh7qHvid3gJuLZ1EmiUKKjbvO5nyloGvz3LvH97BtEQ+vDaUnU59NJp8+86DiNskX0nY2slzNhXpaA0il+FPrd2yNVA81vipUipcdSAd3eYQj10Xraofb03idsXAq6KOzdhn9WlpOjcaazfUtjh8LMFgSMqHNJxs7CbsV9c4ezA1+dN6QLsc5AVSqyNedyDDa7tcbRJv3Cg04H7phARRysqN2Y++I3C5bhfqTqD2N/7tyb78CKiJtviLEvHUCAbJM37fgm6ktPFqaWTsCB19RWxLoNsMvv+qfV/uX+tL3/O+9pykGMufAP9dUL+BKJuJ33ZldhceMe/o5zenAD96cKt+YRcaH+81/z0Uzps0to3egrk4c4BWdH6M7mITF3d994/D4SxHwdHXhYhjT5nLZ09H4LtesZ7+krBKfNs8mV6PirY6DetqlaCc3dnr/XTAjQZeDNmkWog0arbO3C5EKon3CVjnBog+D6n9YII8EK8J/QEbzMsVHt1VPmLwJSyltnB1Ioiv9PztHVFOPn5zP1Rpnao9rktrToUY6Obv9u3Q51Im13z4FErQcqGE4sxit1I+7vLd0mqFBd1r/tEf+zF0nLpyPDNMiQMW32umu3H7Ne9Eo8/lEPTQe117XkiRmrL7Q3HKxDya2x9owsJVlfblssdqnDG3eKNyjQZCTkzySO6NXBVMVa+YjiIv8PFQWWudejq9v5csj4E15DtriL7GuH49msgc9oItCIZtT/sbkapHF10teMY5FfZpKxmWlCX3PY99dUEirmSzIaF2sHTcVRQM2gKVTSmHA5buvCBgTGDgzKDkIi8gn1hvTD7Q56z4P8J9s7bLhc5iXimkf5W+OUchiMOUFiZKAgQur0SV/AL7E6V6hxlw7jv3O+nOzoPhfP3o4j8E7DIPm7HU/QbJDEITTlNQUU8t/3MzgXIvDu9v8PiX++zixoSD2JBXs6bIoR6el9hDgrEfI3lmplPETps6nRPRcdDQ3fOTSSzCOVLfn9aZNNB78SUoz1chKPp7qvkuhwcsBVvVT1SjNLWGc1rZ//3PsaoCvNimN6+JXbcoxQy+wyP27UXo4vGijSxUIFCUsGcg0wJXG5Nb7F/WYP0IpoTTDUlKLLqc93M0AB9Bj7heP1S9L56zDba1QSVwVgWvp1lIPvcO3eqgHpuOC8nVrSVwfU1eW9hahsk7rPtlswsRxOdKTEsuwNteVuUKOkVOHNkRGJ7WxdMG99WFVBzSjlMc9KIoRc8yrLutZzViBjOlXyi0o+eDxeUL/jUYIiLx4WcO4DuL1znf4jX4etHy/f2B0h4S35uPHe4AZ43u55yWZJRxC1kv/dkE4LM63q4iBRciD+c5W/WjPWP849f6g6B0P1iv9+vFvBuzTVjmh/G1aKDbTJZbdAtVjLLixzFVpu+yMmmDvTv4WDLfDqOh0NblBk1urE4qdy+XWkSP744WX1/2AeDC5dX9eSnsRa2t7KRmYjjbkLZ14/PQtO5ZP6HFhlXXb2Y6M78hKhIl0LE/iFk/DgxrMs/h/u0iuvh+qNI1STeIK39L188j/IwT+KO7vKklvsvvFc+umd1fQqaV2yeB+EXCnRCWEWjSv/5w80pcbcRDdXb30rBxs+2rs4biExTsezq0lJcVePkHT8ej75U02+lk6UQKL5tPzOZhoZkr77YI2VQoh9LfaiSg8Qu19o1kzJcqvGNVa0FJFIG379uLINm8fgZQ4FSFHE2l+2ULMe3VaLd7vcVEN2gHNBdVo5iqXzPxzPVYGb0Ybp0twIjJ+++kTCqx50XldfmmCoxuBosrWDfBHr1xns0fZVoM99SfP12C3wsFeUpJVUYTeAtUJJsg8UpP6nPVdW4Knfdj+NCB97z1cTkz9ZgUM/AxFuqC95v8wLSLtVhYS6FMV2/B8qLKgeacuvRobb/clRAH27zvOsIftyIL7P7ucKGCbj69KWXR1wTuihMDLZxRPwa/xSzD81wq/8z7z5JgtYdwRZuan7z0EUIxUhS4M4pbi9r1IYUiXcmLTGDWErM7z3B1oFeftXGPs5hVKHhS9PGLqjMWaduih2Bs8rejfvkeqB3UychWnIMuTqrEhJzfbiy8baz3vg4So7dPKq4nYjxE09bj4ROoioxb+PZMRIoOZ7G1irTuCQgP6TPMQiGtnQ3xcUZ0FwqvdMcM4wlia2vz+TOQqtN7z6L8xhmffhoV1V/4g3j52GVlEk0PpNjXM39icbWEokJ6nXqesWT0wk/0bt9sf8pKv7rkyP5u47KXMmoqgB/rOV2+4DPeDN8uT9isAJrnhMq1uFxYKq9q2O6qxJ7mEN1AhzSMLj5u4KnWiWecokG2pZlwzf/TJLbt0rYcb+t8GEHdnrFNJ1mqwLjN4/rFUUl2PhF7IOtfxUOiEZ+iaCrgCyr+EziyWrEXbcU6lavxpz560MS/dVQiLp6iqagDt5lPaIbv9agVqijU4+zCZukMugFPWrRZ/QyxGS4GQzFFaoGPnUo+eV78mBhKzY1+b0pzqmnntv1z/PR7XjG3HqZdVMjFix2vDYN6cR032q6yaUmSD/XqBCK7YZFzMe2m7uoOaTxlbWntBflmdy/zDa04LTerML6Qj++Wz8dbafmkkohf9pnHiJVb5q+stJtmFxe4Jq4ScLuDRdtosPasdn8byOdNxk3lt0m6U92QkRyOfxELwX0l/VC1391wYw2QffS2SFs3iHk2vG3BwQPOgb9gGGsfQlSDr3djzNfb3wPoh+FVsOum1u2ENERTkd/xWYMro5WvC+CSShipTt+fobaN2QM3i4aUsCa9/xsj9YkHl6Xs418PQSeLxz8V5Kn8NrtZYlc8wh4Db3EOHqmQXOn5uCZe+PYW1Qg/bd8BsbD2Y4fD01hcFB8Px11z9ZsXwnu95wGi+jzTuVrs1CTydXgqK7+pwf90fy1GpZ4y9ZquMvr+JSsBODk+9UrMXPV+NbNNb+LLg4P+lWjBo/U4HPs/rzs7WmIELb7GP2iBvW/i+NMDmWDyexVrUNTDV5rrD7+fb0QGtqyxaqytbCyqd4jf7UEHYW5ibwdtfAVEfDqelGO+09/Hb1gWweOix72Bt+qsOF4pRvf5XpcCeQQEVysBWPc41vtPA14cZhRzl+9EX9br+904m5E8ta09ysPmlE8lRYxz92EXzuSTXqpu2qNo/R5cEkTIj9qHz9C246MexftQq2bcSXGpoBlsAOi4ZaDkrItqFYfeNfa0QXbZ1L+ZwRakbThVbh7Tw+kz82Ixpxug5RGwgbZ2T7wrL94e1+tHTJVtN19ewYwpylG+vmlAxfc59SME4hoGRPqimTqQqLl01enWkj4JfQptDC9G126jqr36Sg4J9cc4+fZi4mruZ9ELg5CdhvbxQ/x/RgVKXu5x3YI+/O3288yEOHe/mpWpXQYNJxnQ/IMSPD7sbA/e/MolqUL1AcnyXgoIiLCLjeGq+da6STjB7HDlPdPjRU1X/qfTEUEDuN0qbqUk9cECvUfWrrlj+KuE31Ksu0k2s8qPPi7bQJpYqcWUs5Pobgk++bC6ynsf583yp00hTPH3A/y5k+j6POUr1H7FGb3e77jLKz7p4ddLeKPw25ctq6uw3WnwvyHNwKwN8VpNGaiDmmhdVxb1L5iTjHn7PChekgQrcduuaRC3SAkMtG8Ht8NLPaHvcyC3MXgnyGEemydT7ojbVSAKk71b156DSA7ybxmKShGuO76jP9aA9Z27tq6h5oj78yOBm7MaITdZQ7RxOlKWCTb/r5wpQlM1r5/a7lqIcXj82GwoAlFTw6km+s2oGbfO+Wjcs3QIpzec6u7CT4j2xt1Bpsh/efCdn7PFgS89d+6x7cFny2WbmUotuHnDSOr87dboTsw4NyzrwNWomv+E+fasLfrjpPJQic6YnRUy0+0w2hf3q1TA9Q9+ir2eLlkB3quytgGd/SiVzT241PzTljcH9m43NcPG/WTadUNXVCW50rgmxuA6Vc9S++bPXC3p2/K30oCye8YFwddH0wWSGv5vGTq60Dfo7F+GG01W2q8RAGvf2+EIz0R+4rde37fHURy7dtW4/MkNCTOnBW2G8KGM82B9R5ksJqoHNMOGQabZDSNw7ZBrJ019En+PoKyprfOD9KGcN/9jWdH8ii+3j1++sO7EZC1GcsG/MdQSnfdwsx6DK9tR62Lbo8j546dTuTHCag6y1LrxziUpS0PKpZMYWt38TtP6s78flp33mpgGvJ/w+L0n06gZHSMwSzk3+eAdhLj+eSfTrM5cY1w+rp3vcrbHw4xPAaR9Y0Iuae6JzsgFm0vPX9n0zeBLq/1gjl19zPy1B9WPN2ETBqnVfP0H7h0X1qH8qQJlEz1b+o2+fiRfdlzKLkJ+cvcFrLHivH6ouw5Al0zzN5voa3hLoOMYS6jxcNmBEu4LBvyVsLq2KN3Vl3NEMwodug/U4OB38ElnVotSLvw6EOdSj3c0hRzU6dbcDdGYHT8chOynKz71fxaUdhg86mY0owLf55MKii3wfj1efVwr1ZcKI9M1T7QjuyLLyv4ZNqR16uS7Le1A0efyg11b+rE9x5+Lx3GThzhLS3e0d6FUNt9k8anuhBrfqiekNSDLXRy32ZfdMNlwc6TzacPRYu6r007enBkOXzfTnsCyO69fL/v9YGPfbt0hTiR6ucvZga5CPjy5nVWuTgJOcafUhw2Eql9514R7SUyNsYx9BzbTcK1FsU9SpcpYH/W0qgmTcbvdLXONMlBsEUlrQh9pOD5SMlRMZEhPL47mneXZgi6183FKIeHceiO5bkXH4ehvWXxgA/NCF42tnW5SY/i0hY1WbmSEQh6hu1a2DUOdjopuN8bhb3nIeO6PxPwfl2bt79/FJ0XmIN2rEzh3THf1c6jY3j2Uu551eo0nl6bIflIjcHxCr8+3+n/PgcWTMrfrJMalX6+Ge83ZtndqPCD8In3afcUmqFY4PrVuS4Gx7XT/W4/a0ZarNCOU/opmJyZTmoIasbhIa7KZ0OZ0NRey2DtbMboth2FzBZ5iD6odTzgQAuM+V3Fd7sV4fHNDJF+mxYcyHLZonO/FMI3Z3e4U1rgcH6fj5pcBSQTjfKf6LTiLdeeJ8YXqtEyMZCkN9yK8557f2wQrkOS4rVwh3dteGfXOKQh1Yj2wt0XWc60wzRefb/8qWbYTHtetF1oh3NUQ2t3eQv4xG03mTR3oEk84LX4gzZwc1/3aS7rhACj7dWkHR0gJV8Qr2npgmzL3zdC6ARJ6f5S/Z9u9O1qkgywpvoFeInBcr14R+FzkBLvRZlHWfylH9S+/2A9UJChH2/js2u2yxNgE3Nep6qdgPJjk25VdERIc2/4oaVKBI2ylNLFZSKYplcY+y+SINZWn31uExlmni9CJQ6TkUv8RdE6RYGpsb5F4iYKtGtm//wxHYTSjtUTQSQKOAy6lrtbhjD6TDo87/sg7D3u+tYojkCt3NJ6t+EQ7G5bNNgPjUIuSdxTgHEYAZYl4/K+49hiIOP02G8YrI4tn6fUJyFxUODqm5VhhDAFXyvhm8aRm9xrzCIj6DRlLZ5mmYFd/IzJFep1yW1u6XQfWv7p4TBFfPNNoR3THi0wU37w+/aMLzi3q3BmR7bA/lFBFWU+Grx+91kcqPd5LuSz79E3ycgToWt5v9CCq4WfBpa3ZqJrvaiKRbgVWctc99UMcnHspsNBL8dWqA487laoBMTLdurZklvhWnLU+3xoCdS/2Eu8pZ5n97t7vPcZlOOW2kNWpcY2bL4R+fa1YBVcvi7a1d1tR8tF7kmjrbXg9Pv8ZmVDB+h4v3AMzdZj9xO9QrrcDtArcvdqmDZh3StBXuJ9J65sS/l59lczNguqHDz1pAtBcj7vz7xqBb9JborH/W6EF56/qbqV6hvE164yNj2w8WaeNAjvQJbXlZ+vE3uxxnMqmyTehaFblhtd//RBIj5FepXQDalcY7Wl5wQYPLxWa+jUi/Uzs2viVP59DR+l6gv0w1EiQHOCRATbJ3XFJ50ETPbTBj5qJ8H55KkLERLUnnjU6PIygYxv73rTyTtIYJvQ0U1dp0D46e4WSgcJBhFeh65cGMLvt7LXPviTYcf+YHOv8zCKc+VWTlyj4O6XFwcpsyMY5xr7WDVOwcTTqEGC8Ri6u+SfMVoO4snGjzVXaCbAIvLgevbYIIrf/Bk0j5gEV2Rmb7/EEN4UtlS03ZjGfbLiFvWXQ3h2cEzr5JEZCMu4fGNyGoLedEcYS+N/fxf60pY/t7maU7mjFXtEA9iVdviib7zk8K7JVoy/PLL3JEs0Kr/Y+b5iaoMOg67XQ77vSNC6N19I7dte1bt93A9noCnK5Oyx9214cnsneVEyB2kMj4K4utog4+j5+/zjQvhfe7htj2g7UszfseQ1FMOGLcnpe2I7zljRcN+h7lT5fSbjPoIduP8rS3SWsxJrandWfSs6kOnfbi7TUg22P5pO0c87sVR/nDpN63BcxDKi9UQXKo1/tVTqNOLGGO2Y5noXRg8NvLuztxklM/SPPo52Q14jR/mXC1X3j2oXacd7IMAk4Ly43oppZsEoDfo+LPrv/0v3qh35j+mUH4r348tm97pnKx2o3Hjx16GPBCQ+2e0k+KYL6l8LtrVuIGLh18vxIdoenJ0gx9I3UHdHl/zrr+97oR35aOLCNxJChSPi/mzuR7bznPbTYDJCNzZfrXYi4C1zhopVGAV6nQHsp5cHwEe4a72eNYhbR9ikFJKp19176Ijs0BAURLsdFDRJmFy3ltrHN4JpV0Uz73kSWr/pMyS+HsVjwfPZTXZkhJjNOZsMjmEjzqvGLZPxVfKkFofOBMRSQ8frtCl4fTJfzY04CaX4CIfr0RTkZhjv/W00DZuKQ7camyjwEnWhvy88A9EQIZGWAQoO5Nv6J23+93fBIhFKviX37ZSzTO3Y53ti4t2BT2C7XbbWeIT6fXrB7QWHoyAuUv4p/Ho7rMJF23hLE6GuFLzUad+OD+zaPB+E0/FcICyioKId5q7WigW7spGUKZnxndrjGC6NrHGs5uP83RNz++07IHFvl63QYhH6tENaAhc7gJuNL6ItSzFEaULp2048k6Tj9Fsvx1W/08+0ObogtflidqhXFV7Y850Qq+pCztKxzWeO1iI0Nv6Clms34l6dLfWsrkfKqfHUev0ezLG70/SJNGHrvoij1Uq90IqR4Dv7iZpjKdy6oWp9GEiwcnk12wLh7p6GPLN+TEbS25642YYAAzGr+O8EVA9Kh9lktGMp7NRHeyr/Kfo/e+X3duI5q0psVj4RDGPXTx2y7cJUK0Nf1XsSPseEvechdeOZqtbF1/pkrMbdaOi/1AsfgTfK9jcoKLoyWPAioA/kuiDiXYVBxJ9IzdMe7IeqiKH6fe0hqBdNrLnzDmCWvkrd5c0wHhx1ak8/QsTrHaXieqkjcH5Z9neKqhe75tt0q79HIbwWni/HT8Lg5Fb23mvjOMB+8yrLJxKap1eWRxMmUNagk5Y7QsKGTw43rNimUOK2+SfNETJ8Uh+81nKZxuYCm60818mwatzMd/bKDHjaJJQXNMlYsRP/I3W241+flOwQL9yiVLhJrAMj45JksZM+qDIPuHlXpQPuLs+39glGguEec94D6w7EfmscUVT5hvnHNlPXUjrwYnC9RuRSGjqSfPkVqfw21Dkp3ez7gUS1FqcqxU78Iiq30LzLQyivr+PBrE7Qb2VC354iPKGPaSMIduGYXXuinkUJJhf4GL3zu/DukzQra08ZKlNfmFC0uvHN8eGRO1cq8UO4zTVkew+e+mhEWGVW46ja9Bf/lh4EXTsR/eJkHfSaYr+XJfXiovQSkTmxAQbv5mYEw/rwfoLy6vnnJoze91Fc/dqPutrqEbv+ZkwrrGarVxOwpF7/9MuhVsR9WlSXovLPU0RzYORZG65tjLr8KImIGyuBpeczqeegmpVU8JCETZzys7vXOmCscM+CLECG+DCXVcKlLrT83LJJn4GCtsAg/HjTDS/5aJWEKQo+V79080jvgUHjeVOGoUFMdm1lZSL0wu3FdhOnySHwesVEpC73wTrJNHaabgQj9q/4HDcQ0NnqdyyFfxRbaS7zKM8T4Ko0fb3p7hiuWdlx59QMQF72uy0lYhzTR1XtsiSJOJk5SL9lZgJfCY8zhFyJCNjXkzR7ZQqjyRoPWIuIWFmLthoPn4bYYdnZFSIR7x8zf2RSnkFK4nzitkkiPnRXsKgp/ntOwK61V9yI6etpJfVOcLuenE08/RF9l4L2yRh3YuehLMdCoQioXD7vdMu/E863CvgKOuJhWP/ueWRDJ3SsxUcJkqkIMnp57zl7F9Y0HC9JhWTCV3YEm0260Eynr9D4NwevRjawNvZ24SuL2Vzck0K48wi/OafZDYLnZXWycTFYHwauFQ9341zpQvZMUilmffPSdjj34E3uZ5sLc+WYuLd1ooza607zhxk1iFUh9691VyzV7+3va7R5u9agfiGr15HShyClWG9eQh08yqe5L3f2o11t0zzhXCNYObco1pEI8D936dPDqiZUD2/eyEzln/XmF47qhWZMbdg4GhFGhNhjWtrTh1vRcqNtaY88CVP7ykstVdpgu9WX5ictGd0i4iECb9rRxJIsHVhJBt2PHNWAqA4IydBXGgVTc16zJzi0tBOnFDNXkt5Qe+Dd/t9dPV0oTE6fMDEeQs123e8O1N+zyoNXL9RoGFG9dppulB7oF+7277IZQd0096PAhl7MOJdn2/iNoqWknZU/ug97bDxeuuaPwWHvZa5vev0Q8UoXS6fuSEN//6+RNASMmv8lXOWfxOsxG6U+OwKcu/VdJIymkNGYlfaZQEDj9w1lNzOnQb/t6Pdc9gFIDtwLV9WaQYmT+cgk/wB4z0gFW93u+tcfmr6Ie08PeYbpd2HHe1W3uVPeqJB15+O06cKRtrIV+zPhcApxiaaP6kLiz90GnCpxsG671/Wurwu3FJtoP4unQLvpjnTNkW5w0lxI+qOaAdErnmd433XjnUnsCamX2dRe1tdHnO4GwwvRZOnP+TA1cDB7ZtQDzbKBrbcNisDicYtH7G8PXF3Ocb17XQLVobLaloheCEzIjGQGlcHYLT2rSrMPztWp9ZT8CvCOmiQXHezHEu+RP9XkKpyqti18SkcA57cD6ZPbasHjTRB9v0TAV5swL/fz9XhD7285TkuE6A7av8H3G6ED3U9Mn4jw1t8x003VB7e127e9giR82b+R4EBoxtozyZ0Pu0mQlzi1/eRcC66xecXf8CLDlb612GpDG2q1zlhwqFBgVKHbWLmjHYkVwZ20vINoq3iXcY6lA5aHhns516n9oOvkjc87OiHn4RJ1d2wI8hilc/jTiTR9oxvexGFYq2n4KA50QfHAl40SgyPQHOEUO5LSTfXf2I05v0Yx4djCuo96n2hyFP23M45jw3eepxd39WJOoJZ97fQEKC8uGeWF9eKQivm2jbqT4NkXakzH1IcfiYL9d72nYDZszzyiT93LnXSfWCunYck9SZf/qQ9sFN0aqXszaJ59k29B1Zuz4en7t1S7/08PPN/mxBcS1l+E6XTDOfT5sQMqnjDK5Hvrat6N+6cqO9TUQ+GdqyJzLKwbDwb3tc2Xx+LWTk/F0O5uXOPsb+98+B27F45O7z3cg4X4KO6GTelQpjyS7LHvgeoBv4JR4g9sOXmN8fhsD65tO/VdF7m4cW1D+IMXvfgsepMUHVWItOU6w9kNfVBfW2JekynGcMayq0Qc1Ud/sUwcvVSKnPXl3w3Uc9Fsvzv1pGg5rGZOKxw/SQCX5RPBlPOVMIucPzS8YwB1xZczuM9XQ1os2aCOyr/W84CGxPO1oNGLWnn+gYgEUusXyQv1mE1z/XV8DwlvrXgrXp5vhOpPXwuFDBLIMh65JYFNmCvKM9ylQ4bMKsnymF8zQhqOGxTupIC1SVAgzr0FKb8YCiwbKUiUEe8qtm2FVOMkvXrwIHT/GpgUP2lDkPzjSFVqjzcat9LIU2iHmLEc4w2NYQT7Mcd95O1AVET+6ozMCO4vWkuvTXRgRe1WZ+HlURAT+Fd8Yzoh+axSmUF2DAsCkdyj17sgM1sa1qM2jpzfKjdcqPrJv2p6kc1oArmdyt7HtbsRm8Kdqu02iRIFmdC8gm7k9G5u40iewoFNjVw29D1YTLcLDOyYBvMwa9Ku0z0ISPmtufPpDDSJfBIZl3pgfLw7KUjiv+fIFFrF5fRVJjfK9mBFkoOfLOYOZdmy7Gt6PZgZ3iYieykE03Vdl1+59kDlwa6XsWdjkC5VXZ5V3gMeI5fnbtcTcd18k9Uxll5c4dJf685MxTPO0+oMVL7lWsrKDl3LREhy7ZDXABXfDmjzJ2bDhWu+mYWq29+P/WTNXueDJdv5zvB8H+68fur5gLEIkkdvBT8N7of8sLVpBrEYUvW6WtvVCLD1Mi8ZySiF7ubNqSLUHrZJ48B7X9dy1JJrGDKp/E8S7mlw6VdiS1OrhttrItaKwp2+UfWBSqtzfbQkNOwc2Hh0Ry1MDQWrSR+pfsCScE+IWIehrl/mzSfIGD1gxXUouQEup1mqu1rI+Hk85u6qWBPmxfK79zpS0CgwnneMsxl1JrXsXyUH0dxs+biUmjea1GFmu3kINB3nsy/XtEDoxcVFx94hxMVHPfnl2woWVte9djnDMCX229CqUXfXpPUCc+QIiHcS+i2ovb5GPHUu3XcUfmLflRSS2vHCU1Qg3XsMhN5AXZPLHTixI1ZRwG8c+t/Ero+UdECk54rk7giqHl6G8Nac6oSQoZDK8fRJar6xk55Qd3fMcsDgtpop6Lmsq9wp74TvFcGo04PTqLm4Mv5pohOMBBNtiskMWlxr1u4tdUKdsfU9857ef/3h7lS+zBPvwU88vWCJ//Wycr8rLiqLiPOI9uJhB+Wa5YEgnOBPl6y43wvbp7Vaj0OisDImvR4W3oujqSaLo8cT8HntcWrzVC/OKr42PKiXAkVTcwuX632QO3Qwmz8pHSUqf6b/5vfhZsGRrzf//MDm0yxpppf6IXqRdjJDORcPfQ75nGzrR5Nr5EaJrwUwLvwb/fYVAdaVTIvX7YqwJ7t+vfDcAHw+GgazC5QgxK1taxSVf/ltG/Nd+krx9bhWywVTIs5O0vLTuJXj1SZGxkZqT+oRYWWLvlAJr4zDZ3xMSFhqWn/+YLAKKospywk0ZPh3b3P67V6D6A1vdogHkbFY3HeiRLAOR5pULt+VoIBuL2XiZ309Eg/4PeSl9kMOUQY3fYNGXKw8mecQOwivghf9zeZNWMh9x2BsOIRDYV0Z9mebkaGU55BxYRhnMy/nl402Q7JKoYqReQRdWevn6Xyp+6bx0DDrrxGQ+bU1vc+2YuvuWnZxwihSvmUmppe3Yu8hRqW25jFqz24nZMi1QVO8nHZT/TievxF5OJbbhm1/Bz50NkxgY8jN0FKOduiJ3bx2uGMSvNn7HpEftsPi2UpXMmkKTIbXLvOEtqNI6OE05+w0Fj5q54b/73nTv+WGA5YzUNmqQddb046Rh3sOfyD+p4eq9nzX0Fk9zrFeOHs8XRBb+wBmdltnItW/t/aEiCv/DUBa+dFJ4sk+aG0oFLlGG4nEUUa/s4Z9OFJdZaDBEAf1zbybFrP6UB5G2cWy8zvGsrSHbrL3owTlExsY0xBzrtTU1LkfC5K0h8CWCUPmW2Me2wjo2XI7Ou9ANt7VMYUOhBPQ5UVOqxDIw6mGgZGvcgNY/KmWYHGxEOHnn+n6UfmfP94a5x1bhB89x9bUHxPRWi74/AlVHw+nFb/kdRFhK5lpEJlZCpVrt+IHVah7YPeEnKFYOWqShT/OdpBQS8fpUlBYAa/YR0tL98kQ/JB0yECiCoZx8V5Mf8jgeZqlMZhdjUivkHDJMApUIg6L+J6oRWPbqxvvqXuRs7G42DGgDqtnLLwaaYeQZmT53X2xHqq+kQxLxUM4ZT2SfU2+ES/WmrLL3YYhWLpT6ZNMEzIa/EeGdEewx3LpRBqhCV1DbYLjoqNgwAmjXc+akRs5fU6GewwkiapgqalmuFiubVXcNg6Nat/ABr0W2A7c+37u7zi+lw+kq5a3IHQxs4h5cQLD1QWB7ntboXczJap7fhJbBxI+qOq0ImD4POPThSkYzh8h5bu3YpD1/EavFWpeqP/UY0loRYK6UTTXmxnIv73LPJzeSu09IjLtAf89ZyyZI54grDVqFd4HjnnLEMt2J1zPKHB4lNuHW9mLNZs7/ZBvLnUscagPpQvlHwXvhSP5izPHbd5+1O55rLu5Pwaet2jtPrzox7yhZKrn7W+oU1xeU2/qxybv8+5uz1Jgd5iXa+oyAUfYv9j3nUwHxbHGz6ycAFbPWuGI2Uzcr4u6xak9AP/7qau+GdnY0KVxxI3KP7f1wfOS1nkQVHrMUKtLhJBT9s9Zqj5amW5YPyonYjO/clg31T9eoWqz9TkS5jr12qJJxXiY/bVzdzIJkislUz6SpaC7yXZM+QwZJ4NujKWElEG0WSNJCGTo0LMv982X4/qWYJ8qNQqyaxQky2Qrsf9Ju/amBQrSAnOcSZ+qcMRvzIEYNghRuVGpN53VcDy9uCR0awhbyzceS2aqRWzRs3dtO4ehe8C9X/9iHfYm/Tn+qn0Ycxl/OjS16yEqo9XZEDWCc6QHnRceN4A73Z3rjs0oUj1Mdn7Sp/ZbjvWkQM0xuK4dfWHN2wS7kb3a/RfH8Xska3XhQxMs87tpHx6bwJCgJvc+YhPKXMtNNDgmMWclatRxjNp/m+iNCnZOwZyHP0j0fjO4476Ue26ehniPupmDZzMEbN69T6WbgV757eT9cc3Q5RPal24/A08fjbnUlGawjQyEEi72/58eJHOH8zdwhE+5y/RDV/limP03B8zvbu0Pova38HiLsNjET3h6Q+Yps3s/ynN8FIqdvmDj+b+a61X9GGSf27uaEgWtMQOCFTsBRClWyd/hcTCeDToVYEnA403Thgs+Sdh+7nG/5RgBWjevfraUTgX94a5N4kYD2HlyS5VMYTpS7u1KdqDyr2vx7s5nkR/grE0xfatGhHKHrdw0VR/7W7URlUbEuU/NFv6n8nCEZtuuDXtJcD+lYhsdV4BjkROz8c4kLNIsZ5nvLEK+tef3r9R9+MSYZuNruWJkpW1jpXMkg5CQDifbEkSe4WsooPZDpU388W4ppXC94ZPeHE3B1n1Fgnf6y0CTLE0SkB5Ex+Sen2F0FcgN88zMGR+EOYFHselAJfSOcf+QCx7CoMLEa2nhKpxo4jGNVB3GkQbHs8cvVWPXM1nLV2wj4OB6IRwsWoOlsnRzM+IIbLsYR72P1OLC5M70PWmjONvSQVmjqYPb2emZPtcxTO1ldW+orEMMsy95v+E4JGeSZWit67HGd/jRYaUJeH+QWExgaUDQzH7bAyKTUBRSYN/2qQFZ9hJ8gnxTyJN3/hq41ICO7MY9d/ZM41Wc+77dUo1wNDkvHbttBs8/12h7Gjei5tDTkRSnGayELm0xfdOIweCTeh4t//TAUzeWfyHipop0bz96NY77nX1rD83fr90e/e7HL78Ma1q7j/BWdOGlPUDACWObgZ+Xg6C+0py0T5eAZ7+1S9bfRWC9iWe8IpGAEtJxf6nhGKxwbBZjo/b9cjve1Ur5BJQkV9DyvB3A2fnPzg9SvlN7qvs9Wyr/YfInCrXcUmEeTNx77zoRNBPTu95T9VGYZz8mH06E5kpjmPDPTLDm8tZf/kOEWZx6QBFvNoL3m+xRuUeCVuSem1dVcnHN76+1QwsJHK3S8t9t81Fv4N/dLU/GhYqloLjYQghk6bGo1pExt03WVfl+EU5MBbmMq1Mg2fYhTs+kGNyl/tEe4xR8OHuHztK6BPJBJw8e+jCIwcf5b+69KoX6CcJDlxNDaFceeN1mXgbXk8ppbj1DkNMu2nX7UTkEf9e/GPcehlNTsNt9xQq8P8UTqKo0gu3XiE9vHavEx1/pUo92j2Knrj2r6UIljOk83yeSqXnBtDyyMaMK2fqZFUE/xjDXODN8U68aUe8GzX29x7G//s/bb8vVuOFBYzH6fAI6cvVCr21rULPaydV6axI0+Z/9VaZqcFGAbfqt2BQmrqVvmLpSC0mx73Pbjkzjk/9XqyuOtSBacpU4s8xgeZ6W7JdQi5PHFeYOucxgzWP3mSdZtZgPPtPVTO3r/5cXjZPimffE6D2o/F4vaTpiE2mLEwYH7hZZEbD/rk6ZR5QnjrxbUDCJJmBX55Y/P1f98Xp7nUcciQDfEMl1i6+hoM97NWQiOAB+mg07L/JH4a3PqEev9wC0fa6FG9F+Rd+rygobKv+FYxanWLoTkKBHNHosTYRgGk1LHlUfX7tlj1/yJCJcNkhrgOof8yVcpQxkIrYR2r18BdKxe9TRr1yChMPXFxX62TOhFfO4+m00Cd+mdVKm/v6AQyyflvAuMkTLiq7WkLNxu6Rz07ArGTYmCdZXSnIh4tXA/34btQ+UCJQsBudDKP/Zs01+FIT/aomUNS6EvOnAYxW+Qag3vzJM3leEYOHoSzIlg5Acz58KWCmCTcjOfVkPhkBiNYvJaymGxrWGnGfMw7Dz5WMgRpbAtPbo5gPlw3iyV3si40kp9Brtum3tRnB2Y+/G7sNlWFnIaeORHsURnm+vB9rL0OL87av39jHQhqrVqVuVI8Zgxk+mbwxO5TpRTdsqkPJ8m0V86jiOSIvGj3tWIIxSe2Kn+wTMUgLqNWgq0fj8c2qB4SQ+MlSVjulWwuZCSdjk9SmMpbaeC4qvRMa9KOXK09N4waV75tNAJfxmfEye75nBTjaL28I0VSj8aebF4T6Dd7YCmrzbqjB/4P7piMp/epDM7Bc3dDQruttEwF3jIWulKzboUUorCRknoEZE/+PaFVfcMNc6pMM2gAqG8ksfMj5h+xfSqehbA4gYdAi/LRsEg+VzAm8iB6j5wHLexigc53TcnltR+a/8zPfHIysKq9V/KY6iRJRXpZnHUPUhXzx6/Al1D7KRbP1vKiSA/3re9LFKImLP3Tad+ZQEbp2krzUcJMgVRj9QNkrB+ZCzP69bksCpIWzMuzMNXTcSspL7SOg5ZlSa/T0dEVeP3pqm+oFCSpe7sEIm3vEeclkvJmPF/94bX8oP6FQ6WbZdpmCI/U5xqEU2nmyZf6pRQ8HOqIjMXxtyoRb27Zuh9iDW93KE0H3IQ1bCxOmR+UGcqTBxZKYtAHPn0rNg/yGEHzcJO/e8EDXT1XTyksMwv9AvQI1/9DXZpMT9HIZW81auNK8icJ82u3g3YQSut3i8bQSKUR/ok8FgOIpZg3j6/pJilD8TbToqOIaOdTq2BIUSHHabLri0RtXHKw/pH5VUfV18KfmDuiek+AhrxHOlOPw3Snw0agKCzFlr+/xKobk2nHXQdhJ3xUtyTQdLwZ2Rdy1BewoXD0r5MR8sg0yb/2qa2DQcNd+qX1cuA4P3R7dnXDMoLHXOjXpaBqG7X38Gec5ApE8nqtqsDD8raVXn+f79n1rRg5Z8WqOKsQzq+X7r9uVy1qA5hjmF/RblByBUFdHPMeQE3ayXtLCg5v2GA2+vPPPG27VUvZ0ZA9AzcpHTG/bD4Hl3eksq/z/ei5tV6wXDUI1gk3KaiPx9Xf3BVH2IzurQZBkSIcyJch+qf3yuTNjhkUhEi0NAG4maL+IOvvuFFolY5xi1WYiIg/ZK4d5EeRKaJsNFd2h/Q8bvixFzsSR0bBMYO878Hbaz4jc2bCeDk8eL31M/BZfNHjY3W5OReoUb6VmpUFV6zHhrloxyTtu6o1vTEabcwWz6nIKAl8leExoZENIp28S6SMF80ZYv7WGZ8LLI4WWj5oPmOcpoMPEH1vbstLvOPQSTk+71jXuy0bv1DpM9hvBwUESYJJOD6128HqaPhnGXsZSZ90kulEXU3rdS+4KuoxVf7Zs8tLyM+3StfgRPxD9rb3DIxxHymR/PXEdR8LT6LtmiAJGehxYrbowh5LpQloNGIewmY5Ta2MYx9X0vgeEgoNuuNrhhaBzK2ycVwjqBiNErIelZE+Die2XZKFKET0J57GKek6Bl3NFRZFOE0LepDJRHUzi5pHNqLakIqdt+d09ITyO3/Omj5ZYiuCnXlHkemMFRUWxTHStCxh+CQ5L3DBLFxB0aZovwgtFkosLiPz1cKxLP8tEPcKbme42mxBdZRxMQx76cbAgcgFQgaZuNoz2W/vy+GVg+gPD5LzEHZF3xNiHwkgWV/7r4+UtbnD5i/EyGytBRIuZM82cc/6ePJ1eV990mQnVyw4M8qn9YFxlK7/UhwpgocO5sXCiO3xR7NNBKxLQm++EQhwhYkffFPNtPAu9Kz+Gx1CiI5X38W2BCwmj8fNwdQgyiTXrJ9fUkhJ1lcT6+OQ4WaqsvPlH3gmnPXi01gQTsP/zwK204GWxO+uQEpUTYHhV2YGCjnv8v205UGH4Hu3gAs6svBS1nVEo2XkpBkz2HpdE+aj+4ISLbSM2nksm8ny4Jg/BmZ7+VJ5UGKTv/3uTLQ9gzSP8sVjQdjOlvM2OIQ/A/wN+gyp+BvaU/L/N9GMaGCq+rZkyZqLXIzBkUGkENqxxr+EQmGm7YHno3NgIrjzsTD3J/oGIfA9E/ehQlV07PKrzKwm3zjXpJD8Yw8141u/NENvgcnASlT4yD2beU1b8+G8bNnlLVi+PY9y3z5UGdHEwi3jm6cgIXMj56CXXm4P6LgSGB4EmUbv+xK+xiLrQMdhlnm0zhO1m/Q98zF9P3k+ec5aexK61QOrEqF4VJTi2kQzOYZ2gMHJzIhS3xwUlNnxmI01wPT1zMxffcCW+2jH960DP8mU9jczi3K38ALLekRetkDEGza9qBu3MAYfTHs3/L2MBwxYVsTuW/W9HK4U6mA1aeSN7k5iGiwS+8TZqqj/n4DH93WSL4t/CpmVP9o290y1SBJRFFTH95xqn5oqP4fiQthQjR+8HNIWv+OENzl8Fwjgiv/Bvsf6WCoKvZIkwUJ0HoKn+BwocvoDcmX2LyJsFr5/PJK/fDQZBT5pybIME4Z+1U4YZIOA2RpGyUyViP4Y9TDo2Cp7cT4XMBGbGf4pp3C8fA7o370oWzFCQ/zY1dr4qFVdFawKVUCj5Pacl3qsXhw1Jft+v5QQws8HwQ742H6JZbm36WD2JP7o7rsRrf4GpjknxJZwgNEi3K96sToS/rWSu6MgRzvfhNkQLfQfQf/h0cNgwlzrFVQ6dkaMwx5jFfH0HW6XrNrgcpKAs8vPPS+gjW3Aro49dTwF/cPJObPQqLh4y7WbxS4d/1Q4zJagyvL1go/2VJwxoTS2WB+DiC9Xfd03VNg0pgcrowwwRcJhSHBObTMD5/hYnSOoGzVssBtkrpeCTqV8QTM4k3+/q9zYLTcdo3LfGK9RTWift2q3emQzzF4EyI0jSW/X73um3IwKEV7SZjvhm8V4qyjOXMQIKYrVPYpxloztEn6B/OgCKnFid54J8e7IyJ+T8kpIyjRwfAlfcxN9JRH9Y+kUX/4389Omj/d6pfiB/vr9NmJ2Iy4dABRlkbxNGeM9soTESfVcneCKp/sEQ88X1+hwi6hDHZF9R8Gd1+zD3Eg4jLgYUpn6+6ou0iF8WlggjHTV2ZHtGe2FAVQBDaTMLmUxTxMup+iXtwlC5QhYTb1hVySPqEeKE4h+QYEtRUv7u5d/lhiT8g1Ji6H421ffnUaT5j0y62iaoH1L3wgXXHnkNBkDugFJDYQgbxenBxiFQIpBSLgzbIURD9gabfRisUpOgvA0nlFJxpeJG15Vw4XL6fe+SqMIi2raOoOh+B07Qaxs+6B7FZh9bUXCgSOaLKAdzPh/DNevDaz2NRaOQ5zfNw+zAkWx7qMLBHg73V+8RSyjBi5MvkDVeiUef/nd9BZwTRBh0CJS0xcDFj0IlnHEUfzRP7r6GxEDl9x4W2chRpMUKzSTpf0c5aPsbmNIauX3u3GW6LQ2wgj/QTuXEYe0hEWMXHQcL3+DUjpgk8G8sKeH8uHr5O+yQ8+iYgRaj6rfk9HlIFDu3HEichnFba2caSgL8mFhrBtlOwC66cYzBIQLEKzfqZW9Ow4iDe5A9OQESOrLv4MWo+HFddTipIwAXXv8vvfWdQR1HjE6lOQMiGjn3HVv7pYfRbd/6BOaWc//Gfn2PWn+B4EyzFtmKBjERkfN9qm0jVh0hj+ta3h4gwVHWV/UX1D55KroUj1H1Afjfw1Iqql1dBTObe1H7AoZDXsWfIHCqRxQa5QUTsaOU7vELtp6+NEmIjqfnwdue7365RtqBvq6++upuEi/HslHXqvrVIlbMOvUvCgGhjX3SiA+b6uhai00jg2ppGZOh0wmsbbStFRjLkKxMFFf9+AA9XmqW7KRkMgfSHzA+4Qnd242+FATJSohc3XL3kjq7khSQzVQo0VTRCbql7giC+68ZoM4WaP/fIdme84U/73s9eYxBlpEvWBUIfId8057VraBCKO4ZZewV90NbMXvjYZgizTGrC+Yc/wexKiLYU+zBsNh8zOsHiCx+vHekaecN4JNmmTJn3hftFa0+1hyNYuGFm5Fznh9A3YrX97KMgT1ccyg7wx+6mZ3s+NI7i20q35Ra1AOyxrBzucx+DM3NwOCvdZ2jO9nrdVRoH2zrllk34Z8hLLGU7sk1A3/+excTxQDxJvTqfRZqAbl0O7enoQKgfeRCunTIJk2Ry3dnNQUjVUt1YYz+F8Cfpg3QaQVC3ro1UV5vGg08k7gGvIAyopxco8c+gNFk/VDgtCH2mEhNMfjNUfbFPpxUEYfjd1v/j/396sDrxK5+OP4zSQU/1/3EGKv+XYWbwSneFhQjaywPi8VR9pA0IMxOp/eDZM+fTZBktNEwFpX2g5sOLAO9ofape8k/3BcwaU/diOfOu8cEHCAt6+GJvGBHaShD7esUQkdfvTq93EtEmvOEOfZQxyDQ8OuF7STAdKFQJfmsCycb93L/vkbCwnESzO/ElphMWjy78IIHi1Hf0coc5JKyjNT6zkFFIfMUt+8cS/y9QAa/+gz5fxsx24r9bA9txICB6v4BoMBGA4+K/lFmdzJU9fL/wjMiNmVjjv4ElU508VH6/dEb2fPrV479YIz3POjGAv7HByTh9W+S/ZWkuQDYzgb/yqj9k9Ojkv/xFHXEML4K/MFISIip+5b+gyJtmoCODv/oZBVnfGua/sVpdbr4PhL8ALpMMy77mv/hdBhsf8oS/hixO0Zlp5784jywxa8mFv7OhWWTtGui//O19hUCUhr880jttXNLovyCCPsA3UYe/4WSlcHKP6b/o9f7s6v6Hv+zdzvmvUeq/vmw1r/ybiL/RQpr+ihjrv9O8deYfJ4m/lX7PgG/j67+epTeAH5+Jv/2Jm2zAsey/jjxcK+YCir/0aA6y2ILtvy+EaZuFUYq/+sTNlAxW7r8/Z2MIPYqKvy9suSqrKu+/lk2PnH6sir8S7vm+PpXvv+gibovzt4q/sjvsYg==AQAAAAAAAAAAgAAAAAAAABhmAAAAAAAAswcAAAAAAAA=eF593DkOq9gCRdGa2ZtqDYGQkMABAQFCCCGEEH1nMAyhIhL+W39nXnZ65A7uP//8v/7983cP4CE8gn/gMTyBp/AMnsMLeAmv4DW8gbfwDt7DB/gIn+AzfIGv8A2+w7/wA37Cf/ALfsOf3jt4Hr89gIfwCP6Bx/AEnsIzeA4v4CW8gtfwBt7CO3gPH+AjfILP8AW+wjf4Dv/CD/gJ/8Ev+A1/er8PPK97ewAP4RH8A4/hCTyFZ/AcXsBLeAWv4Q28hXfwHj7AR/gEn+ELfIVv8B3+hR/wE/6DX/Ab/vT+HPTs4u0BPIRH8A88hifwFJ7Bc3gBL+EVvIY38BbewXv4AB/hE3yGL/AVvsF3+Bd+wE/4D37Bb/hT9PJnF28P4CE8gn/gMTyBp/AMnsMLeAmv4DW8gbfwDt7DB/gIn+AzfIGv8A2+w7/wA37Cf/ALfsOf3t+D//3zdw/gITyCf+AxPIGn8Ayewwt4Ca/gNbyBt/AO3sMH+Aif4DN8ga/wDb7Dv/ADfsJ/8At+w5/evwM9u3h7AA/hEfwDj+EJPIVn8BxewEt4Ba/hDbyFd/AePsBH+ASf4Qt8hW/wHf6FH/AT/oNf8Bv+9P4d9NnF2wN4CI/gH3gMT+ApPIPn8AJewit4DW/gLbyD9/ABPsIn+Axf4Ct8g+/wL/yAn/Af/ILf8Kf05c8u3h7AQ3gE/8BjeAJP4Rk8hxfwEl7Ba3gDb+EdvIcP8BE+wWf4Al/hG3yHf+EH/IT/4Bf8hj+9/wd7dvH2AB7CI/gHHsMTeArP4Dm8gJfwCl7DG3gL7+A9fICP8Ak+wxf4Ct/gO/wLP+An/Ae/4Df86f0/8LOLtwfwEB7BP/AYnsBTeAbP4QW8hFfwGt7AW3gH7+EDfIRP8Bm+wFf4Bt/hX/gBP+E/+AW/4U/v6yCeXbw9gIfwCP6Bx/AEnsIzeA4v4CW8gtfwBt7CO3gPH+AjfILP8AW+wjf4Dv/CD/gJ/8Ev+A1/Kl/+7OLtATyER/APPIYn8BSewXN4AS/hFbyGN/AW3sF7+AAf4RN8hi/wFb7Bd/gXfsBP+A9+wW/40/s6uGcXbw/gITyCf+AxPIGn8Ayewwt4Ca/gNbyBt/AO3sMH+Aif4DN8ga/wDb7Dv/ADfsJ/8At+w5/e14E+u3h7AA/hEfwDj+EJPIVn8BxewEt4Ba/hDbyFd/AePsBH+ASf4Qt8hW/wHf6FH/AT/oNf8Bv+9L4O+tnF2wN4CI/gH3gMT+ApPIPn8AJewit4DW/gLbyD9/ABPsIn+Axf4Ct8g+/wL/yAn/Af/ILf8Kf25c8u3h7AQ3gE/8BjeAJP4Rk8hxfwEl7Ba3gDb+EdvIcP8BE+wWf4Al/hG3yHf+EH/IT/4Bf8hj+974P598/fPYCH8Aj+gcfwBJ7CM3gOL+AlvILX8Abewjt4Dx/gI3yCz/AFvsI3+A7/wg/4Cf/BL/gNf3rfB/bs4u0BPIRH8A88hifwFJ7Bc3gBL+EVvIY38BbewXv4AB/hE3yGL/AVvsF3+Bd+wE/4D37Bb/jT+z7IZxdvD+AhPIJ/4DE8gafwDJ7DC3gJr+A1vIG38A7ewwf4CJ/gM3yBr/ANvsO/8AN+wn/wC37Dn8aXP7t4ewAP4RH8A4/hCTyFZ/AcXsBLeAWv4Q28hXfwHj7AR/gEn+ELfIVv8B3+hR/wE/6DX/Ab/vS+D/7ZxdsDeAiP4B94DE/gKTyD5/ACXsIreA1v4C28g/fwAT7CJ/gMX+ArfIPv8C/8gJ/wH/yC3/Cn9zkQzy7eHsBDeAT/wGN4Ak/hGTyHF/ASXsFreANv4R28hw/wET7BZ/gCX+EbfId/4Qf8hP/gF/yGP73PQXl28fYAHsIj+AcewxN4Cs/gObyAl/AKXsMbeAvv4D18gI/wCT7DF/gK3+A7/As/4Cf8B7/gN/xpffmzi7cH8BAewT/wGJ7AU3gGz+EFvIRX8BrewFt4B+/hA3yET/AZvsBX+Abf4V/4AT/hP/gFv+FP73Ownl28PYCH8Aj+gcfwBJ7CM3gOL+AlvILX8Abewjt4Dx/gI3yCz/AFvsI3+A7/wg/4Cf/BL/gNf3qfA/fs4u0BPIRH8A88hifwFJ7Bc3gBL+EVvIY38BbewXv4AB/hE3yGL/AVvsF3+Bd+wE/4D37Bb/jT+xzEZxdvD+AhPIJ/4DE8gafwDJ7DC3gJr+A1vIG38A7ewwf4CJ/gM3yBr/ANvsO/8AN+wn/wC37Dn46XP7t4ewAP4RH8A4/hCTyFZ/AcXsBLeAWv4Q28hXfwHj7AR/gEn+ELfIVv8B3+hR/wE/6DX/Ab/vQ+B/fZxdsDeAiP4B94DE/gKTyD5/ACXsIreA1v4C28g/fwAT7CJ/gMX+ArfIPv8C/8gJ/wH/yC3/Cn9znQzy7eHsBDeAT/wGN4Ak/hGTyHF/ASXsFreANv4R28hw/wET7BZ/gCX+EbfId/4Qf8hP/gF/yGP73PQX928fYAHsIj+AcewxN4Cs/gObyAl/AKXsMbeAvv4D18gI/wCT7DF/gK3+A7/As/4Cf8B7/gN/z9/NOzi7cH8BAewT/wGJ7AU3gGz+EFvIRX8BrewFt4B+/hA3yET/AZvsBX+Abf4V/4AT/hP/gFv//H/wPNRQ7IAQAAAAAAAAAAgAAAAAAAAAAAAAAAAAAAqRMAAAAAAAA=eF5d3UUURlXBRmG/oLu7u7vTTixUTMouREXFBAW7BVsUsLABOwFFscDuALGxsZN/8N9ncPadvGvtu6d7dNa952Y3+/9nNu22024TPp92u3j4Ytrt4+HLaXeIh68w7Y7x8BWn3SkevtK0O8fDV552l3j4KtPuGg9fddrd4uGrTbt7PHz1afeIh68x7Z7x8DWn3Ssevta0e8fD1552n3j4OtPuGw9fd9r94uHrTbt/PHz9aQ+Ih28w7YHx8A2nPSgevtG0B8fDN572kHj4JtMeGg/fdNrD4uGbTXt4PHzzaY+Ih28x7ZHx8C2nPSoevtW0R8fDt5725vF0rudbTnuLcD3fKh6u51vHw/V8m3i4nm8bD9fz7eLher59PFzPd4iH6/mO8XA93ykerudj4uF6vnM8XM93iYfr+a7xcD3fLR6u57vHw/V8bDxcz/eIh+v5nvFwPd8rHq7n4+Lher53PFzP94mH6/m+8XA93y8eruf7x8P1/IB4uJ6Pj4fr+YR4uJ5PjIcfPe1J8XDdnxxP53p+0LQPDNfzg+Phen5IPFzPD42H6/lh8XA9PzwerudHxMP1/Mh4uJ4fFQ/X86Pj4Xo+JR6u58fEw/V8ajxcz4+Nh+v5cfFwPT8+Hq7n0+Lhen5CPFzPT4yH6/lJ8XA9nx4P1/OT4+F6fko8XM9PjYfr+WnxcD0/PR6u52fEw/V8Rjxcz2fGw0+c9pnxcD0/Kx6u+7Pi6VzPz5727HA9PycerufnxsP1/Lx4uJ6fHw/X8wvi4Xp+YTxczy+Kh+v5xfFwPb8kHq7nl8bD9fyyeLieXx4P1/Mr4uF6Picerudz4+F6fmU8XM+viofr+dXxcD2/Jh6u59fGw/X8uni4nl8fD9fzG+Lhej4vHq7nN8bD9fymeLiez4+HnzntBfFwPV8YD9fzm+Phun9LPJ3r+W3TvjVcz2+Ph+v5oni4nt8RD9fzO+Phen5XPFzP746H6/k98XA9vzceruf3xcP1fHE8XM+XxMP1fGk8XM/vj4fr+QPxcD1/MB6u5w/Fw/X84Xi4nj8SD9fzR+Phev5YPFzPH4+H6/kT8XA9fzIerudPxcP1fFk8XM+Xx8PPn/aKeLiePx0P1/Nn4uF6vjIervvPxtO5nq+a9nPhev58PFzPX4iH6/mL8XA9fykerucvx8P1fHU8XM/XxMP1/JV4uJ6/Gg/X89fi4Xr+ejxcz9+Ih+v5m/FwPX8rHq7nb8fD9fydeLievxsP1/P34uF6/n48XM8/iIfr+YfxcD3/KB6u5x/Hw/V8bTxcz9fFwy+f9ifxcD1fHw/X80/j4Xr+WTxczz+Ph+v+F/F0rudfTfvLcD3/Oh6u5xvi4Xr+TTxcz7+Nh+v5d/FwPf8+Hq7nP8TD9fzHeLie/xQP1/ON8XA9/zkerue/xMP1/Nd4uJ7/Fg/X89/j4Xr+Rzxcz/+Mh+v5X/FwPf87Hq7n/8TD9fzfeLie/xcP1/NN8XA9O/Di4dd57UBsNnI9z2ejh+t5MRs9XM/L2ejhel5hNnq4nlecjR6u+5Vmo6dzPa8y8ZVnI9fzqrPRw/W82mz0cD2vPhs9XM9rzEYP1/Oas9HD9bzWbPRwPa89Gz1cz+vMRg/X87qz0cP1vN5s9HA9rz8bPVzPG8xGD9fzhrPRw/W80Wz0cD1vPBs9XM+bzEYP1/Oms9HD9bzZbPRwPW8+Gz1cz1vMRg/X85az0cP1vNVs9HA9bz0bPVzP28xGD9fztvFwPW8XD9fz9vFwPe8QD9fzjvFwPe8UD9f9zvF0ruddp90lXM+7xcP1vHs8XM97xMP1vGc8XM97xcP1vHc8XM/7xMP1vG88XM/7xcP1vH88XM8HxMP1fGA8XM8HxcP1fHA8XM+HxMP1fGg8XM+HxcP1fHg8XM9HxMP1fGQ8XM9HxcP1fHQ8XM83j4fr+RbxcD3fMh6u51vFw/V863i4nm8TD9fzbePher5dPFz3t4+ncz3fcdo7hOv5TvFwPR8TD9fznePher5LPFzPd42H6/lu8XA93z0erudj4+F6vkc8XM/3jIfr+V7xcD0fFw/X873j4Xq+Tzxcz/eNh+v5fvFwPd8/Hq7nB8TD9Xx8PFzPJ8TD9XxiPFzPJ8XD9XxyPFzPD4yH6/lB8XA9PzgerueHxMP1/NB4uJ4fFg/X88Pj4bp/RDyd6/lR0z4yXM+Pjofr+ZR4uJ4fEw/X86nxcD0/Nh6u58fFw/X8+Hi4nk+Lh+v5CfFwPT8xHq7nJ8XD9Xx6PFzPT46H6/kp8XA9PzUeruenxcP1/PR4uJ6fEQ/X8xnxcD2fGQ/X8zPj4Xp+Vjxcz2fFw/V8djxcz8+Oh+v5OfFwPT83Hq7n58XD9fz8eLieXxAP1/0L4+lczy+e9kXhen5JPFzPL42H6/ll8XA9vzwerudXxMP1fE48XM/nxsP1/Mp4uJ5fFQ/X86vj4Xp+TTxcz6+Nh+v5dfFwPb8+Hq7nN8TD9XxePFzPb4yH6/lN8XA9nx8P1/MF8XA9XxgP1/Ob4+F6fks8XM9vjYfr+W3xcD2/PR6u54vi4Xp+Rzxcz++Mh+v5XfFw3b87ns71/N5p3xOu5/fFw/V8cTxcz5fEw/V8aTxcz++Ph+v5A/FwPX8wHq7nD8XD9fzheLiePxIP1/NH4+F6/lg8XM8fj4fr+RPxcD1/Mh6u50/Fw/V8WTxcz5fHw/V8RTxcz5+Oh+v5M/FwPV8ZD9fzZ+Phev5cPFzPV8XD9fz5eLievxAP1/MX4+F6/lI8XM9fjofr/up4OtfzV6a9JlzPX42H6/lr8XA9fz0erudvxMP1/M14uJ6/FQ/X87fj4Xr+Tjxcz9+Nh+v5e/FwPX8/Hq7nH8TD9fzDeLiefxQP1/OP4+F6vjYerufr4uF6/kk8XM/Xx8P1/NN4uJ5/Fg/X88/j4Xr+RTxcz7+Mh+v5V/FwPf86Hq7nG+Lhev5NPFzPv42H6/l38XDd/z6ezvX8x2n/EK7nP8XD9XxjPFzPf46H6/kv8XA9/zUerue/xcP1/Pd4uJ7/EQ/X8z/j4Xr+Vzxcz/+Oh+v5P/FwPf83Hq7n/8XD9XxTPFzPPui8KVzPMx98zkeu5/l89HA9L+ajh+t5OR89XM8rzEcP1/OK89HD9bzSfPRwPa88Hz1cz6vMRw/X86rz0cP1vNp89HA9rz4fPVzPa8xHD9fzmvPRw3W/1nz0dK7ndSa+9nzkel53Pnq4ntebjx6u5/Xno4freYP56OF63nA+erieN5qPHq7njeejh+t5k/no4XredD56uJ43m48erufN56OH63mL+ejhet5yPnq4nreajx6u563no4freZv56OF63jYeruft4uF63j4erucd4uF63jEerued4uF63jkerudd4uF63jUerufd4uF63j0eruc94uF63jMerue94uG63zuezvW877T7hOt5v3i4nvePh+v5gHi4ng+Mh+v5oHi4ng+Oh+v5kHi4ng+Nh+v5sHi4ng+Ph+v5iHi4no+Mh+v5qHi4no+Oh+v55vFwPd8iHq7nW8bD9XyreLiebx0P1/Nt4uF6vm08XM+3i4fr+fbxcD3fIR6u5zvGw/V8p3i4no+Jh+v5zvFwPd8lHq7nu8bDdX+3eDrX87HT3j1cz/eIh+v5nvFwPd8rHq7n4+Lher53PFzP94mH6/m+8XA93y8eruf7x8P1/IB4uJ6Pj4fr+YR4uJ5PjIfr+aR4uJ5Pjofr+YHxcD0/KB6u5wfHw/X8kHi4nh8aD9fzw+Lhen54PFzPj4iH6/mR8XA9PyoerudHx8P1fEo8XM+PiYfr+dR4uJ4fGw/X/ePi6VzPp037+HA9PyEerucnxsP1/KR4uJ5Pj4fr+cnxcD0/JR6u56fGw/X8tHi4np8eD9fzM+Lhej4jHq7nM+Phen5mPFzPz4qH6/mseLiez46H6/nZ8XA9PycerufnxsP1/Lx4uJ6fHw/X8wvi4Xp+YTxczy+Kh+v5xfFwPb8kHq7nl8bD9fyyeLieXx4P1/Mr4uG6PyeezvX8ymnPDdfzq+Lhen51PFzPr4mH6/m18XA9vy4erufXx8P1/IZ4uJ7Pi4fr+Y3xcD2/KR6u5/Pj4Xq+IB6u5wvj4Xp+czxcz2+Jh+v5rfFwPb8tHq7nt8fD9XxRPFzP74iH6/md8XA9vyserud3x8P1/J54uJ7fGw/X8/vi4Xq+OB6u50vi4Xq+NB6u5/fHw3X/gXg61/OHpv1guJ4/HA/X80fi4Xr+aDxczx+Lh+v54/FwPX8iHq7nT8bD9fypeLieL4uH6/nyeLier4iH6/nT8XA9fyYerucr4+F6/mw8XM+fi4fr+ap4uJ4/Hw/X8xfi4Xr+Yjxcz1+Kh+v5y/FwPV8dD9fzNfFwPX8lHq7nr8bD9fy1eLievx4P1/M34uF6/mY8XPffiqdzPX9n2m+H6/m78XA9fy8erufvx8P1/IN4uJ5/GA/X84/i4Xr+cTxcz9fGw/V8XTxczz+Jh+v5+ni4nn8aD9fzz+Lhev55PFzPv4iH6/mX8XA9/yoerudfx8P1fEM8XM+/iYfr+bfxcD3/Lh6u59/Hw/X8h3i4nv8YD9fzn+Lher4xHq7nP8fD9fyXeLie/xoP1/3f4ulcz/+Y9u/hev5nPFzP/4qH6/nf8XA9/yceruf/xsP1/L94uJ5viofr2Q+LbwrX88wPjRcj1/N8MXq4nheL0cP1vFyMHq7nFRajh+t5xcXo4XpeaTF6uJ5XXowerudVFqOH63nVxejhel5tMXq4nldfjB6u5zUWo4frec3F6OF6Xmsxerie116MHq7ndRajh+t53cXo4XpebzF6uJ7XX4werucNFqOH63nDxejhut9oMXo61/MmE994MXI9b7oYPVzPmy1GD9fz5ovRw/W8xWL0cD1vuRg9XM9bLUYP1/PWi9HD9bzNYvRwPW8bD9fzdvFwPW8fD9fzDvFwPe8YD9fzTvFwPe8cD9fzLvFwPe8aD9fzbvFwPe8eD9fzHvFwPe8ZD9fzXvFwPe8dD9fzPvFwPe8bD9fzfvFwPe8fD9fzAfFwPR8YD9fzQfFw3R8cT+d6PnTaQ8L1fFg8XM+Hx8P1fEQ8XM9HxsP1fFQ8XM9Hx8P1fPN4uJ5vEQ/X8y3j4Xq+VTxcz7eOh+v5NvFwPd82Hq7n28XD9Xz7eLie7xAP1/Md4+F6vlM8XM/HxMP1fOd4uJ7vEg/X813j4Xq+Wzxcz3ePh+v52Hi4nu8RD9fzPePher5XPFzPx8XD9XzveLju7xNP53q+37T3Ddfz/ePhen5APFzPx8fD9XxCPFzPJ8bD9XxSPFzPJ8fD9fzAeLieHxQP1/OD4+F6fkg8XM8PjYfr+WHxcD0/PB6u50fEw/X8yHi4nh8VD9fzo+Phej4lHq7nx8TD9XxqPFzPj42H6/lx8XA9Pz4erufT4uF6fkI8XM9PjIfr+UnxcD2fHg/X85Pj4bp/Sjyd6/lp0z41XM9Pj4fr+RnxcD2fEQ/X85nxcD0/Mx6u52fFw/V8Vjxcz2fHw/X87Hi4np8TD9fzc+Phen5ePFzPz4+H6/kF8XA9vzAerucXxcP1/OJ4uJ5fEg/X80vj4Xp+WTxczy+Ph+v5FfFwPZ8TD9fzufFwPb8yHq7nV8XD9fzqeLieXxMP1/Nr4+F6fl08XPevj6dzPZ837RvC9fzGeLie3xQP1/P58XA9XxAP1/OF8XA9vzkerue3xMP1/NZ4uJ7fFg/X89vj4Xq+KB6u53fEw/X8zni4nt8VD9fzu+Phen5PPFzP742H6/l98XA9XxwP1/Ml8XA9XxoP1/P74+F6/kA8XM8fjIfr+UPxcD1/OB6u54/Ew/X80Xi4nj8WD9fzx+Phuv9EPJ3r+VPTfjJcz5fFw/V8eTxcz1fEw/X86Xi4nj8TD9fzlfFwPX82Hq7nz8XD9XxVPFzPn4+H6/kL8XA9fzEerucvxcP1/OV4uJ6vjofr+Zp4uJ6/Eg/X81fj4Xr+Wjxcz1+Ph+v5G/FwPX8zHq7nb8XD9fzteLievxMP1/N34+F6/l48XM/fj4fr+QfxcD3/MB6u+x/F07mer532x+F6vi4eruefxMP1fH08XM8/jYfr+WfxcD3/PB6u51/Ew/X8y3i4nn8VD9fzr+Pher4hHq7n38TD9fzbeLiefxcP1/Pv4+F6/kM8XM9/jIfr+U/xcD3fGA/X85/j4Xr+Szxcz3+Nh+v5b/FwPf89Hq7nf8TD9fzPeLie/xUP1/O/4+F6/k88XM//jYfr/n/xdK5nF/LeFK7nmQt7lyPX83w5erieF8vRw/W8XI4erucVlqOH63nF5ejhel5pOXq4nldejh6u51WWo4fredXl6OF6Xm05erieV1+OHq7nNZajh+t5zeXo4Xpeazl6uJ7XXo4erud1lqOH63nd5ejhel5vOXq4ntdfjh6u5w2Wo4frecPl6OF63mg5erieN16OHq7nTZajh+t50+Xo4XrebDl6uJ43X44eructlqOH63nL5ejhut9qOXo61/M20269HLmet42H63m7eLiet4+H63mHeLied4yH63mneLied46H63mXeLied42H63m3eLied4+H63mPeLie94yH63mveLie946H63mfeLie942H63m/eLie94+H6/mAeLieD4yH6/mgeLieD46H6/mQeLieD42H6/mweLieD4+H6/mIeLiej4yH6/moeLjuj473f/HX4qI=AQAAAAAAAAAAgAAAAAAAAAAgAAAAAAAAdAUAAAAAAAA=eF412RGjAooBBeC7rW1ta3tBEARBEARBEARBEARBEARBEARBEARBEFwIgiAIgiAIgiAIgiAIgiAIggtBEATB4H7vyPcbzjmBr98EGWKYEUYZY5wJJplimhlmmWOeBRZZYpkVVlljnQ022WKbHXbZY58DDvnNEceccMoZ51xwyRXX3HDLHfc88MgTz7zwyhvv/OGDT7745odff/klwCBDDDPCKGOMM8EkU0wzwyxzzLPAIksss8Iqa6yzwSZbbLPDLnvsc8AhvznimBNOOeOcCy654pobbrnjngceeeKZF155450/fPDJF9/88OuvvwQYZIhhRhhljHEmmGSKaWaYZY55FlhkiWVWWGWNdTbYZIttdthlj30OOOQ3RxxzwilnnHPBJVdcc8Mtd9zzwCNPPPPCK2+884cPPvnimx9+/e2XAIMMMcwIo4wxzgSTTDHNDLPMMc8CiyyxzAqrrLHOBptssc0Ou+yxzwGH/OaIY0445YxzLrjkimtuuOWOex545IlnXnjljXf+8MEnX3zzw68AGGSIYUYYZYxxJphkimlmmGWOeRZYZIllVlhljXU22GSLbXbYZY99DjjkN0ccc8IpZ5xzwSVXXHPDLXfc88AjTzzzwitvvPOHDz754psffv39lwCDDDHMCKOMMc4Ek0wxzQyzzDHPAossscwKq6yxzgabbLHNDrvssc8Bh/zmiGNOOOWMcy645IprbrjljnseeOSJZ1545Y13/vDBJ19888Ovf/wSYJAhhhlhlDHGmWCSKaaZYZY55llgkSWWWWGVNdbZYJMtttlhlz32OeCQ3xxxzAmnnHHOBZdccc0Nt9xxzwOPPPHMC6+88c4fPvjki29++PXPXwIMMsQwI4wyxjgTTDLFNDPMMsc8CyyyxDIrrLLGOhtsssU2O+yyxz4HHPKbI4454ZQzzrngkiuuueGWO+554JEnnnnhlTfe+cMHn3zxzQ//LPwBBhlimBFGGWOcCSaZYpoZZpljngUWWWKZFVZZY50NNtlimx122WOfAw75zRHHnHDKGedccMkV19xwyx33PPDIE8+88Mob7/zhg0+++OaHX//6JcAgQwwzwihjjDPBJFNMM8Msc8yzwCJLLLPCKmuss8EmW2yzwy577HPAIb854pgTTjnjnAsuueKaG265454HHnnimRdeeeOdP3zwyRff/PDr378EGGSIYUYYZYxxJphkimlmmGWOeRZYZIllVlhljXU22GSLbXbYZY99DjjkN0ccc8IpZ5xzwSVXXHPDLXfc88AjTzzzwitvvPOHDz754psffv3nlwCDDDHMCKOMMc4Ek0wxzQyzzDHPAossscwKq6yxzgabbLHNDrvssc8Bh/zmiGNOOOWMcy645IprbrjljnseeOSJZ1545Y13/vDBJ19888M/h/4AgwwxzAijjDHOBJNMMc0Ms8wxzwKLLLHMCqussc4Gm2yxzQ677LHPAYf85ohjTjjljHMuuOSKa2645Y57HnjkiWdeeOWNd/7wwSdffPPDr//+EmCQIYYZYZQxxplgkimmmWGWOeZZYJEllllhlTXW2WCTLbbZYZc99jngkN8cccwJp5xxzgWXXHHNDbfccc8DjzzxzAuvvPHOHz745Itvfvj1v18CDDLEMCOMMsY4E0wyxTQzzDLHPAssssQyK6yyxjobbLLFNjvsssc+BxzymyOOOeGUM8654JIrrrnhljvueeCRJ5554ZU33vnDB5988c0Pv/74JcAgQwwzwihjjDPBJFNMM8Msc8yzwCJLLLPCKmuss8EmW2yzwy577HPAIb854pgTTjnjnAsuueKaG265454HHnnimRdeeeOdP3zwyRff/PDPg///TKoWLw==AQAAAAAAAAAAgAAAAAAAAAAEAAAAAAAAEQAAAAAAAAA=eF7j5BwFo2AUjFQAABo4JAE= - </AppendedData> -</VTKFile> diff --git a/ogstools/studies/convergence/examples/steady_state_diffusion.prj b/ogstools/studies/convergence/examples/steady_state_diffusion.prj new file mode 100644 index 0000000000000000000000000000000000000000..c9f784e5a69de88215a3bf2154a22bbe7d687cbb --- /dev/null +++ b/ogstools/studies/convergence/examples/steady_state_diffusion.prj @@ -0,0 +1,132 @@ +<?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> diff --git a/ogstools/studies/convergence/examples/steady_state_diffusion.py b/ogstools/studies/convergence/examples/steady_state_diffusion.py new file mode 100644 index 0000000000000000000000000000000000000000..1144b180a7936481dba93f68e63c5311fe3ee776 --- /dev/null +++ b/ogstools/studies/convergence/examples/steady_state_diffusion.py @@ -0,0 +1,33 @@ +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 diff --git a/ogstools/studies/convergence/generate_report.py b/ogstools/studies/convergence/generate_report.py deleted file mode 100644 index 9ed2ecc40777addbdf4a47ff810a6a1ab271517b..0000000000000000000000000000000000000000 --- a/ogstools/studies/convergence/generate_report.py +++ /dev/null @@ -1,41 +0,0 @@ -import argparse -import tempfile -from pathlib import Path - -import jupytext -import papermill as pm - - -def execute_convergence_study( - output_path: str, mesh_paths: list[str], property_name: str, timestep: int -) -> None: - params = { - "mesh_paths": mesh_paths, - "property_name": property_name, - "timestep": timestep, - } - 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, - ) - 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("--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 - ) diff --git a/ogstools/studies/convergence/study.py b/ogstools/studies/convergence/study.py new file mode 100644 index 0000000000000000000000000000000000000000..710f286835f49d52616308ad0bb0c8b8a1fd8a3e --- /dev/null +++ b/ogstools/studies/convergence/study.py @@ -0,0 +1,56 @@ +import argparse +from pathlib import Path +from typing import Optional + +from ogstools.definitions import ROOT_DIR +from ogstools.workflow import jupytext_to_jupyter + + +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": [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, + } + 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, + ) + + +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() + + run_convergence_study( + args.output, + args.mesh_paths, + args.topology_path, + args.property_name, + args.timestep, + ) diff --git a/ogstools/studies/templates/convergence_study.py b/ogstools/studies/templates/convergence_study.py new file mode 100644 index 0000000000000000000000000000000000000000..3044b5cf661d9845f5c335f7a9c671fe45664b82 --- /dev/null +++ b/ogstools/studies/templates/convergence_study.py @@ -0,0 +1,128 @@ +# --- +# jupyter: +# kernelspec: +# display_name: .venv +# language: python +# name: python3 +# --- +# %% tags=["remove_cell"] +# mypy: ignore-errors +# +# This script is intended to be run via papermill with custom parameters. + +# %% tags=["parameters", "remove_cell"] +# These are placeholders and get replaced with injected parameters. +mesh_paths = [""] +topology_path = None +timestep: int = 0 +property_name: str = "" +refinement_ratio: float = None +reference_solution_path = None + +# %% tags=["remove_input"] +# Import required modules and customize plot settings. +from copy import deepcopy # 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 + +# %% tags=["remove_input"] +# Here, the meshes are read, a Property object is created from the property +# name, the topology is read and the Richardson extrapolation calculated. +# The 3 finest meshes of those provided will be used for the Richardson +# extrapolation. + +meshes: list[pv.DataSet] = [] +for mesh_path in mesh_paths: + if mesh_path.split(".")[-1] in ["xdmf", "pvd"]: + meshes += [meshlib.MeshSeries(mesh_path).read(timestep)] + else: + meshes += [pv.read(mesh_path)] +topology: pv.DataSet = pv.read(topology_path) +if property_name in meshes[0].point_data: + data_shape = meshes[0][property_name].shape +else: + data_shape = None +mesh_property = propertylib.presets.resolve_property(property_name, data_shape) +richardson = studies.convergence.richardson_extrapolation( + meshes, mesh_property, topology, refinement_ratio +) + +# %% [markdown] +# ## Grid convergence +# +# If the shown values are approximately 1, this means that the results are in +# asymptotic range of convergence. + +# %% tags=["remove_input"] +fig = meshplotlib.plot(richardson, "grid_convergence") + +# %% [markdown] +# ## Grid comparison +# +# Visualizing the requested mesh property on the 3 finest discretizations: + +# %% tags=["remove_input"] +fig = meshplotlib.plot(meshes[-3:], mesh_property) + +# %% [markdown] +# ## Richardson extrapolation. +# +# Visualizing the Richardson extrapolation of the requested mesh property on +# the chosen topology. 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. + +# %% tags=["remove_input"] +fig = meshplotlib.plot(richardson, mesh_property) + +if reference_solution_path is None: + reference_solution = richardson + diff = deepcopy(topology.sample(meshes[-1])) +else: + if reference_solution_path.split(".")[-1] in ["xdmf", "pvd"]: + reference_solution = meshlib.MeshSeries(reference_solution_path).read( + timestep + ) + else: + reference_solution = pv.read(reference_solution_path) + diff = deepcopy(richardson) +diff["difference"] = ( + reference_solution[mesh_property.data_name] - diff[mesh_property.data_name] +) +diff_unit = (mesh_property(1) - mesh_property(1)).units +diff_property = type(mesh_property)( + data_name="difference", + data_unit=diff_unit, + output_unit=diff_unit, + output_name=mesh_property.output_name + " difference", +) +fig = meshplotlib.plot(diff, diff_property) + +# %% [markdown] +# ## Convergence metrics + +# %% tags=["remove_input"] +metrics = studies.convergence.convergence_metrics( + meshes, richardson, mesh_property +) +metrics.style.format("{:,.5g}").hide() + +# %% [markdown] +# ## Relative errors + +# %% tags=["remove_input"] +meshplotlib.core.plt.rcdefaults() +fig = studies.convergence.plot_convergence_errors(metrics) + +# %% [markdown] +# ## Absolute values + +# %% tags=["remove_input"] +fig = studies.convergence.plot_convergence(metrics, mesh_property) diff --git a/ogstools/workflow/__init__.py b/ogstools/workflow/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..a21b980e4402f3f176eff4dda202e713093be225 --- /dev/null +++ b/ogstools/workflow/__init__.py @@ -0,0 +1,9 @@ +# Author: Florian Zill (Helmholtz Centre for Environmental Research GmbH - UFZ) +"""utility functions for workflows.""" + +from .jupyter_conversion import jupyter_to_html, jupytext_to_jupyter + +__all__ = [ + "jupytext_to_jupyter", + "jupyter_to_html", +] diff --git a/ogstools/workflow/jupyter_conversion.py b/ogstools/workflow/jupyter_conversion.py new file mode 100644 index 0000000000000000000000000000000000000000..b9c75a215cd25abdb1dcb437a784b6542fb8de08 --- /dev/null +++ b/ogstools/workflow/jupyter_conversion.py @@ -0,0 +1,54 @@ +import os +import sys +import tempfile +import warnings +from pathlib import Path + +from nbconvert import HTMLExporter +from nbconvert.preprocessors import TagRemovePreprocessor +from traitlets.config import Config + +# suppress output of version incompatibility between sphinx-gallery and jupytext +# it works anyways +with warnings.catch_warnings(): + warnings.filterwarnings("ignore") + sys.stdout = open(os.devnull, "w") # noqa: SIM115, PTH123 + import jupytext + + sys.stdout = sys.__stdout__ +import papermill as pm + + +def jupyter_to_html(input_path: Path, show_input: bool = False) -> str: + conf = Config() + conf.TagRemovePreprocessor.remove_cell_tags = ("remove_cell",) + hide_input_tags = ["injected-parameters"] + if not show_input: + hide_input_tags += ["remove_input"] + conf.TagRemovePreprocessor.remove_input_tags = hide_input_tags + exporter = HTMLExporter(config=conf) + exporter.register_preprocessor(TagRemovePreprocessor(config=conf), True) + exporter.exclude_input_prompt = True + exporter.exclude_output_prompt = True + report_html, _ = exporter.from_filename(str(input_path)) + return report_html + + +def jupytext_to_jupyter( + template_path: Path, + output_name: Path, + params: dict, + prepare_only: bool = False, + progress_bar: bool = False, +) -> None: + nb = jupytext.read(template_path) + with tempfile.NamedTemporaryFile(delete=False, suffix=".ipynb") as temp: + jupytext.write(nb, temp.name) + pm.execute_notebook( + input_path=temp.name, + output_path=output_name, + parameters=params, + prepare_only=prepare_only, + progress_bar=progress_bar, + ) + Path(temp.name).unlink() diff --git a/pyproject.toml b/pyproject.toml index e4b4fe7b3db4fc96c82dc39f51ceb8d72bc5cf98..2e9d974038e4e03cc2901dc20743659f332b898e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,6 +15,7 @@ dependencies = [ "jupytext>=1.15.2", "meshio==5.3.4", "matplotlib>=3.7.1", + "nbconvert>=7.9.2", "pyvista>=0.39.1", "rich>=13.4.2", "scipy>=1.10.1", diff --git a/tests/test_studies_convergence.py b/tests/test_studies_convergence.py index f8e7859e40ecc07a9b5aa6117e4a97ad1c1caa1c..6ca47f9e71da2717c0fa037044aadc964d867c90 100644 --- a/tests/test_studies_convergence.py +++ b/tests/test_studies_convergence.py @@ -1,50 +1,83 @@ """Unit tests for meshplotlib.""" import unittest +from pathlib import Path +from shutil import rmtree +from tempfile import mkdtemp import numpy as np +import pyvista as pv +from ogs6py import ogs -from ogstools.propertylib import Scalar -from ogstools.studies.convergence import ( - convergence_metrics, - grid_convergence, - log_fit, - plot_convergence, - plot_convergence_errors, - richardson_extrapolation, +from ogstools import meshlib, msh2vtu, propertylib +from ogstools.studies import convergence +from ogstools.studies.convergence.examples import ( + steady_state_diffusion_analytical_solution, ) -from ogstools.studies.convergence.examples import analytical_solution, meshes class ConvergenceTest(unittest.TestCase): """Test case for convergent meshes.""" def test_square_neumann_benchmark(self): - topology = meshes[-3] - mesh_property = Scalar("pressure", "m", "m") - conv = grid_convergence( - meshes, mesh_property, topology, refinement_ratio=2.0 + temp_dir = Path(mkdtemp()) + sim_results = [] + for i in range(3, 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) + model = ogs.OGS( + PROJECT_FILE=temp_dir / "default.prj", + INPUT_FILE=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) + sim_results += [pv.read(result_path)] + + topology = sim_results[-3] + spacing = convergence.add_grid_spacing(topology)["grid_spacing"] + np.testing.assert_array_less(0.0, spacing) + mesh_property = propertylib.Scalar("pressure", "m", "m") + conv = convergence.grid_convergence( + sim_results, mesh_property, topology, refinement_ratio=2.0 + ) + richardson = convergence.richardson_extrapolation( + sim_results, mesh_property, topology, refinement_ratio=2.0 ) - richardson = richardson_extrapolation(meshes, mesh_property, topology) np.testing.assert_allclose(conv["r"], richardson["r"], rtol=1e-10) - analytical = analytical_solution(topology) + analytical = steady_state_diffusion_analytical_solution(topology) np.testing.assert_allclose( richardson[mesh_property.data_name], analytical[mesh_property.data_name], rtol=2e-3, verbose=True, ) - metrics = convergence_metrics(meshes, richardson, mesh_property) + metrics = convergence.convergence_metrics( + sim_results, richardson, mesh_property + ) el_len = metrics["mean element length"].to_numpy()[:-1] re_max = metrics["rel. error (max)"].to_numpy()[:-1] ratio = re_max[:-1] / re_max[1:] - np.testing.assert_array_less(np.ones(ratio.shape) * 2.0, ratio) + np.testing.assert_array_less(2.0, ratio) - order, _ = log_fit(el_len, re_max) + order, _ = convergence.log_fit(el_len, re_max) self.assertGreater(order, 2.0) - _ = plot_convergence(metrics, mesh_property) - _ = plot_convergence_errors(metrics) + _ = convergence.plot_convergence(metrics, mesh_property) + _ = convergence.plot_convergence_errors(metrics) + + rmtree(temp_dir) if __name__ == "__main__":