From 80b125da4b8ebd286c2756bcb500d5ed883a2db2 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 5 Oct 2023 13:45:08 -0500 Subject: [PATCH 01/17] Use numpy linalg norm for convergence error --- polaris/ocean/convergence/spherical/analysis.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/polaris/ocean/convergence/spherical/analysis.py b/polaris/ocean/convergence/spherical/analysis.py index f2bdebc39..6bda1ee39 100644 --- a/polaris/ocean/convergence/spherical/analysis.py +++ b/polaris/ocean/convergence/spherical/analysis.py @@ -286,6 +286,7 @@ def compute_error(self, mesh_name, variable_name, zidx=None, error : float The error of the variable given by variable_name """ + norm_type = {'l2': None, 'inf': np.inf} ds_mesh = xr.open_dataset(f'{mesh_name}_mesh.nc') config = self.config section = config['spherical_convergence'] @@ -302,14 +303,14 @@ def compute_error(self, mesh_name, variable_name, zidx=None, if error_type == 'l2': area = area_for_field(ds_mesh, diff) - total_area = np.sum(area) - den_l2 = np.sum(field_exact**2 * area) / total_area - num_l2 = np.sum(diff**2 * area) / total_area - error = np.sqrt(num_l2) / np.sqrt(den_l2) - elif error_type == 'inf': - error = np.amax(diff) / np.amax(np.abs(field_exact)) - else: - raise ValueError(f'Unsupported error type {error_type}') + diff = diff * area + + error = np.linalg.norm(diff, ord=norm_type[error_type]) + + if error_type == 'l2': + field_exact = field_exact * area + den_l2 = np.linalg.norm(field_exact, ord=norm_type[error_type]) + error = np.divide(error, den_l2) return error From 2a4090871e12aee8babd7c1dddd08bc2fa1612a6 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 5 Oct 2023 15:39:21 -0500 Subject: [PATCH 02/17] Remove icos attribute from analysis step --- polaris/ocean/convergence/spherical/analysis.py | 11 +---------- polaris/ocean/tasks/cosine_bell/__init__.py | 2 +- polaris/ocean/tasks/cosine_bell/analysis.py | 8 +------- 3 files changed, 3 insertions(+), 18 deletions(-) diff --git a/polaris/ocean/convergence/spherical/analysis.py b/polaris/ocean/convergence/spherical/analysis.py index 6bda1ee39..38ebc900a 100644 --- a/polaris/ocean/convergence/spherical/analysis.py +++ b/polaris/ocean/convergence/spherical/analysis.py @@ -20,10 +20,6 @@ class SphericalConvergenceAnalysis(Step): resolutions : list of float The resolutions of the meshes that have been run - icosahedral : bool - Whether to use icosahedral, as opposed to less regular, JIGSAW - meshes - dependencies_dict : dict of dict of polaris.Steps The dependencies of this step must be given as separate keys in the dict: @@ -58,7 +54,7 @@ class SphericalConvergenceAnalysis(Step): The z-index to use for variables that have an nVertLevels dimension, which should be None for variables that don't """ - def __init__(self, component, resolutions, icosahedral, subdir, + def __init__(self, component, resolutions, subdir, dependencies, convergence_vars): """ Create the step @@ -71,10 +67,6 @@ def __init__(self, component, resolutions, icosahedral, subdir, resolutions : list of float The resolutions of the meshes that have been run - icosahedral : bool - Whether to use icosahedral, as opposed to less regular, JIGSAW - meshes - subdir : str The subdirectory that the step resides in @@ -114,7 +106,6 @@ def __init__(self, component, resolutions, icosahedral, subdir, """ super().__init__(component=component, name='analysis', subdir=subdir) self.resolutions = resolutions - self.icosahedral = icosahedral self.dependencies_dict = dependencies self.convergence_vars = convergence_vars diff --git a/polaris/ocean/tasks/cosine_bell/__init__.py b/polaris/ocean/tasks/cosine_bell/__init__.py index c3144c534..ce3182b18 100644 --- a/polaris/ocean/tasks/cosine_bell/__init__.py +++ b/polaris/ocean/tasks/cosine_bell/__init__.py @@ -197,7 +197,7 @@ def _setup_steps(self): step.dependencies_dict = analysis_dependencies else: step = Analysis(component=component, resolutions=resolutions, - icosahedral=icosahedral, subdir=subdir, + subdir=subdir, dependencies=analysis_dependencies) step.set_shared_config(config, link=config_filename) self.add_step(step, symlink=symlink) diff --git a/polaris/ocean/tasks/cosine_bell/analysis.py b/polaris/ocean/tasks/cosine_bell/analysis.py index 2d77ded8e..d205e65ba 100644 --- a/polaris/ocean/tasks/cosine_bell/analysis.py +++ b/polaris/ocean/tasks/cosine_bell/analysis.py @@ -11,8 +11,7 @@ class Analysis(SphericalConvergenceAnalysis): """ A step for analyzing the output from the cosine bell test case """ - def __init__(self, component, resolutions, icosahedral, subdir, - dependencies): + def __init__(self, component, resolutions, subdir, dependencies): """ Create the step @@ -24,10 +23,6 @@ def __init__(self, component, resolutions, icosahedral, subdir, resolutions : list of float The resolutions of the meshes that have been run - icosahedral : bool - Whether to use icosahedral, as opposed to less regular, JIGSAW - meshes - subdir : str The subdirectory that the step resides in @@ -39,7 +34,6 @@ def __init__(self, component, resolutions, icosahedral, subdir, 'zidx': 0}] super().__init__(component=component, subdir=subdir, resolutions=resolutions, - icosahedral=icosahedral, dependencies=dependencies, convergence_vars=convergence_vars) From c1e713bf8e330672c5bed73ccbb85b8dd1868a41 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Thu, 5 Oct 2023 16:53:53 -0500 Subject: [PATCH 03/17] Try using spherical convergence analysis step in IGW --- .../tasks/inertial_gravity_wave/__init__.py | 24 ++-- .../tasks/inertial_gravity_wave/analysis.py | 105 ++++++++---------- 2 files changed, 66 insertions(+), 63 deletions(-) diff --git a/polaris/ocean/tasks/inertial_gravity_wave/__init__.py b/polaris/ocean/tasks/inertial_gravity_wave/__init__.py index ab61bae4c..e20d9a025 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/__init__.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/__init__.py @@ -1,4 +1,6 @@ -from polaris import Task +from typing import Dict + +from polaris import Step, Task from polaris.ocean.tasks.inertial_gravity_wave.analysis import Analysis from polaris.ocean.tasks.inertial_gravity_wave.forward import Forward from polaris.ocean.tasks.inertial_gravity_wave.init import Init @@ -34,15 +36,23 @@ def __init__(self, component): super().__init__(component=component, name=name, subdir=subdir) self.resolutions = [200., 100., 50., 25.] - for res in self.resolutions: - self.add_step(Init(component=component, resolution=res, - taskdir=self.subdir)) - self.add_step(Forward(component=component, resolution=res, - taskdir=self.subdir)) + analysis_dependencies: Dict[str, Dict[float, Step]] = ( + dict(mesh=dict(), init=dict(), forward=dict())) + for resolution in self.resolutions: + init_step = Init(component=component, resolution=resolution, + taskdir=self.subdir) + self.add_step(init_step) + forward_step = Forward(component=component, resolution=resolution, + taskdir=self.subdir) + self.add_step(forward_step) + analysis_dependencies['mesh'][resolution] = init_step + analysis_dependencies['init'][resolution] = init_step + analysis_dependencies['forward'][resolution] = forward_step self.add_step(Analysis(component=component, resolutions=self.resolutions, - taskdir=self.subdir)) + subdir=f'{subdir}/analysis', + dependencies=analysis_dependencies)) self.add_step(Viz(component=component, resolutions=self.resolutions, taskdir=self.subdir), run_by_default=False) diff --git a/polaris/ocean/tasks/inertial_gravity_wave/analysis.py b/polaris/ocean/tasks/inertial_gravity_wave/analysis.py index b5b344d11..23c70bf0f 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/analysis.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/analysis.py @@ -1,18 +1,12 @@ -import datetime -import warnings - -import cmocean # noqa: F401 -import numpy as np import xarray as xr -from polaris import Step -from polaris.ocean.resolution import resolution_to_subdir +from polaris.ocean.convergence.spherical import SphericalConvergenceAnalysis from polaris.ocean.tasks.inertial_gravity_wave.exact_solution import ( ExactSolution, ) -class Analysis(Step): +class Analysis(SphericalConvergenceAnalysis): """ A step for analysing the output from the inertial gravity wave test case @@ -22,7 +16,7 @@ class Analysis(Step): resolutions : list of float The resolutions of the meshes that have been run """ - def __init__(self, component, resolutions, taskdir): + def __init__(self, component, resolutions, subdir, dependencies): """ Create the step @@ -34,54 +28,53 @@ def __init__(self, component, resolutions, taskdir): resolutions : list of float The resolutions of the meshes that have been run - taskdir : str - The subdirectory that the task belongs to - """ - super().__init__(component=component, name='analysis', indir=taskdir) - self.resolutions = resolutions - - for resolution in resolutions: - mesh_name = resolution_to_subdir(resolution) - self.add_input_file( - filename=f'init_{mesh_name}.nc', - target=f'../init/{mesh_name}/initial_state.nc') - self.add_input_file( - filename=f'output_{mesh_name}.nc', - target=f'../forward/{mesh_name}/output.nc') + subdir : str + The subdirectory that the step resides in - def run(self): + dependencies : dict of dict of polaris.Steps + The dependencies of this step """ - Run this step of the test case + convergence_vars = [{'name': 'ssh', + 'title': 'SSH', + 'zidx': None}] + super().__init__(component=component, subdir=subdir, + resolutions=resolutions, + dependencies=dependencies, + convergence_vars=convergence_vars) + + # TODO move conv thresh to convergence section + # section = config['inertial_gravity_wave'] + # conv_thresh = section.getfloat('conv_thresh') + + def exact_solution(self, mesh_name, field_name, time, zidx=None): """ - config = self.config - resolutions = self.resolutions - - section = config['inertial_gravity_wave'] - conv_thresh = section.getfloat('conv_thresh') - conv_max = section.getfloat('conv_max') - - rmse = [] - for i, res in enumerate(resolutions): - mesh_name = f'{res:g}km' - init = xr.open_dataset(f'init_{mesh_name}.nc') - ds = xr.open_dataset(f'output_{mesh_name}.nc') - exact = ExactSolution(init, config) - - t0 = datetime.datetime.strptime(ds.xtime.values[0].decode(), - '%Y-%m-%d_%H:%M:%S') - tf = datetime.datetime.strptime(ds.xtime.values[-1].decode(), - '%Y-%m-%d_%H:%M:%S') - t = (tf - t0).total_seconds() - ssh_model = ds.ssh.values[-1, :] - rmse.append(np.sqrt(np.mean((ssh_model - exact.ssh(t).values)**2))) + Get the exact solution - p = np.polyfit(np.log10(resolutions), np.log10(rmse), 1) - conv = p[0] - - if conv < conv_thresh: - raise ValueError(f'order of convergence ' - f' {conv} < min tolerence {conv_thresh}') - - if conv > conv_max: - warnings.warn(f'order of convergence ' - f'{conv} > max tolerence {conv_max}') + Parameters + ---------- + mesh_name : str + The mesh name which is the prefix for the initial condition file + + field_name : str + The name of the variable of which we evaluate convergence + For the default method, we use the same convergence rate for all + fields + + time : float + The time at which to evaluate the exact solution in seconds. + For the default method, we always use the initial state. + + zidx : int, optional + The z-index for the vertical level at which to evaluate the exact + solution + + Returns + ------- + solution : xarray.DataArray + The exact solution as derived from the initial condition + """ + init = xr.open_dataset(f'{mesh_name}_init.nc') + exact = ExactSolution(init, self.config) + if field_name != 'ssh': + raise ValueError(f'{field_name} is not currently supported') + return exact.ssh(time) From e6667334f881588f84e5b496591c655bad7e5cf2 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 6 Oct 2023 10:12:46 -0500 Subject: [PATCH 04/17] Fixup IGW viz step --- polaris/ocean/tasks/inertial_gravity_wave/viz.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/polaris/ocean/tasks/inertial_gravity_wave/viz.py b/polaris/ocean/tasks/inertial_gravity_wave/viz.py index 73dbf27a0..768ba101d 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/viz.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/viz.py @@ -73,9 +73,9 @@ def run(self): rmse = [] error_range = None for i, res in enumerate(resolutions): - ds_mesh = xr.open_dataset(f'mesh_{res}km.nc') - ds_init = xr.open_dataset(f'init_{res}km.nc') - ds = xr.open_dataset(f'output_{res}km.nc') + ds_mesh = xr.open_dataset(f'mesh_{res:g}km.nc') + ds_init = xr.open_dataset(f'init_{res:g}km.nc') + ds = xr.open_dataset(f'output_{res:g}km.nc') ds['maxLevelCell'] = ds_init.maxLevelCell exact = ExactSolution(ds_init, config) From 8e14679f49ca4bdddbd275272bb27c5065f3a81e Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 6 Oct 2023 10:34:20 -0500 Subject: [PATCH 05/17] Reorg convergence module --- polaris/ocean/convergence/__init__.py | 1 + .../convergence/{spherical => }/analysis.py | 12 ++-- polaris/ocean/convergence/convergence.cfg | 32 +++++++++++ .../ocean/convergence/spherical/__init__.py | 3 - .../ocean/convergence/spherical/forward.py | 2 +- .../ocean/convergence/spherical/spherical.cfg | 32 ----------- polaris/ocean/tasks/cosine_bell/__init__.py | 2 + polaris/ocean/tasks/cosine_bell/analysis.py | 4 +- .../ocean/tasks/cosine_bell/cosine_bell.cfg | 8 +-- .../tasks/inertial_gravity_wave/__init__.py | 2 + .../tasks/inertial_gravity_wave/analysis.py | 8 +-- .../tasks/inertial_gravity_wave/forward.py | 56 ++++++++++++++----- .../tasks/inertial_gravity_wave/forward.yaml | 7 ++- .../inertial_gravity_wave.cfg | 35 +++++++++++- 14 files changed, 129 insertions(+), 75 deletions(-) rename polaris/ocean/convergence/{spherical => }/analysis.py (97%) create mode 100644 polaris/ocean/convergence/convergence.cfg diff --git a/polaris/ocean/convergence/__init__.py b/polaris/ocean/convergence/__init__.py index e69de29bb..0602f4ef3 100644 --- a/polaris/ocean/convergence/__init__.py +++ b/polaris/ocean/convergence/__init__.py @@ -0,0 +1 @@ +from polaris.ocean.convergence.analysis import ConvergenceAnalysis diff --git a/polaris/ocean/convergence/spherical/analysis.py b/polaris/ocean/convergence/analysis.py similarity index 97% rename from polaris/ocean/convergence/spherical/analysis.py rename to polaris/ocean/convergence/analysis.py index 38ebc900a..0b82f1851 100644 --- a/polaris/ocean/convergence/spherical/analysis.py +++ b/polaris/ocean/convergence/analysis.py @@ -11,9 +11,9 @@ from polaris.viz import use_mplstyle -class SphericalConvergenceAnalysis(Step): +class ConvergenceAnalysis(Step): """ - A step for analyzing the output from the geostrophic convergence test + A step for analyzing the output from convergence tests Attributes ---------- @@ -54,8 +54,8 @@ class SphericalConvergenceAnalysis(Step): The z-index to use for variables that have an nVertLevels dimension, which should be None for variables that don't """ - def __init__(self, component, resolutions, subdir, - dependencies, convergence_vars): + def __init__(self, component, resolutions, subdir, dependencies, + convergence_vars): """ Create the step @@ -280,7 +280,7 @@ def compute_error(self, mesh_name, variable_name, zidx=None, norm_type = {'l2': None, 'inf': np.inf} ds_mesh = xr.open_dataset(f'{mesh_name}_mesh.nc') config = self.config - section = config['spherical_convergence'] + section = config['convergence'] eval_time = section.getfloat('convergence_eval_time') s_per_day = 86400.0 @@ -393,7 +393,7 @@ def convergence_parameters(self, field_name=None): The error norm to compute """ config = self.config - section = config['spherical_convergence'] + section = config['convergence'] conv_thresh = section.getfloat('convergence_thresh') error_type = section.get('error_type') return conv_thresh, error_type diff --git a/polaris/ocean/convergence/convergence.cfg b/polaris/ocean/convergence/convergence.cfg new file mode 100644 index 000000000..3b9e45ad1 --- /dev/null +++ b/polaris/ocean/convergence/convergence.cfg @@ -0,0 +1,32 @@ +[convergence] + +# Evaluation time for convergence analysis (in days) +convergence_eval_time = 1.0 + +# Convergence threshold below which a test fails +convergence_thresh = 1.0 + +# Type of error to compute +error_type = l2 + +# config options for convergence forward steps +[convergence_forward] + +# time integrator: {'split_explicit', 'RK4'} +time_integrator = RK4 + +# RK4 time step per resolution (s/km), since dt is proportional to resolution +rk4_dt_per_km = 3.0 + +# split time step per resolution (s/km), since dt is proportional to resolution +split_dt_per_km = 30.0 + +# the barotropic time step (s/km) for simulations using split time stepping, +# since btr_dt is proportional to resolution +btr_dt_per_km = 1.5 + +# Run duration in days +run_duration = ${convergence:convergence_eval_time} + +# Output interval in days +output_interval = ${run_duration} diff --git a/polaris/ocean/convergence/spherical/__init__.py b/polaris/ocean/convergence/spherical/__init__.py index 418252ded..702a1352f 100644 --- a/polaris/ocean/convergence/spherical/__init__.py +++ b/polaris/ocean/convergence/spherical/__init__.py @@ -1,6 +1,3 @@ -from polaris.ocean.convergence.spherical.analysis import ( - SphericalConvergenceAnalysis, -) from polaris.ocean.convergence.spherical.forward import ( SphericalConvergenceForward, ) diff --git a/polaris/ocean/convergence/spherical/forward.py b/polaris/ocean/convergence/spherical/forward.py index aa14f26d8..6bb7f9fae 100644 --- a/polaris/ocean/convergence/spherical/forward.py +++ b/polaris/ocean/convergence/spherical/forward.py @@ -111,7 +111,7 @@ def dynamic_model_config(self, at_setup): if not at_setup and vert_levels == 1: self.add_yaml_file('polaris.ocean.config', 'single_layer.yaml') - section = config['spherical_convergence_forward'] + section = config['convergence_forward'] time_integrator = section.get('time_integrator') diff --git a/polaris/ocean/convergence/spherical/spherical.cfg b/polaris/ocean/convergence/spherical/spherical.cfg index 630250da4..69ea2dd22 100644 --- a/polaris/ocean/convergence/spherical/spherical.cfg +++ b/polaris/ocean/convergence/spherical/spherical.cfg @@ -6,35 +6,3 @@ icos_resolutions = 60, 120, 240, 480 # a list of quasi-uniform mesh resolutions (km) to test qu_resolutions = 60, 90, 120, 150, 180, 210, 240 - -# Evaluation time for convergence analysis (in days) -convergence_eval_time = 1.0 - -# Convergence threshold below which a test fails -convergence_thresh = 1.0 - -# Type of error to compute -error_type = l2 - - -# config options for spherical convergence forward steps -[spherical_convergence_forward] - -# time integrator: {'split_explicit', 'RK4'} -time_integrator = RK4 - -# RK4 time step per resolution (s/km), since dt is proportional to resolution -rk4_dt_per_km = 3.0 - -# split time step per resolution (s/km), since dt is proportional to resolution -split_dt_per_km = 30.0 - -# the barotropic time step (s/km) for simulations using split time stepping, -# since btr_dt is proportional to resolution -btr_dt_per_km = 1.5 - -# Run duration in days -run_duration = ${spherical_convergence:convergence_eval_time} - -# Output interval in days -output_interval = ${run_duration} diff --git a/polaris/ocean/tasks/cosine_bell/__init__.py b/polaris/ocean/tasks/cosine_bell/__init__.py index ce3182b18..34cc3e052 100644 --- a/polaris/ocean/tasks/cosine_bell/__init__.py +++ b/polaris/ocean/tasks/cosine_bell/__init__.py @@ -21,6 +21,8 @@ def add_cosine_bell_tasks(component): filepath = f'spherical/{prefix}/cosine_bell/cosine_bell.cfg' config = PolarisConfigParser(filepath=filepath) + config.add_from_package('polaris.ocean.convergence', + 'convergence.cfg') config.add_from_package('polaris.ocean.convergence.spherical', 'spherical.cfg') config.add_from_package('polaris.ocean.tasks.cosine_bell', diff --git a/polaris/ocean/tasks/cosine_bell/analysis.py b/polaris/ocean/tasks/cosine_bell/analysis.py index d205e65ba..5200e2587 100644 --- a/polaris/ocean/tasks/cosine_bell/analysis.py +++ b/polaris/ocean/tasks/cosine_bell/analysis.py @@ -3,11 +3,11 @@ from mpas_tools.transects import lon_lat_to_cartesian from mpas_tools.vector import Vector -from polaris.ocean.convergence.spherical import SphericalConvergenceAnalysis +from polaris.ocean.convergence import ConvergenceAnalysis from polaris.ocean.tasks.cosine_bell.init import cosine_bell -class Analysis(SphericalConvergenceAnalysis): +class Analysis(ConvergenceAnalysis): """ A step for analyzing the output from the cosine bell test case """ diff --git a/polaris/ocean/tasks/cosine_bell/cosine_bell.cfg b/polaris/ocean/tasks/cosine_bell/cosine_bell.cfg index c2b8b5ff7..5486e33a5 100644 --- a/polaris/ocean/tasks/cosine_bell/cosine_bell.cfg +++ b/polaris/ocean/tasks/cosine_bell/cosine_bell.cfg @@ -20,8 +20,8 @@ partial_cell_type = None min_pc_fraction = 0.1 -# config options for spherical convergence tests -[spherical_convergence] +# config options for convergence tests +[convergence] # Evaluation time for convergence analysis (in days) convergence_eval_time = ${cosine_bell:vel_pd} @@ -33,8 +33,8 @@ convergence_thresh = ${cosine_bell:convergence_thresh} error_type = l2 -# config options for spherical convergence tests -[spherical_convergence_forward] +# config options for convergence forward steps +[convergence_forward] # time integrator: {'split_explicit', 'RK4'} time_integrator = RK4 diff --git a/polaris/ocean/tasks/inertial_gravity_wave/__init__.py b/polaris/ocean/tasks/inertial_gravity_wave/__init__.py index e20d9a025..628a1fae4 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/__init__.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/__init__.py @@ -57,6 +57,8 @@ def __init__(self, component): taskdir=self.subdir), run_by_default=False) + self.config.add_from_package('polaris.ocean.convergence', + 'convergence.cfg') self.config.add_from_package( 'polaris.ocean.tasks.inertial_gravity_wave', 'inertial_gravity_wave.cfg') diff --git a/polaris/ocean/tasks/inertial_gravity_wave/analysis.py b/polaris/ocean/tasks/inertial_gravity_wave/analysis.py index 23c70bf0f..c9383358e 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/analysis.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/analysis.py @@ -1,12 +1,12 @@ import xarray as xr -from polaris.ocean.convergence.spherical import SphericalConvergenceAnalysis +from polaris.ocean.convergence import ConvergenceAnalysis from polaris.ocean.tasks.inertial_gravity_wave.exact_solution import ( ExactSolution, ) -class Analysis(SphericalConvergenceAnalysis): +class Analysis(ConvergenceAnalysis): """ A step for analysing the output from the inertial gravity wave test case @@ -42,10 +42,6 @@ def __init__(self, component, resolutions, subdir, dependencies): dependencies=dependencies, convergence_vars=convergence_vars) - # TODO move conv thresh to convergence section - # section = config['inertial_gravity_wave'] - # conv_thresh = section.getfloat('conv_thresh') - def exact_solution(self, mesh_name, field_name, time, zidx=None): """ Get the exact solution diff --git a/polaris/ocean/tasks/inertial_gravity_wave/forward.py b/polaris/ocean/tasks/inertial_gravity_wave/forward.py index 9727b2f78..121f350f9 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/forward.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/forward.py @@ -1,9 +1,7 @@ -import time - import numpy as np from polaris.mesh.planar import compute_planar_hex_nx_ny -from polaris.ocean.model import OceanModelStep +from polaris.ocean.model import OceanModelStep, get_time_interval_string from polaris.ocean.resolution import resolution_to_subdir @@ -62,10 +60,8 @@ def __init__(self, component, resolution, taskdir, filename='output.nc', validate_vars=['layerThickness', 'normalVelocity']) - self.add_yaml_file('polaris.ocean.config', - 'single_layer.yaml') - self.add_yaml_file('polaris.ocean.tasks.inertial_gravity_wave', - 'forward.yaml') + self.package = 'polaris.ocean.tasks.inertial_gravity_wave' + self.yaml_filename = 'forward.yaml' def compute_cell_count(self): """ @@ -96,12 +92,42 @@ def dynamic_model_config(self, at_setup): """ super().dynamic_model_config(at_setup=at_setup) - # dt is proportional to resolution config = self.config - section = config['inertial_gravity_wave'] - dt_per_km = section.getfloat('dt_per_km') - dt = dt_per_km * self.resolution - # https://stackoverflow.com/a/1384565/7728169 - dt_str = time.strftime('%H:%M:%S', time.gmtime(dt)) - options = {'config_dt': dt_str} - self.add_model_config_options(options) + + vert_levels = config.getfloat('vertical_grid', 'vert_levels') + if not at_setup and vert_levels == 1: + self.add_yaml_file('polaris.ocean.config', 'single_layer.yaml') + + # dt is proportional to resolution + section = config['convergence_forward'] + + time_integrator = section.get('time_integrator') + + # dt is proportional to resolution: default 30 seconds per km + if time_integrator == 'RK4': + dt_per_km = section.getfloat('rk4_dt_per_km') + else: + dt_per_km = section.getfloat('split_dt_per_km') + dt_str = get_time_interval_string(seconds=dt_per_km * self.resolution) + + # btr_dt is also proportional to resolution: default 1.5 seconds per km + btr_dt_per_km = section.getfloat('btr_dt_per_km') + btr_dt_str = get_time_interval_string( + seconds=btr_dt_per_km * self.resolution) + + run_duration = section.getfloat('run_duration') + run_duration_str = get_time_interval_string(days=run_duration) + + output_interval = section.getfloat('output_interval') + output_interval_str = get_time_interval_string(days=output_interval) + + replacements = dict( + time_integrator=time_integrator, + dt=dt_str, + btr_dt=btr_dt_str, + run_duration=run_duration_str, + output_interval=output_interval_str, + ) + + self.add_yaml_file(self.package, self.yaml_filename, + template_replacements=replacements) diff --git a/polaris/ocean/tasks/inertial_gravity_wave/forward.yaml b/polaris/ocean/tasks/inertial_gravity_wave/forward.yaml index be242a139..bc3c2335e 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/forward.yaml +++ b/polaris/ocean/tasks/inertial_gravity_wave/forward.yaml @@ -1,8 +1,9 @@ omega: time_management: - config_run_duration: 10:00:00 + config_run_duration: {{ run_duration }} time_integration: - config_time_integrator: RK4 + config_dt: {{ dt }} + config_time_integrator: {{ time_integrator }} bottom_drag: config_bottom_drag_mode: implicit config_implicit_bottom_drag_type: constant @@ -21,7 +22,7 @@ omega: output: type: output filename_template: output.nc - output_interval: 0000_00:20:00 + output_interval: {{ output_interval }} clobber_mode: truncate contents: - xtime diff --git a/polaris/ocean/tasks/inertial_gravity_wave/inertial_gravity_wave.cfg b/polaris/ocean/tasks/inertial_gravity_wave/inertial_gravity_wave.cfg index f08ad30a0..84e8d4534 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/inertial_gravity_wave.cfg +++ b/polaris/ocean/tasks/inertial_gravity_wave/inertial_gravity_wave.cfg @@ -19,12 +19,12 @@ n_wavelengths_y = 2 # Convergence threshold below which the test fails conv_thresh = 1.8 -# Convergence rate above which a warning is issued -conv_max = 2.2 - # time step per resolution (s/km), since dt is proportional to resolution dt_per_km = 3.0 +# Run duration in days +run_duration = 0.4166667 + [vertical_grid] # The type of vertical grid @@ -45,6 +45,35 @@ partial_cell_type = None # The minimum fraction of a layer for partial cells min_pc_fraction = 0.1 +# config options for spherical convergence tests +[convergence] + +# Evaluation time for convergence analysis (in days) +convergence_eval_time = ${inertial_gravity_wave:run_duration} + +# Convergence threshold below which a test fails +convergence_thresh = ${inertial_gravity_wave:conv_thresh} + +# Type of error to compute +error_type = l2 + + +# config options for spherical convergence tests +[convergence_forward] + +# time integrator: {'split_explicit', 'RK4'} +time_integrator = RK4 + +# RK4 time step per resolution (s/km), since dt is proportional to resolution +rk4_dt_per_km = ${inertial_gravity_wave:dt_per_km} + +# Run duration in days +run_duration = ${inertial_gravity_wave:run_duration} + +# Output interval in days +output_interval = ${inertial_gravity_wave:run_duration} + + [ocean] # the number of cells per core to aim for goal_cells_per_core = 200 From 4222d42f3d39840c6c2542d0ff105cf97655fdcc Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 6 Oct 2023 11:17:52 -0500 Subject: [PATCH 06/17] Replace g with constants in IGW --- polaris/ocean/tasks/inertial_gravity_wave/exact_solution.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/polaris/ocean/tasks/inertial_gravity_wave/exact_solution.py b/polaris/ocean/tasks/inertial_gravity_wave/exact_solution.py index 5fe57c214..d117f1d05 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/exact_solution.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/exact_solution.py @@ -1,4 +1,5 @@ import numpy as np +from mpas_tools.cime.constants import constants class ExactSolution(): @@ -64,7 +65,7 @@ def __init__(self, ds, config): npx = section.getfloat('n_wavelengths_x') npy = section.getfloat('n_wavelengths_y') - self.g = 9.80616 + self.g = constants['SHR_CONST_G'] ly = np.sqrt(3.0) / 2.0 * lx self.kx = npx * 2.0 * np.pi / (lx * 1e3) self.ky = npy * 2.0 * np.pi / (ly * 1e3) From ae4daeeb1d30f9189773631e7c0aef90bf88f8f5 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 6 Oct 2023 11:19:30 -0500 Subject: [PATCH 07/17] Replace s_per_day with constants in cosine_bell --- polaris/ocean/tasks/cosine_bell/init.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/polaris/ocean/tasks/cosine_bell/init.py b/polaris/ocean/tasks/cosine_bell/init.py index 251918a04..5eb41aa80 100644 --- a/polaris/ocean/tasks/cosine_bell/init.py +++ b/polaris/ocean/tasks/cosine_bell/init.py @@ -1,5 +1,6 @@ import numpy as np import xarray as xr +from mpas_tools.cime.constants import constants from mpas_tools.io import write_netcdf from mpas_tools.transects import lon_lat_to_cartesian from mpas_tools.vector import Vector @@ -97,7 +98,7 @@ def run(self): ds['tracer3'] = ds.tracer1 # Initialize velocity - seconds_per_day = 86400.0 + seconds_per_day = constants['SHR_CONST_CDAY'] velocity = (2.0 * np.pi * np.cos(angleEdge) * sphere_radius * np.cos(latEdge) / (seconds_per_day * vel_pd)) velocity_array, _ = xr.broadcast(velocity, ds.refZMid) From 870ee9f0bd0c9e28d57d24764ab5bbd509a3db1b Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 6 Oct 2023 14:33:26 -0500 Subject: [PATCH 08/17] Reorg convergence forward step --- polaris/ocean/convergence/__init__.py | 1 + polaris/ocean/convergence/forward.py | 144 ++++++++++++++++++ .../ocean/convergence/spherical/forward.py | 81 +--------- .../tasks/inertial_gravity_wave/__init__.py | 9 +- .../tasks/inertial_gravity_wave/forward.py | 96 ++---------- .../tasks/inertial_gravity_wave/forward.yaml | 4 +- .../ocean/tasks/inertial_gravity_wave/init.py | 2 + 7 files changed, 172 insertions(+), 165 deletions(-) create mode 100644 polaris/ocean/convergence/forward.py diff --git a/polaris/ocean/convergence/__init__.py b/polaris/ocean/convergence/__init__.py index 0602f4ef3..7fb71d53f 100644 --- a/polaris/ocean/convergence/__init__.py +++ b/polaris/ocean/convergence/__init__.py @@ -1 +1,2 @@ from polaris.ocean.convergence.analysis import ConvergenceAnalysis +from polaris.ocean.convergence.forward import ConvergenceForward diff --git a/polaris/ocean/convergence/forward.py b/polaris/ocean/convergence/forward.py new file mode 100644 index 000000000..edd866ece --- /dev/null +++ b/polaris/ocean/convergence/forward.py @@ -0,0 +1,144 @@ +from polaris.ocean.model import OceanModelStep, get_time_interval_string + + +class ConvergenceForward(OceanModelStep): + """ + A step for performing forward ocean component runs as part of a spherical + convergence test + + Attributes + ---------- + resolution : float + The resolution of the (uniform) mesh in km + + package : Package + The package name or module object that contains the YAML file + + yaml_filename : str + A YAML file that is a Jinja2 template with model config options + + """ + + def __init__(self, component, name, subdir, resolution, base_mesh, init, + package, yaml_filename='forward.yaml', options=None, + output_filename='output.nc', validate_vars=None): + """ + Create a new step + + Parameters + ---------- + component : polaris.Component + The component the step belongs to + + name : str + The name of the step + + subdir : str + The subdirectory for the step + + resolution : float + The resolution of the (uniform) mesh in km + + package : Package + The package name or module object that contains the YAML file + + yaml_filename : str, optional + A YAML file that is a Jinja2 template with model config options + + options : dict, optional + A dictionary of options and value to replace model config options + with new values + + output_filename : str, optional + The output file that will be written out at the end of the forward + run + + validate_vars : list of str, optional + A list of variables to validate against a baseline if requested + """ + super().__init__(component=component, name=name, subdir=subdir, + openmp_threads=1) + + self.resolution = resolution + self.package = package + self.yaml_filename = yaml_filename + + # make sure output is double precision + self.add_yaml_file('polaris.ocean.config', 'output.yaml') + + if options is not None: + self.add_model_config_options(options=options) + + self.add_input_file( + filename='init.nc', + work_dir_target=f'{init.path}/initial_state.nc') + self.add_input_file( + filename='graph.info', + work_dir_target=f'{base_mesh.path}/graph.info') + + self.add_output_file(filename=output_filename, + validate_vars=validate_vars) + + def compute_cell_count(self): + """ + Compute the approximate number of cells in the mesh, used to constrain + resources + + Returns + ------- + cell_count : int or None + The approximate number of cells in the mesh + """ + raise ValueError('compute_cell_count method must be overridden by ' + 'spherical or planar method') + + def dynamic_model_config(self, at_setup): + """ + Set the model time step from config options at setup and runtime + + Parameters + ---------- + at_setup : bool + Whether this method is being run during setup of the step, as + opposed to at runtime + """ + super().dynamic_model_config(at_setup=at_setup) + + config = self.config + + vert_levels = config.getfloat('vertical_grid', 'vert_levels') + if not at_setup and vert_levels == 1: + self.add_yaml_file('polaris.ocean.config', 'single_layer.yaml') + + section = config['convergence_forward'] + + time_integrator = section.get('time_integrator') + + # dt is proportional to resolution: default 30 seconds per km + if time_integrator == 'RK4': + dt_per_km = section.getfloat('rk4_dt_per_km') + else: + dt_per_km = section.getfloat('split_dt_per_km') + dt_str = get_time_interval_string(seconds=dt_per_km * self.resolution) + + # btr_dt is also proportional to resolution: default 1.5 seconds per km + btr_dt_per_km = section.getfloat('btr_dt_per_km') + btr_dt_str = get_time_interval_string( + seconds=btr_dt_per_km * self.resolution) + + run_duration = section.getfloat('run_duration') + run_duration_str = get_time_interval_string(days=run_duration) + + output_interval = section.getfloat('output_interval') + output_interval_str = get_time_interval_string(days=output_interval) + + replacements = dict( + time_integrator=time_integrator, + dt=dt_str, + btr_dt=btr_dt_str, + run_duration=run_duration_str, + output_interval=output_interval_str, + ) + + self.add_yaml_file(self.package, self.yaml_filename, + template_replacements=replacements) diff --git a/polaris/ocean/convergence/spherical/forward.py b/polaris/ocean/convergence/spherical/forward.py index 6bb7f9fae..b7c55cbb0 100644 --- a/polaris/ocean/convergence/spherical/forward.py +++ b/polaris/ocean/convergence/spherical/forward.py @@ -1,7 +1,7 @@ -from polaris.ocean.model import OceanModelStep, get_time_interval_string +from polaris.ocean.convergence import ConvergenceForward -class SphericalConvergenceForward(OceanModelStep): +class SphericalConvergenceForward(ConvergenceForward): """ A step for performing forward ocean component runs as part of a spherical convergence test @@ -57,27 +57,11 @@ def __init__(self, component, name, subdir, resolution, base_mesh, init, A list of variables to validate against a baseline if requested """ super().__init__(component=component, name=name, subdir=subdir, - openmp_threads=1) - - self.resolution = resolution - self.package = package - self.yaml_filename = yaml_filename - - # make sure output is double precision - self.add_yaml_file('polaris.ocean.config', 'output.yaml') - - if options is not None: - self.add_model_config_options(options=options) - - self.add_input_file( - filename='init.nc', - work_dir_target=f'{init.path}/initial_state.nc') - self.add_input_file( - filename='graph.info', - work_dir_target=f'{base_mesh.path}/graph.info') - - self.add_output_file(filename=output_filename, - validate_vars=validate_vars) + resolution=resolution, base_mesh=base_mesh, + init=init, package=package, + yaml_filename=yaml_filename, + output_filename=output_filename, + validate_vars=validate_vars) def compute_cell_count(self): """ @@ -92,54 +76,3 @@ def compute_cell_count(self): # use a heuristic based on QU30 (65275 cells) and QU240 (10383 cells) cell_count = 6e8 / self.resolution**2 return cell_count - - def dynamic_model_config(self, at_setup): - """ - Set the model time step from config options at setup and runtime - - Parameters - ---------- - at_setup : bool - Whether this method is being run during setup of the step, as - opposed to at runtime - """ - super().dynamic_model_config(at_setup=at_setup) - - config = self.config - - vert_levels = config.getfloat('vertical_grid', 'vert_levels') - if not at_setup and vert_levels == 1: - self.add_yaml_file('polaris.ocean.config', 'single_layer.yaml') - - section = config['convergence_forward'] - - time_integrator = section.get('time_integrator') - - # dt is proportional to resolution: default 30 seconds per km - if time_integrator == 'RK4': - dt_per_km = section.getfloat('rk4_dt_per_km') - else: - dt_per_km = section.getfloat('split_dt_per_km') - dt_str = get_time_interval_string(seconds=dt_per_km * self.resolution) - - # btr_dt is also proportional to resolution: default 1.5 seconds per km - btr_dt_per_km = section.getfloat('btr_dt_per_km') - btr_dt_str = get_time_interval_string( - seconds=btr_dt_per_km * self.resolution) - - run_duration = section.getfloat('run_duration') - run_duration_str = get_time_interval_string(days=run_duration) - - output_interval = section.getfloat('output_interval') - output_interval_str = get_time_interval_string(days=output_interval) - - replacements = dict( - time_integrator=time_integrator, - dt=dt_str, - btr_dt=btr_dt_str, - run_duration=run_duration_str, - output_interval=output_interval_str, - ) - - self.add_yaml_file(self.package, self.yaml_filename, - template_replacements=replacements) diff --git a/polaris/ocean/tasks/inertial_gravity_wave/__init__.py b/polaris/ocean/tasks/inertial_gravity_wave/__init__.py index 628a1fae4..a3023a74e 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/__init__.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/__init__.py @@ -1,6 +1,7 @@ from typing import Dict from polaris import Step, Task +from polaris.ocean.resolution import resolution_to_subdir from polaris.ocean.tasks.inertial_gravity_wave.analysis import Analysis from polaris.ocean.tasks.inertial_gravity_wave.forward import Forward from polaris.ocean.tasks.inertial_gravity_wave.init import Init @@ -39,11 +40,15 @@ def __init__(self, component): analysis_dependencies: Dict[str, Dict[float, Step]] = ( dict(mesh=dict(), init=dict(), forward=dict())) for resolution in self.resolutions: + mesh_name = resolution_to_subdir(resolution) init_step = Init(component=component, resolution=resolution, taskdir=self.subdir) self.add_step(init_step) - forward_step = Forward(component=component, resolution=resolution, - taskdir=self.subdir) + forward_step = Forward(component=component, + name=f'forward_{mesh_name}', + resolution=resolution, + subdir=f'{self.subdir}/forward/{mesh_name}', + init=init_step) self.add_step(forward_step) analysis_dependencies['mesh'][resolution] = init_step analysis_dependencies['init'][resolution] = init_step diff --git a/polaris/ocean/tasks/inertial_gravity_wave/forward.py b/polaris/ocean/tasks/inertial_gravity_wave/forward.py index 121f350f9..480328db0 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/forward.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/forward.py @@ -1,11 +1,10 @@ import numpy as np from polaris.mesh.planar import compute_planar_hex_nx_ny -from polaris.ocean.model import OceanModelStep, get_time_interval_string -from polaris.ocean.resolution import resolution_to_subdir +from polaris.ocean.convergence import ConvergenceForward -class Forward(OceanModelStep): +class Forward(ConvergenceForward): """ A step for performing forward ocean component runs as part of inertial gravity wave test cases. @@ -15,8 +14,7 @@ class Forward(OceanModelStep): resolution : float The resolution of the test case in km """ - def __init__(self, component, resolution, taskdir, - ntasks=None, min_tasks=None, openmp_threads=1): + def __init__(self, component, name, resolution, subdir, init): """ Create a new test case @@ -30,38 +28,14 @@ def __init__(self, component, resolution, taskdir, taskdir : str The subdirectory that the task belongs to - - ntasks : int, optional - the number of tasks the step would ideally use. If fewer tasks - are available on the system, the step will run on all available - tasks as long as this is not below ``min_tasks`` - - min_tasks : int, optional - the number of tasks the step requires. If the system has fewer - than this number of tasks, the step will fail - - openmp_threads : int, optional - the number of OpenMP threads the step will use """ - self.resolution = resolution - mesh_name = resolution_to_subdir(resolution) super().__init__(component=component, - name=f'forward_{mesh_name}', - subdir=f'{taskdir}/forward/{mesh_name}', - ntasks=ntasks, min_tasks=min_tasks, - openmp_threads=openmp_threads) - - self.add_input_file(filename='initial_state.nc', - target=f'../../init/{mesh_name}/initial_state.nc') - self.add_input_file(filename='graph.info', - target=f'../../init/{mesh_name}/culled_graph.info') - - self.add_output_file( - filename='output.nc', - validate_vars=['layerThickness', 'normalVelocity']) - - self.package = 'polaris.ocean.tasks.inertial_gravity_wave' - self.yaml_filename = 'forward.yaml' + name=name, subdir=subdir, + resolution=resolution, base_mesh=init, init=init, + package='polaris.ocean.tasks.inertial_gravity_wave', + yaml_filename='forward.yaml', + output_filename='output.nc', + validate_vars=['layerThickness', 'normalVelocity']) def compute_cell_count(self): """ @@ -79,55 +53,3 @@ def compute_cell_count(self): nx, ny = compute_planar_hex_nx_ny(lx, ly, self.resolution) cell_count = nx * ny return cell_count - - def dynamic_model_config(self, at_setup): - """ - Set the model time step from config options at setup and runtime - - Parameters - ---------- - at_setup : bool - Whether this method is being run during setup of the step, as - opposed to at runtime - """ - super().dynamic_model_config(at_setup=at_setup) - - config = self.config - - vert_levels = config.getfloat('vertical_grid', 'vert_levels') - if not at_setup and vert_levels == 1: - self.add_yaml_file('polaris.ocean.config', 'single_layer.yaml') - - # dt is proportional to resolution - section = config['convergence_forward'] - - time_integrator = section.get('time_integrator') - - # dt is proportional to resolution: default 30 seconds per km - if time_integrator == 'RK4': - dt_per_km = section.getfloat('rk4_dt_per_km') - else: - dt_per_km = section.getfloat('split_dt_per_km') - dt_str = get_time_interval_string(seconds=dt_per_km * self.resolution) - - # btr_dt is also proportional to resolution: default 1.5 seconds per km - btr_dt_per_km = section.getfloat('btr_dt_per_km') - btr_dt_str = get_time_interval_string( - seconds=btr_dt_per_km * self.resolution) - - run_duration = section.getfloat('run_duration') - run_duration_str = get_time_interval_string(days=run_duration) - - output_interval = section.getfloat('output_interval') - output_interval_str = get_time_interval_string(days=output_interval) - - replacements = dict( - time_integrator=time_integrator, - dt=dt_str, - btr_dt=btr_dt_str, - run_duration=run_duration_str, - output_interval=output_interval_str, - ) - - self.add_yaml_file(self.package, self.yaml_filename, - template_replacements=replacements) diff --git a/polaris/ocean/tasks/inertial_gravity_wave/forward.yaml b/polaris/ocean/tasks/inertial_gravity_wave/forward.yaml index bc3c2335e..13dae39fa 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/forward.yaml +++ b/polaris/ocean/tasks/inertial_gravity_wave/forward.yaml @@ -15,9 +15,9 @@ omega: config_disable_vel_hmix: true streams: mesh: - filename_template: initial_state.nc + filename_template: init.nc input: - filename_template: initial_state.nc + filename_template: init.nc restart: {} output: type: output diff --git a/polaris/ocean/tasks/inertial_gravity_wave/init.py b/polaris/ocean/tasks/inertial_gravity_wave/init.py index 4935e912a..544e5176b 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/init.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/init.py @@ -5,6 +5,7 @@ from mpas_tools.planar_hex import make_planar_hex_mesh from polaris import Step +from polaris.io import symlink from polaris.mesh.planar import compute_planar_hex_nx_ny from polaris.ocean.resolution import resolution_to_subdir from polaris.ocean.tasks.inertial_gravity_wave.exact_solution import ( @@ -72,6 +73,7 @@ def run(self): ds_mesh = cull(ds_mesh, logger=logger) ds_mesh = convert(ds_mesh, graphInfoFileName='culled_graph.info', logger=logger) + symlink('culled_graph.info', 'graph.info') write_netcdf(ds_mesh, 'culled_mesh.nc') bottom_depth = config.getfloat('vertical_grid', 'bottom_depth') From 1d766ca76688c2690e3fa2ed94beba6955e041e4 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 6 Oct 2023 15:41:20 -0500 Subject: [PATCH 09/17] Use shared convergence for manufactured solution --- .../tasks/manufactured_solution/__init__.py | 30 ++++-- .../tasks/manufactured_solution/analysis.py | 101 ++++++++---------- .../tasks/manufactured_solution/forward.py | 60 ++--------- .../tasks/manufactured_solution/forward.yaml | 11 +- .../ocean/tasks/manufactured_solution/init.py | 5 + .../manufactured_solution.cfg | 32 +++++- 6 files changed, 120 insertions(+), 119 deletions(-) diff --git a/polaris/ocean/tasks/manufactured_solution/__init__.py b/polaris/ocean/tasks/manufactured_solution/__init__.py index 056e89047..2c95b634e 100644 --- a/polaris/ocean/tasks/manufactured_solution/__init__.py +++ b/polaris/ocean/tasks/manufactured_solution/__init__.py @@ -1,4 +1,7 @@ -from polaris import Task +from typing import Dict + +from polaris import Step, Task +from polaris.ocean.resolution import resolution_to_subdir from polaris.ocean.tasks.manufactured_solution.analysis import Analysis from polaris.ocean.tasks.manufactured_solution.forward import Forward from polaris.ocean.tasks.manufactured_solution.init import Init @@ -39,19 +42,32 @@ def __init__(self, component): super().__init__(component=component, name=name, subdir=subdir) self.resolutions = [200., 100., 50., 25.] - for res in self.resolutions: - self.add_step(Init(component=component, resolution=res, - taskdir=self.subdir)) - self.add_step(Forward(component=component, resolution=res, - taskdir=self.subdir)) + analysis_dependencies: Dict[str, Dict[float, Step]] = ( + dict(mesh=dict(), init=dict(), forward=dict())) + for resolution in self.resolutions: + mesh_name = resolution_to_subdir(resolution) + init_step = Init(component=component, resolution=resolution, + taskdir=self.subdir) + self.add_step(init_step) + forward_step = Forward(component=component, resolution=resolution, + name=f'forward_{mesh_name}', + subdir=f'{self.subdir}/forward/{mesh_name}', + init=init_step) + self.add_step(forward_step) + analysis_dependencies['mesh'][resolution] = init_step + analysis_dependencies['init'][resolution] = init_step + analysis_dependencies['forward'][resolution] = forward_step self.add_step(Analysis(component=component, resolutions=self.resolutions, - taskdir=self.subdir)) + subdir=f'{self.subdir}/analysis', + dependencies=analysis_dependencies)) self.add_step(Viz(component=component, resolutions=self.resolutions, taskdir=self.subdir), run_by_default=False) + self.config.add_from_package('polaris.ocean.convergence', + 'convergence.cfg') self.config.add_from_package( 'polaris.ocean.tasks.manufactured_solution', 'manufactured_solution.cfg') diff --git a/polaris/ocean/tasks/manufactured_solution/analysis.py b/polaris/ocean/tasks/manufactured_solution/analysis.py index a1748049e..25a17e717 100644 --- a/polaris/ocean/tasks/manufactured_solution/analysis.py +++ b/polaris/ocean/tasks/manufactured_solution/analysis.py @@ -1,18 +1,12 @@ -import datetime -import warnings - -import cmocean # noqa: F401 -import numpy as np import xarray as xr -from polaris import Step -from polaris.ocean.resolution import resolution_to_subdir +from polaris.ocean.convergence import ConvergenceAnalysis from polaris.ocean.tasks.manufactured_solution.exact_solution import ( ExactSolution, ) -class Analysis(Step): +class Analysis(ConvergenceAnalysis): """ A step for analysing the output from the manufactured solution test case @@ -22,7 +16,7 @@ class Analysis(Step): resolutions : list of float The resolutions of the meshes that have been run """ - def __init__(self, component, resolutions, taskdir): + def __init__(self, component, resolutions, subdir, dependencies): """ Create the step @@ -34,54 +28,49 @@ def __init__(self, component, resolutions, taskdir): resolutions : list of float The resolutions of the meshes that have been run - taskdir : str - The subdirectory that the task belongs to - """ - super().__init__(component=component, name='analysis', indir=taskdir) - self.resolutions = resolutions - - for resolution in resolutions: - mesh_name = resolution_to_subdir(resolution) - self.add_input_file( - filename=f'init_{mesh_name}.nc', - target=f'../init/{mesh_name}/initial_state.nc') - self.add_input_file( - filename=f'output_{mesh_name}.nc', - target=f'../forward/{mesh_name}/output.nc') + subdir : str + The subdirectory that the step resides in - def run(self): + dependencies : dict of dict of polaris.Steps + The dependencies of this step """ - Run this step of the test case + convergence_vars = [{'name': 'ssh', + 'title': 'SSH', + 'zidx': None}] + super().__init__(component=component, subdir=subdir, + resolutions=resolutions, + dependencies=dependencies, + convergence_vars=convergence_vars) + + def exact_solution(self, mesh_name, field_name, time, zidx=None): """ - config = self.config - resolutions = self.resolutions - - section = config['manufactured_solution'] - conv_thresh = section.getfloat('conv_thresh') - conv_max = section.getfloat('conv_max') - - rmse = [] - for i, res in enumerate(resolutions): - mesh_name = f'{res:g}km' - init = xr.open_dataset(f'init_{mesh_name}.nc') - ds = xr.open_dataset(f'output_{mesh_name}.nc') - exact = ExactSolution(config, init) - - t0 = datetime.datetime.strptime(ds.xtime.values[0].decode(), - '%Y-%m-%d_%H:%M:%S') - tf = datetime.datetime.strptime(ds.xtime.values[-1].decode(), - '%Y-%m-%d_%H:%M:%S') - t = (tf - t0).total_seconds() - ssh_model = ds.ssh.values[-1, :] - rmse.append(np.sqrt(np.mean((ssh_model - exact.ssh(t).values)**2))) + Get the exact solution - p = np.polyfit(np.log10(resolutions), np.log10(rmse), 1) - conv = p[0] - - if conv < conv_thresh: - raise ValueError(f'order of convergence ' - f' {conv} < min tolerence {conv_thresh}') - - if conv > conv_max: - warnings.warn(f'order of convergence ' - f'{conv} > max tolerence {conv_max}') + Parameters + ---------- + mesh_name : str + The mesh name which is the prefix for the initial condition file + + field_name : str + The name of the variable of which we evaluate convergence + For the default method, we use the same convergence rate for all + fields + + time : float + The time at which to evaluate the exact solution in seconds. + For the default method, we always use the initial state. + + zidx : int, optional + The z-index for the vertical level at which to evaluate the exact + solution + + Returns + ------- + solution : xarray.DataArray + The exact solution as derived from the initial condition + """ + init = xr.open_dataset(f'{mesh_name}_init.nc') + exact = ExactSolution(self.config, init) + if field_name != 'ssh': + raise ValueError(f'{field_name} is not currently supported') + return exact.ssh(time) diff --git a/polaris/ocean/tasks/manufactured_solution/forward.py b/polaris/ocean/tasks/manufactured_solution/forward.py index d8f56c92b..198d38c1b 100644 --- a/polaris/ocean/tasks/manufactured_solution/forward.py +++ b/polaris/ocean/tasks/manufactured_solution/forward.py @@ -1,16 +1,13 @@ -import time - import numpy as np from polaris.mesh.planar import compute_planar_hex_nx_ny -from polaris.ocean.model import OceanModelStep -from polaris.ocean.resolution import resolution_to_subdir +from polaris.ocean.convergence import ConvergenceForward from polaris.ocean.tasks.manufactured_solution.exact_solution import ( ExactSolution, ) -class Forward(OceanModelStep): +class Forward(ConvergenceForward): """ A step for performing forward ocean component runs as part of manufactured solution test cases. @@ -20,8 +17,7 @@ class Forward(OceanModelStep): resolution : float The resolution of the test case in km """ - def __init__(self, component, resolution, taskdir, - ntasks=None, min_tasks=None, openmp_threads=1): + def __init__(self, component, name, resolution, subdir, init): """ Create a new test case @@ -35,40 +31,14 @@ def __init__(self, component, resolution, taskdir, taskdir : str The subdirectory that the task belongs to - - ntasks : int, optional - the number of tasks the step would ideally use. If fewer tasks - are available on the system, the step will run on all available - tasks as long as this is not below ``min_tasks`` - - min_tasks : int, optional - the number of tasks the step requires. If the system has fewer - than this number of tasks, the step will fail - - openmp_threads : int, optional - the number of OpenMP threads the step will use """ - self.resolution = resolution - mesh_name = resolution_to_subdir(resolution) super().__init__(component=component, - name=f'forward_{mesh_name}', - subdir=f'{taskdir}/forward/{mesh_name}', - ntasks=ntasks, min_tasks=min_tasks, - openmp_threads=openmp_threads) - - self.add_input_file(filename='initial_state.nc', - target=f'../../init/{mesh_name}/initial_state.nc') - self.add_input_file(filename='graph.info', - target=f'../../init/{mesh_name}/culled_graph.info') - - self.add_output_file( - filename='output.nc', - validate_vars=['layerThickness', 'normalVelocity']) - - self.add_yaml_file('polaris.ocean.config', - 'single_layer.yaml') - self.add_yaml_file('polaris.ocean.tasks.manufactured_solution', - 'forward.yaml') + name=name, subdir=subdir, + resolution=resolution, base_mesh=init, init=init, + package='polaris.ocean.tasks.manufactured_solution', + yaml_filename='forward.yaml', + output_filename='output.nc', + validate_vars=['layerThickness', 'normalVelocity']) def compute_cell_count(self): """ @@ -100,16 +70,8 @@ def dynamic_model_config(self, at_setup): """ super().dynamic_model_config(at_setup=at_setup) - # dt is proportional to resolution - config = self.config - section = config['manufactured_solution'] - dt_per_km = section.getfloat('dt_per_km') - dt = dt_per_km * self.resolution - # https://stackoverflow.com/a/1384565/7728169 - dt_str = time.strftime('%H:%M:%S', time.gmtime(dt)) - exact_solution = ExactSolution(config) - options = {'config_dt': dt_str, - 'config_manufactured_solution_amplitude': + exact_solution = ExactSolution(self.config) + options = {'config_manufactured_solution_amplitude': exact_solution.eta0, 'config_manufactured_solution_wavelength_x': exact_solution.lambda_x, diff --git a/polaris/ocean/tasks/manufactured_solution/forward.yaml b/polaris/ocean/tasks/manufactured_solution/forward.yaml index 6b3502fce..aa40ed4b1 100644 --- a/polaris/ocean/tasks/manufactured_solution/forward.yaml +++ b/polaris/ocean/tasks/manufactured_solution/forward.yaml @@ -1,8 +1,9 @@ omega: time_management: - config_run_duration: 10:00:00 + config_run_duration: {{ run_duration }} time_integration: - config_time_integrator: RK4 + config_dt: {{ dt }} + config_time_integrator: {{ time_integrator }} bottom_drag: config_bottom_drag_mode: implicit config_implicit_bottom_drag_type: constant @@ -13,14 +14,14 @@ omega: config_disable_vel_hmix: true streams: mesh: - filename_template: initial_state.nc + filename_template: init.nc input: - filename_template: initial_state.nc + filename_template: init.nc restart: {} output: type: output filename_template: output.nc - output_interval: 0000_00:20:00 + output_interval: {{ output_interval }} clobber_mode: truncate contents: - xtime diff --git a/polaris/ocean/tasks/manufactured_solution/init.py b/polaris/ocean/tasks/manufactured_solution/init.py index 73c997d48..7942b10b9 100644 --- a/polaris/ocean/tasks/manufactured_solution/init.py +++ b/polaris/ocean/tasks/manufactured_solution/init.py @@ -5,6 +5,7 @@ from mpas_tools.planar_hex import make_planar_hex_mesh from polaris import Step +from polaris.io import symlink from polaris.mesh.planar import compute_planar_hex_nx_ny from polaris.ocean.resolution import resolution_to_subdir from polaris.ocean.tasks.manufactured_solution.exact_solution import ( @@ -43,6 +44,9 @@ def __init__(self, component, resolution, taskdir): name=f'init_{mesh_name}', subdir=f'{taskdir}/init/{mesh_name}') self.resolution = resolution + for filename in ['culled_mesh.nc', 'initial_state.nc', + 'culled_graph.info']: + self.add_output_file(filename=filename) def run(self): """ @@ -69,6 +73,7 @@ def run(self): ds_mesh = cull(ds_mesh, logger=logger) ds_mesh = convert(ds_mesh, graphInfoFileName='culled_graph.info', logger=logger) + symlink('culled_graph.info', 'graph.info') write_netcdf(ds_mesh, 'culled_mesh.nc') bottom_depth = config.getfloat('vertical_grid', 'bottom_depth') diff --git a/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg b/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg index 14d42ca80..f241fcb8e 100644 --- a/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg +++ b/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg @@ -30,8 +30,8 @@ dt_per_km = 3.0 # Convergence threshold below which the test fails conv_thresh = 1.8 -# Convergence rate above which a warning is issued -conv_max = 2.2 +# Run duration in days +run_duration = 0.4166667 [vertical_grid] @@ -52,3 +52,31 @@ partial_cell_type = None # The minimum fraction of a layer for partial cells min_pc_fraction = 0.1 + +# config options for spherical convergence tests +[convergence] + +# Evaluation time for convergence analysis (in days) +convergence_eval_time = ${manufactured_solution:run_duration} + +# Convergence threshold below which a test fails +convergence_thresh = ${manufactured_solution:conv_thresh} + +# Type of error to compute +error_type = l2 + + +# config options for spherical convergence tests +[convergence_forward] + +# time integrator: {'split_explicit', 'RK4'} +time_integrator = RK4 + +# RK4 time step per resolution (s/km), since dt is proportional to resolution +rk4_dt_per_km = ${manufactured_solution:dt_per_km} + +# Run duration in days +run_duration = ${manufactured_solution:run_duration} + +# Output interval in days +output_interval = ${manufactured_solution:run_duration} From 0ebf264fde7b6a16e78abccb4981281ba1ddcca4 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Sun, 8 Oct 2023 17:03:14 -0600 Subject: [PATCH 10/17] Fixup mesh parameters in convergence steps Co-authored-by: Xylar Asay-Davis --- polaris/ocean/convergence/forward.py | 7 +-- .../ocean/convergence/spherical/forward.py | 44 ------------------- polaris/ocean/tasks/cosine_bell/__init__.py | 2 +- polaris/ocean/tasks/cosine_bell/forward.py | 4 +- .../tasks/inertial_gravity_wave/forward.py | 3 +- .../ocean/tasks/inertial_gravity_wave/init.py | 2 - .../ocean/tasks/inertial_gravity_wave/viz.py | 7 +-- .../tasks/manufactured_solution/forward.py | 3 +- .../ocean/tasks/manufactured_solution/init.py | 2 - 9 files changed, 15 insertions(+), 59 deletions(-) diff --git a/polaris/ocean/convergence/forward.py b/polaris/ocean/convergence/forward.py index edd866ece..c482952cd 100644 --- a/polaris/ocean/convergence/forward.py +++ b/polaris/ocean/convergence/forward.py @@ -19,9 +19,10 @@ class ConvergenceForward(OceanModelStep): """ - def __init__(self, component, name, subdir, resolution, base_mesh, init, + def __init__(self, component, name, subdir, resolution, mesh, init, package, yaml_filename='forward.yaml', options=None, - output_filename='output.nc', validate_vars=None): + graph_filename='graph.info', output_filename='output.nc', + validate_vars=None): """ Create a new step @@ -74,7 +75,7 @@ def __init__(self, component, name, subdir, resolution, base_mesh, init, work_dir_target=f'{init.path}/initial_state.nc') self.add_input_file( filename='graph.info', - work_dir_target=f'{base_mesh.path}/graph.info') + work_dir_target=f'{mesh.path}/{graph_filename}') self.add_output_file(filename=output_filename, validate_vars=validate_vars) diff --git a/polaris/ocean/convergence/spherical/forward.py b/polaris/ocean/convergence/spherical/forward.py index b7c55cbb0..88be45217 100644 --- a/polaris/ocean/convergence/spherical/forward.py +++ b/polaris/ocean/convergence/spherical/forward.py @@ -19,50 +19,6 @@ class SphericalConvergenceForward(ConvergenceForward): """ - def __init__(self, component, name, subdir, resolution, base_mesh, init, - package, yaml_filename='forward.yaml', options=None, - output_filename='output.nc', validate_vars=None): - """ - Create a new step - - Parameters - ---------- - component : polaris.Component - The component the step belongs to - - name : str - The name of the step - - subdir : str - The subdirectory for the step - - resolution : float - The resolution of the (uniform) mesh in km - - package : Package - The package name or module object that contains the YAML file - - yaml_filename : str, optional - A YAML file that is a Jinja2 template with model config options - - options : dict, optional - A dictionary of options and value to replace model config options - with new values - - output_filename : str, optional - The output file that will be written out at the end of the forward - run - - validate_vars : list of str, optional - A list of variables to validate against a baseline if requested - """ - super().__init__(component=component, name=name, subdir=subdir, - resolution=resolution, base_mesh=base_mesh, - init=init, package=package, - yaml_filename=yaml_filename, - output_filename=output_filename, - validate_vars=validate_vars) - def compute_cell_count(self): """ Compute the approximate number of cells in the mesh, used to constrain diff --git a/polaris/ocean/tasks/cosine_bell/__init__.py b/polaris/ocean/tasks/cosine_bell/__init__.py index 34cc3e052..64433d1c1 100644 --- a/polaris/ocean/tasks/cosine_bell/__init__.py +++ b/polaris/ocean/tasks/cosine_bell/__init__.py @@ -162,7 +162,7 @@ def _setup_steps(self): else: forward_step = Forward(component=component, name=name, subdir=subdir, resolution=resolution, - base_mesh=base_mesh_step, + mesh=base_mesh_step, init=init_step) forward_step.set_shared_config(config, link=config_filename) self.add_step(forward_step, symlink=symlink) diff --git a/polaris/ocean/tasks/cosine_bell/forward.py b/polaris/ocean/tasks/cosine_bell/forward.py index 5212332d0..638b58134 100644 --- a/polaris/ocean/tasks/cosine_bell/forward.py +++ b/polaris/ocean/tasks/cosine_bell/forward.py @@ -7,7 +7,7 @@ class Forward(SphericalConvergenceForward): bell test case """ - def __init__(self, component, name, subdir, resolution, base_mesh, init): + def __init__(self, component, name, subdir, resolution, mesh, init): """ Create a new step @@ -34,7 +34,7 @@ def __init__(self, component, name, subdir, resolution, base_mesh, init): package = 'polaris.ocean.tasks.cosine_bell' validate_vars = ['normalVelocity', 'tracer1'] super().__init__(component=component, name=name, subdir=subdir, - resolution=resolution, base_mesh=base_mesh, + resolution=resolution, mesh=mesh, init=init, package=package, yaml_filename='forward.yaml', output_filename='output.nc', diff --git a/polaris/ocean/tasks/inertial_gravity_wave/forward.py b/polaris/ocean/tasks/inertial_gravity_wave/forward.py index 480328db0..599c22c56 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/forward.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/forward.py @@ -31,9 +31,10 @@ def __init__(self, component, name, resolution, subdir, init): """ super().__init__(component=component, name=name, subdir=subdir, - resolution=resolution, base_mesh=init, init=init, + resolution=resolution, mesh=init, init=init, package='polaris.ocean.tasks.inertial_gravity_wave', yaml_filename='forward.yaml', + graph_filename='culled_graph.info', output_filename='output.nc', validate_vars=['layerThickness', 'normalVelocity']) diff --git a/polaris/ocean/tasks/inertial_gravity_wave/init.py b/polaris/ocean/tasks/inertial_gravity_wave/init.py index 544e5176b..4935e912a 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/init.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/init.py @@ -5,7 +5,6 @@ from mpas_tools.planar_hex import make_planar_hex_mesh from polaris import Step -from polaris.io import symlink from polaris.mesh.planar import compute_planar_hex_nx_ny from polaris.ocean.resolution import resolution_to_subdir from polaris.ocean.tasks.inertial_gravity_wave.exact_solution import ( @@ -73,7 +72,6 @@ def run(self): ds_mesh = cull(ds_mesh, logger=logger) ds_mesh = convert(ds_mesh, graphInfoFileName='culled_graph.info', logger=logger) - symlink('culled_graph.info', 'graph.info') write_netcdf(ds_mesh, 'culled_mesh.nc') bottom_depth = config.getfloat('vertical_grid', 'bottom_depth') diff --git a/polaris/ocean/tasks/inertial_gravity_wave/viz.py b/polaris/ocean/tasks/inertial_gravity_wave/viz.py index 768ba101d..93ea2d8fb 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/viz.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/viz.py @@ -73,9 +73,10 @@ def run(self): rmse = [] error_range = None for i, res in enumerate(resolutions): - ds_mesh = xr.open_dataset(f'mesh_{res:g}km.nc') - ds_init = xr.open_dataset(f'init_{res:g}km.nc') - ds = xr.open_dataset(f'output_{res:g}km.nc') + mesh_name = resolution_to_subdir(res) + ds_mesh = xr.open_dataset(f'mesh_{mesh_name}.nc') + ds_init = xr.open_dataset(f'init_{mesh_name}.nc') + ds = xr.open_dataset(f'output_{mesh_name}.nc') ds['maxLevelCell'] = ds_init.maxLevelCell exact = ExactSolution(ds_init, config) diff --git a/polaris/ocean/tasks/manufactured_solution/forward.py b/polaris/ocean/tasks/manufactured_solution/forward.py index 198d38c1b..3890f2895 100644 --- a/polaris/ocean/tasks/manufactured_solution/forward.py +++ b/polaris/ocean/tasks/manufactured_solution/forward.py @@ -34,9 +34,10 @@ def __init__(self, component, name, resolution, subdir, init): """ super().__init__(component=component, name=name, subdir=subdir, - resolution=resolution, base_mesh=init, init=init, + resolution=resolution, mesh=init, init=init, package='polaris.ocean.tasks.manufactured_solution', yaml_filename='forward.yaml', + graph_filename='culled_graph.info', output_filename='output.nc', validate_vars=['layerThickness', 'normalVelocity']) diff --git a/polaris/ocean/tasks/manufactured_solution/init.py b/polaris/ocean/tasks/manufactured_solution/init.py index 7942b10b9..af6353a1d 100644 --- a/polaris/ocean/tasks/manufactured_solution/init.py +++ b/polaris/ocean/tasks/manufactured_solution/init.py @@ -5,7 +5,6 @@ from mpas_tools.planar_hex import make_planar_hex_mesh from polaris import Step -from polaris.io import symlink from polaris.mesh.planar import compute_planar_hex_nx_ny from polaris.ocean.resolution import resolution_to_subdir from polaris.ocean.tasks.manufactured_solution.exact_solution import ( @@ -73,7 +72,6 @@ def run(self): ds_mesh = cull(ds_mesh, logger=logger) ds_mesh = convert(ds_mesh, graphInfoFileName='culled_graph.info', logger=logger) - symlink('culled_graph.info', 'graph.info') write_netcdf(ds_mesh, 'culled_mesh.nc') bottom_depth = config.getfloat('vertical_grid', 'bottom_depth') From 74d52255867a40ccf926e671be69407d1269c5da Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Sun, 8 Oct 2023 17:55:58 -0500 Subject: [PATCH 11/17] Normalize all error types --- polaris/ocean/convergence/analysis.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/polaris/ocean/convergence/analysis.py b/polaris/ocean/convergence/analysis.py index 0b82f1851..0821bb2ff 100644 --- a/polaris/ocean/convergence/analysis.py +++ b/polaris/ocean/convergence/analysis.py @@ -292,16 +292,16 @@ def compute_error(self, mesh_name, variable_name, zidx=None, zidx=zidx) diff = field_exact - field_mpas + # Only the L2 norm is area-weighted if error_type == 'l2': area = area_for_field(ds_mesh, diff) + field_exact = field_exact * area diff = diff * area - error = np.linalg.norm(diff, ord=norm_type[error_type]) - if error_type == 'l2': - field_exact = field_exact * area - den_l2 = np.linalg.norm(field_exact, ord=norm_type[error_type]) - error = np.divide(error, den_l2) + # Normalize the error norm by the vector norm of the exact solution + den = np.linalg.norm(field_exact, ord=norm_type[error_type]) + error = np.divide(error, den) return error From 5a145cc1bc87bc71b0170e904e5f6fc642b2eaca Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Fri, 6 Oct 2023 16:07:37 -0500 Subject: [PATCH 12/17] Update docs --- docs/developers_guide/ocean/api.md | 35 +++++--- docs/developers_guide/ocean/framework.md | 79 +++++++++--------- .../ocean/tasks/cosine_bell.md | 8 +- .../ocean/tasks/inertial_gravity_wave.md | 19 +++-- .../ocean/tasks/manufactured_solution.md | 35 ++++---- docs/users_guide/ocean/tasks/cosine_bell.md | 23 ++--- .../inertial_gravity_wave_convergence.png | Bin 45708 -> 43697 bytes .../ocean/tasks/inertial_gravity_wave.md | 7 +- .../ocean/tasks/manufactured_solution.md | 9 +- polaris/ocean/tasks/cosine_bell/forward.py | 2 +- .../tasks/inertial_gravity_wave/forward.py | 7 +- .../tasks/manufactured_solution/forward.py | 5 +- 12 files changed, 125 insertions(+), 104 deletions(-) diff --git a/docs/developers_guide/ocean/api.md b/docs/developers_guide/ocean/api.md index c4ee9fd5c..7043f7642 100644 --- a/docs/developers_guide/ocean/api.md +++ b/docs/developers_guide/ocean/api.md @@ -66,7 +66,7 @@ InertialGravityWave analysis.Analysis - analysis.Analysis.run + analysis.Analysis.exact_solution exact_solution.ExactSolution exact_solution.ExactSolution.ssh @@ -74,7 +74,6 @@ forward.Forward forward.Forward.compute_cell_count - forward.Forward.dynamic_model_config init.Init init.Init.run @@ -127,7 +126,7 @@ ManufacturedSolution analysis.Analysis - analysis.Analysis.run + analysis.Analysis.exact_solution exact_solution.ExactSolution exact_solution.ExactSolution.ssh @@ -135,7 +134,6 @@ forward.Forward forward.Forward.compute_cell_count - forward.Forward.dynamic_model_config init.Init init.Init.run @@ -169,6 +167,27 @@ ## Ocean Framework +### Convergence Tests + +```{eval-rst} +.. currentmodule:: polaris.ocean.convergence + +.. autosummary:: + :toctree: generated/ + + ConvergenceForward + ConvergenceForward.compute_cell_count + ConvergenceForward.dynamic_model_config + + ConvergenceAnalysis + ConvergenceAnalysis.compute_error + ConvergenceAnalysis.convergence_parameters + ConvergenceAnalysis.exact_solution + ConvergenceAnalysis.get_output_field + ConvergenceAnalysis.plot_convergence + ConvergenceAnalysis.run + ConvergenceAnalysis.setup +``` ### Spherical Convergence Tests ```{eval-rst} @@ -178,13 +197,7 @@ :toctree: generated/ SphericalConvergenceForward - SphericalConvergenceAnalysis - SphericalConvergenceAnalysis.compute_error - SphericalConvergenceAnalysis.convergence_parameters - SphericalConvergenceAnalysis.exact_solution - SphericalConvergenceAnalysis.get_output_field - SphericalConvergenceAnalysis.run - SphericalConvergenceAnalysis.plot_convergence + SphericalConvergenceForward.compute_cell_count ``` ### Ocean Model diff --git a/docs/developers_guide/ocean/framework.md b/docs/developers_guide/ocean/framework.md index ab0b524d3..dfcea88ca 100644 --- a/docs/developers_guide/ocean/framework.md +++ b/docs/developers_guide/ocean/framework.md @@ -136,14 +136,14 @@ The function {py:func}`polaris.ocean.mesh.spherical.add_spherical_base_mesh_step returns a step for for a spherical `qu` or `icos` mesh of a given resolution (in km). The step can be shared between tasks. -(dev-ocean-spherical-convergence)= +(dev-ocean-convergence)= -## Spherical Convergence Tests +## Convergence Tests Several tests that are in Polaris or which we plan to add are convergence -tests on {ref}`dev-ocean-spherical-meshes`. The ocean framework includes -shared config options and a base class for forward steps that are expected -to be useful across these tests. +tests on {ref}`dev-ocean-spherical-meshes` and planar meshes. +The ocean framework includes shared config options and base classes for +forward and analysis steps that are expected to be useful across these tests. The shared config options are: ```cfg @@ -156,8 +156,10 @@ icos_resolutions = 60, 120, 240, 480 # a list of quasi-uniform mesh resolutions (km) to test qu_resolutions = 60, 90, 120, 150, 180, 210, 240 -# Evaluation time for convergence analysis (in days) -convergence_eval_time = 1.0 +[convergence] + +# Evaluation time for convergence analysis (in hours) +convergence_eval_time = 24.0 # Convergence threshold below which a test fails convergence_thresh = 1.0 @@ -165,9 +167,8 @@ convergence_thresh = 1.0 # Type of error to compute error_type = l2 - -# config options for spherical convergence forward steps -[spherical_convergence_forward] +# config options for convergence forward steps +[convergence_forward] # time integrator: {'split_explicit', 'RK4'} time_integrator = RK4 @@ -182,17 +183,19 @@ split_dt_per_km = 30.0 # since btr_dt is proportional to resolution btr_dt_per_km = 1.5 -# Run duration in days -run_duration = ${spherical_convergence:convergence_eval_time} +# Run duration in hours +run_duration = ${convergence:convergence_eval_time} -# Output interval in days +# Output interval in hours output_interval = ${run_duration} ``` The first 2 are the default resolutions for icosahedral and quasi-uniform -base meshes, respectively. The `convergence_eval_time` will generally be -modified by each test case. The `convergence_thresh` will also be modified by -each test case, and will depend on the numerical methods being tested. The -`error_type` is the L2 norm by default. +base meshes, respectively. + +The `convergence_eval_time` will generally be modified by each test case. The +`convergence_thresh` will also be modified by each test case, and will depend +on the numerical methods being tested. The `error_type` is the L2 norm by +default. The L-infty norm, `inf`, is also supported. `time_integrator` will typically be overridden by the specific convergence task's config options, and indicates which time integrator to use for the @@ -201,7 +204,7 @@ forward run. Depending on the time integrator, either `rk4_dt_per_km` or mesh resolution (proportional to the cell size). For split time integrators, `btr_dt_per_km` will be used to compute the barotropic time step in a similar way. The `run_duration` and `output_interval` are typically the same, and -they are given in days. +they are given in hours. Each convergence test can override these defaults with its own defaults by defining them in its own config file. Convergence tests should bring in this @@ -218,6 +221,8 @@ def add_cosine_bell_tasks(component): filepath = f'spherical/{prefix}/cosine_bell/cosine_bell.cfg' config = PolarisConfigParser(filepath=filepath) + config.add_from_package('polaris.ocean.convergence', + 'convergence.cfg') config.add_from_package('polaris.ocean.convergence.spherical', 'spherical.cfg') config.add_from_package('polaris.ocean.tasks.cosine_bell', @@ -228,9 +233,13 @@ In addition, the {py:class}`polaris.ocean.convergence.spherical.SphericalConverg step can serve as a parent class for forward steps in convergence tests. This parent class takes care of setting the time step based on the `dt_per_km` config option and computes the approximate number of cells in the mesh, used -for determining the computational resources required, using a heuristic -appropriate for approximately uniform spherical meshes. A convergence test's -`Forward` step should descend from this class like in this example: +for determining the computational resources required. When convergence tests +are run on spherical meshes, +the {py:class}`polaris.ocean.convergence.spherical.SphericalConvergenceForward` +should be invoked and overrides the `compute_cell_count` method with a +heuristic appropriate for approximately uniform spherical meshes. A +convergence test's `Forward` step should descend from this class like in this +example: ```python from polaris.ocean.convergence.spherical import SphericalConvergenceForward @@ -242,7 +251,7 @@ class Forward(SphericalConvergenceForward): bell test case """ - def __init__(self, component, name, subdir, resolution, base_mesh, init): + def __init__(self, component, name, subdir, resolution, mesh, init): """ Create a new step @@ -260,7 +269,7 @@ class Forward(SphericalConvergenceForward): resolution : float The resolution of the (uniform) mesh in km - base_mesh : polaris.Step + mesh : polaris.Step The base mesh step init : polaris.Step @@ -269,7 +278,7 @@ class Forward(SphericalConvergenceForward): package = 'polaris.ocean.tasks.cosine_bell' validate_vars = ['normalVelocity', 'tracer1'] super().__init__(component=component, name=name, subdir=subdir, - resolution=resolution, base_mesh=base_mesh, + resolution=resolution, mesh=mesh, init=init, package=package, yaml_filename='forward.yaml', output_filename='output.nc', @@ -286,9 +295,9 @@ analyze. The `validate_vars` are a list of variables to compare against a baseline (if one is provided), and can be `None` if baseline validation should not be performed. -The `base_mesh` step should be created with the function described in +The `mesh` step should be created with the function described in {ref}`dev-ocean-spherical-meshes`, and the `init` step should produce a file -`initial_state.nc` that will be the initial condition for the forward run. +`init.nc` that will be the initial condition for the forward run. The `forward.yaml` file should be a YAML file with Jinja templating for the time integrator, time step, run duration and output interval, e.g.: @@ -318,11 +327,11 @@ omega: - normalVelocity - layerThickness ``` -`SphericalConvergenceForward` takes care of filling in the template based +`ConvergenceForward` takes care of filling in the template based on the associated config options (first at setup and again at runtime in case the config options have changed). -In addition, the {py:class}`polaris.ocean.convergence.spherical.SphericalConvergenceAnalysis` +In addition, the {py:class}`polaris.ocean.convergence.ConvergenceAnalysis` step can serve as a parent class for analysis steps in convergence tests. This parent class computes the error norm for the output from each resolution's forward step. It also produces the convergence plot. @@ -331,12 +340,11 @@ This is an example of how a task's analysis step can descend from the parent class: ```python -class Analysis(SphericalConvergenceAnalysis): +class Analysis(ConvergenceAnalysis): """ A step for analyzing the output from the cosine bell test case """ - def __init__(self, component, resolutions, icosahedral, subdir, - dependencies): + def __init__(self, component, resolutions, subdir, dependencies): """ Create the step @@ -348,10 +356,6 @@ class Analysis(SphericalConvergenceAnalysis): resolutions : list of float The resolutions of the meshes that have been run - icosahedral : bool - Whether to use icosahedral, as opposed to less regular, JIGSAW - meshes - subdir : str The subdirectory that the step resides in @@ -363,18 +367,17 @@ class Analysis(SphericalConvergenceAnalysis): 'zidx': 0}] super().__init__(component=component, subdir=subdir, resolutions=resolutions, - icosahedral=icosahedral, dependencies=dependencies, convergence_vars=convergence_vars) ``` Many tasks will also need to override the -{py:meth}`polaris.ocean.convergence.spherical.SphericalConvergenceAnalysis.exact_solution()` +{py:meth}`polaris.ocean.convergence.ConvergenceAnalysis.exact_solution()` method. If not overridden, the analysis step will compute the difference of the output from the initial state. In some cases, the child class will also need to override the -{py:meth}`polaris.ocean.convergence.spherical.SphericalConvergenceAnalysis.get_output_field()` +{py:meth}`polaris.ocean.convergence.ConvergenceAnalysis.get_output_field()` method if the requested field is not available directly from the output put rather needs to be computed. The default behavior is to read the requested variable (the value associate the `'name'` key) at the time index closest to diff --git a/docs/developers_guide/ocean/tasks/cosine_bell.md b/docs/developers_guide/ocean/tasks/cosine_bell.md index e825e39f4..4d4cdef2f 100644 --- a/docs/developers_guide/ocean/tasks/cosine_bell.md +++ b/docs/developers_guide/ocean/tasks/cosine_bell.md @@ -18,7 +18,7 @@ model config options related to drag and default horizontal and vertical momentum and tracer diffusion, as well as defining `mesh`, `input`, `restart`, and `output` streams. This file has Jinja templating that is used to update model config options based on Polaris config options, see -{ref}`dev-ocean-spherical-convergence`. +{ref}`dev-ocean-convergence`. ### base_mesh @@ -36,16 +36,16 @@ tracer distributed in a cosine-bell shape. The class {py:class}`polaris.ocean.tasks.cosine_bell.forward.Forward` descends from {py:class}`polaris.ocean.convergence.spherical.SphericalConvergenceForward`, and defines a step for running MPAS-Ocean from an initial condition produced in -an `init` step. See {ref}`dev-ocean-spherical-convergence` for some relevant +an `init` step. See {ref}`dev-ocean-convergence` for some relevant discussion of the parent class. The time step is determined from the resolution -based on the `dt_per_km` config option in the `[spherical_convergences]` +based on the `dt_per_km` config option in the `[convergence_forward]` section. Other model config options are taken from `forward.yaml`. ### analysis The class {py:class}`polaris.ocean.tasks.cosine_bell.analysis.Analysis` descends from -{py:class}`polaris.ocean.convergence.spherical.SphericalConvergenceAnalysis`, +{py:class}`polaris.ocean.convergence.ConvergenceAnalysis`, and defines a step for computing the error norm (L2) for the results at each resolution, saving them in `convergence_tracer1.csv` and plotting them in `convergence_tracer1.png`. diff --git a/docs/developers_guide/ocean/tasks/inertial_gravity_wave.md b/docs/developers_guide/ocean/tasks/inertial_gravity_wave.md index 100576758..751be935f 100644 --- a/docs/developers_guide/ocean/tasks/inertial_gravity_wave.md +++ b/docs/developers_guide/ocean/tasks/inertial_gravity_wave.md @@ -48,10 +48,11 @@ the exact solution. The tracer and coriolis fields are uniform in space. ### forward The class {py:class}`polaris.ocean.tasks.inertial_gravity_wave.forward.Forward` -defines a step for running MPAS-Ocean from the initial condition produced in the -`init` step. Namelist and streams files are updated in -{py:meth}`polaris.ocean.tasks.inertial_gravity_wave.forward.Forward.dynamic_model_config()` -with time steps determined algorithmically based on config options. The number +descends from + {py:class}`polaris.ocean.convergence.ConvergenceForward` +and runs MPAS-Ocean from the initial condition produced in the `init` step. +Namelist and streams files are updated by the parent class with time steps +determined algorithmically based on config options. The number of cells is approximated from config options in {py:meth}`polaris.ocean.tasks.inertial_gravity_wave.forward.Forward.compute_cell_count()` so that this can be used to constrain the number of MPI tasks that Polaris @@ -65,14 +66,14 @@ performed against a baseline if one is provided when calling ### analysis The class -{py:class}`polaris.ocean.tasks.inertial_gravity_wave.analysis.Analysis` defines -a step for computing the root mean-square-error from the final simulated field +{py:class}`polaris.ocean.tasks.inertial_gravity_wave.analysis.Analysis` +descends from {py:class}`polaris.ocean.convergence.ConvergenceAnalysis` +a step for computing the error from the final simulated field and the exact solution. It uses the config options to determine whether the convergence rate falls within acceptable bounds. ### viz The class {py:class}`polaris.ocean.tasks.inertial_gravity_wave.viz.Viz` defines -a step for visualization. It produces two plots: the convergence of the RMSE -with resolution and a plan-view of the simulated, exact, and (simulated - exact) -SSH fields. +a step for visualization. It produces a plan-view figure of the simulated, +exact, and (simulated - exact) SSH fields. diff --git a/docs/developers_guide/ocean/tasks/manufactured_solution.md b/docs/developers_guide/ocean/tasks/manufactured_solution.md index 14f00e49a..0e5bee95d 100644 --- a/docs/developers_guide/ocean/tasks/manufactured_solution.md +++ b/docs/developers_guide/ocean/tasks/manufactured_solution.md @@ -45,31 +45,32 @@ the exact solution. The tracer and coriolis fields are uniform in space. ### forward The class {py:class}`polaris.ocean.tasks.manufactured_solution.forward.Forward` -defines a step for running MPAS-Ocean from the initial condition produced in -the `init` step. Namelist and streams files are updated in -{py:meth}`polaris.ocean.tasks.manufactured_solution.forward.Forward.dynamic_model_config()` -with time steps determined algorithmically based on config options. The -number of cells is approximated from config options in -{py:meth}`polaris.ocean.tasks.manufactured_solution.forward.Forward.compute_cell_count()` -so that this can be used to constrain the number of MPI tasks that Polaris tasks -have as their target and minimum (if the resources are not explicitly -prescribed). For MPAS-Ocean, PIO namelist options are modified and a -graph partition is generated as part of `runtime_setup()`. Next, the ocean -model is run. Finally, validation of `temperature`, `layerThickness` and +descends from + {py:class}`polaris.ocean.convergence.ConvergenceForward` +and runs MPAS-Ocean from the initial condition produced in the `init` step. +Namelist and streams files are updated by the parent class with time steps +determined algorithmically based on config options. The number +of cells is approximated from config options in +{py:meth}`polaris.ocean.tasks.inertial_gravity_wave.forward.Forward.compute_cell_count()` +so that this can be used to constrain the number of MPI tasks that Polaris +tasks have as their target and minimum (if the resources are not explicitly +prescribed). For MPAS-Ocean, PIO namelist options are modified and a graph +partition is generated as part of `runtime_setup()`. Then, the ocean model +is run. Finally, validation of `temperature`, `layerThickness` and `normalVelocity` are performed against a baseline if one is provided when calling {ref}`dev-polaris-setup`. ### analysis The class {py:class}`polaris.ocean.tasks.manufactured_solution.analysis.Analysis` -defines a step for computing the root mean-square-error from the final -simulated field and the exact solution. It uses the config options to determine -whether the convergence rate falls within acceptable bounds. +descends from {py:class}`polaris.ocean.convergence.ConvergenceAnalysis` +a step for computing the error from the final simulated field +and the exact solution. It uses the config options to determine whether the +convergence rate falls within acceptable bounds. ### viz The class {py:class}`polaris.ocean.tasks.manufactured_solution.viz.Viz` -defines a step for visualization. It produces two plots: the convergence of the -RMSE with resolution and a plan-view of the simulated, exact, and (simulated - -exact) SSH fields. +defines a step for visualization. It produces a plan-view figure of the +simulated, exact, and (simulated - exact) SSH fields. diff --git a/docs/users_guide/ocean/tasks/cosine_bell.md b/docs/users_guide/ocean/tasks/cosine_bell.md index 18b1898d7..7955e139f 100644 --- a/docs/users_guide/ocean/tasks/cosine_bell.md +++ b/docs/users_guide/ocean/tasks/cosine_bell.md @@ -167,8 +167,8 @@ alter this before setup (in a user config file) or before running the task (in the config file in the work directory). ```cfg -# config options for spherical convergence tests -[spherical_convergence_forward] +# config options for convergence tests +[convergence_forward] # time integrator: {'split_explicit', 'RK4'} time_integrator = RK4 @@ -182,12 +182,12 @@ period for advection to make a full rotation around the globe, 24 days: ```cfg # config options for spherical convergence tests -[spherical_convergence_forward] +[convergence_forward] -# Run duration in days +# Run duration in hours run_duration = ${cosine_bell:vel_pd} -# Output interval in days +# Output interval in hours output_interval = ${cosine_bell:vel_pd} ``` @@ -220,8 +220,8 @@ radius = 2123666.6667 # hill max of tracer psi0 = 1.0 -# time (days) for bell to transit equator once -vel_pd = 24.0 +# time (hours) for bell to transit equator once +vel_pd = 576.0 # convergence threshold below which the test fails convergence_thresh = 1.8 @@ -254,8 +254,8 @@ norm_args = {'vmin': 0., 'vmax': 1.} The 7 options from `temperature` to `vel_pd` are used to control properties of the cosine bell and the rest of the sphere, as well as the advection. -The options `qu_conv_thresh` to `icos_conv_max` are thresholds for determining -when the convergence rates are not within the expected range. +The option `convergence_thresh` is a threshold for determining +when the convergence rates are not above a minimum convergence rate. The options in the `cosine_bell_viz` section are used in visualizing the initial and final states on a lon-lat grid for `cosine_bell/with_viz` tasks. @@ -266,11 +266,14 @@ of these config options can be changed here: ```cfg # config options for spherical convergence tests -[spherical_convergence] +[convergence] # Evaluation time for convergence analysis (in days) convergence_eval_time = ${cosine_bell:vel_pd} +# Convergence threshold below which a test fails +convergence_thresh = ${cosine_bell:convergence_thresh} + # Type of error to compute error_type = l2 ``` diff --git a/docs/users_guide/ocean/tasks/images/inertial_gravity_wave_convergence.png b/docs/users_guide/ocean/tasks/images/inertial_gravity_wave_convergence.png index 320ad561cb3e88d45d62bbc974905259e922804c..a07e766243f47a88084ecb14fd537c49d0564e31 100644 GIT binary patch literal 43697 zcmbSz2RxQ--~UO4N}^QCrX($jtRk|~LP*Lkd+%h1jAT?IMOj%1Wfa+~r6Hq|l~G8@ zND2S%q5FN__j%v{_P-|xGAM~K={w zrYM@X40QO3;WbWC{EwuwqK@-%J5y&jBgeDU5hG`NYddFa3u7MFvyM&{cDB2BNbC^a z#$)d6Z0{tsbEnO}Ua-T?(QKz!UO^V#WU0N%NhgY0X+-|eq|2sRP*f(@K?ONY_gmlE z-8_#^{-z&&{mioXXUb|_^Fy1@9yh0MTu3^>s#&v@=YXQUeeIG*a!cqHm(VjQe4pe? z61=%qFl+xh&ixjL4*j$pP>u8S@v*;Za6sd$%#q(u404|2fAMtsu`%`8)i0j6N4nWe z7z3E)8F&{{JYCWh{v{x?Bj5`E!oRE)l;k6SY2o+(`-6DLUq8-%@|mBtcXg$!cNm?T z>X+UqFBORQ`_8j@Gd)Lr``FsGYwO%%;~l%JmTc3>rM$hpg~i2#GBUQ%F|zrd8}ILa zlX$>0NZCU&RU?F+%blWU;;8WY*&ZZj)nZs1y2h=vGxbo6u+N;Q`}yfsv#x+v($1>N@kz0|IBa&3ou&x*3@O{JQyu6r8ay~}eS zGxq$|7h_ZWI`Ke|ut^=GsJOVv=b8}rE#Sk+5NfukHDS6uZRU;K)7Xn|~~@ z`o)xdsC4f4(EEIg)m2qhqeD66gPysfqN0Z5{SAUUcYfI-Z_Kc4*)l4BY+qU zA&2QvqcCn+eJgGgM*5Xp0cwIx)-%0Vk44_P6}DvgTEn(N8>!}mZ{Ms%#l&1bpJaAq zS7NZx!Z%w~arf|z#{L<~P+eX9Tis=L>$yL_UViv+!-1ndoKw1v{GoF3v~rK99nE`o zb@Og3E31|t$MWvq56?dN@ZyRMqQN&}V}(zhI(a(`tXl3FQSb9i>fhsr zeu(5f5Q_zIDKxa$BOz(k)KJ!28#x<&2DH1Ox>|;e#M`&IFNKG< z%skN&TGLkMCCx7Fy6s_Mq3PYbcVD))#%$Mr602A4U0R*@>EJ1!dVPI;M&DQ2eM9F; zaM^Ueffp|pN+wvmyDLI%hr&)CGfWBlUc>!>$PO#vaRH-CP4V>jEb zGYs1Ibo{PdS@rE>!E&sa8tkOeF8`H-dpm--y@O5a!@IF9s#6YL6;0#UedsmgDj2Q0 zP5%kot5>hwIIN6~lO~@ig>!_we0liknO%5GdPzyi>B&}uJiAUA4k>4bWTnfEl65b3 zmWRd0vb~HIS#)*tvATqP=W|_$b)WV*KAHJmLiy5jdomx2*-?vCHjz*_Gmv(hqOiu_ zS(msDZ#c7S{g3f+k#j%Jyo?iDzG>5@=@Fm#ON;24x}HYyvxyi7v^d1Y#R(f#EpBRQ zG46kPBT{8!Tx(Z@cxO$UUHOU=Cr<30op{F_qA>m1XU;z?jIlK5h@#>Wtnr|xCe1KT z>9xo29+q=)5%Zp(b;IW2n*a6MN7$-mTl`Ezy!1M=*9o=x7VmgimyO|b*N+7#+RF(a|v%KRFyPwI(4UAwIOvb@&%&j!qtJ zNJz-^+{_5enl(#q9Hyh(ZUvdy*sMrY3$Ch?Z)$F?{PgM1eUT3zKG0CRcJE%{9nLMn z>}h5eFkw3Q{4Gkv`#3SHiH^Cs8N2M6E}zZ>a+MGI`%SaG2A{J~Q~0`H_ja7JE3Kt+ z>$%Q{WoQaDOr@l3JhCR^?p>vadU}$(#Kp~?J$q*3t-tF=7>!}jN^aY6FE*2J8EFZ| z)ce8atb)Rk&hpzdhA7^;%{(rBPwVT0*08dg7sgz_zKflm-PFS38Sb=b8~6EdB6Pm- z(r!nGrX3&O!Y1Cma^=eG*0U2|-yYt$D$k#OrRC*v&tJP9=I5WOtgO^+YmN>_&3O3x zw@3U8L!f;Bz`(QT&!3C7p@@9ws`Q`K4(IZkP^a_tni($AZ9|EO4GN-na&|tmV#SKM z;X%u~U>1HI7ncIl^GZr|t3wp5-`tVYtq%zgH>KX!gsgV?oSlsx%UN{q zVqx>WlMmS_(w7RYKUA$+z52%K9d<5%u=Lm{)3ax5=O%M}^7xg*b3;QzrDZPuF7O_! zU3Y4>t)QqV?K#WaOLJPxj`V}QUyl}jx$UDrBtyk2@EN71rb^Wx_lXP$3L5fyGB-Q3 zJ5T53j~_prr^rsGF_fK~PBy(MC8_(^@oY/$*wie0;Qee}G2_pX4crRDQ)!^0=_ z^x|L5Z8?5dn9@zNZ+iB2nEUapk(`fsgE+{gG9Rjl2n zb4R>moxZyt%^RuGyGO&XT)E<0m!ns5oyHJ-D&tid9>(BCWnA82LDmP@-Gg3#Udznj zG5h=(3o^a=d(q6-x5-jH%`>_J<$ivC~sV^McnywJ@n63+FnFU^VgJ!$>oa%9_E3o-eRhR?5U zd0+11bF+SCdb*;rGO>BKI%s9<2kqm>{qNn|ws`U4iuQIrv+~0N)h7(R0&8j%;$KO6 z)el*&U%%em+}!uDfd7pfH#(-7I3#@Wm*|$3!|c1x1ix~hbh?+B*=T)J)MC}?)2ByK zM+|%FF2@~wj62Ow4Ij5Odv%M3LeC$F--=4XzoQpozF|IQy40^t>%H*vFkwFq7HYm6es%HtYU7JQ`o<&mT|54Wj8sehwV=1TKFOF!efto*?|h*^I3DbpQ-#)*7nUrKynTCRh(hD0 zyr{CvYyCE+uz#{KlJuI%zITtJYHMpfMv^=3bG0Vam$b);@_vu<m}c3Z=vPivmunnTBC!%wYVyY}L1*?A7i9)+nSyr{6SuD{`CcVA!Ng9p1( zQc|1`olY4ch)cf;4HFMe+0^uI$I6=b?{z4*?~5(7HsAJKC7;NAs>rtaCS!hzzHD7* zC$~lh>)6;>dQQ$r;dE7XH47VC@V@g?ckbV3lAWEndN--*O7ZgHiGg?V9(xPkyiucy zq?X#H*|p^wRaZEij@dUiylvYyX8aORSc64gi?i0YeUs8m@}4%j zCFQb=+j~0kRo@P(>)VHi9#_K(J^F8$-O&`vv;MFPSk-*k3efJ>(V#+~Sx2$g_L*c=*;c zZq0)qvb3|S%DiSqdm?1{_U+?bz4@qug~b-v;g3rK<^6CKG}y&@KXgW(j5BRuXD8bN z+lFDaz~TLBYRu@_4<0@w`1o2(Oi#R9$JUFJG>e9(h8xvr$bjOPrYf%a<=l z+ub_lc)tE|yN^!Fm{ay!sO7%i^19epTf_=bXpR3&kHzjbW0LvRyIMg(VZTr`A8=LL zvD-?x9DH)yM65sEimI(SPcHoIIKL7N5Iba#0w0U(NXKPVrPIaE2fsYF`S4&#%ua*4 zOsy>Pjf|Vp)$weV^!2%RoPOk+bs{GzH#fIB`Z9PC-<~~d=jVQwqc8;lFRj$p)_&U1 zQ0updeiUur{7T>b(9{y2`elK%IpuRwxYkSR$x6JUqO5qK;Eo-2=u=VHlHAy9QI(V2bN*Is;>`AmVi=*f?Z%Oj!o}LfKkL{~P1-${BzgtAa zsH3C9`unHHQbY72u30P2|N7DQ^2LiwA+^T+U0oGi9$$G&w`!2zEzOJImQA2BM7Q>6 zd)_Q0B=jC+=hcM^7YJi`n43#5$+x%4+-FPNMp`yx>GRGUC0(7CmX>to@87?FOO4St zK6maMBULg#H{(**+S>YX1^b@d;9v$*TiaKMj~uxX5wY^7nAO^u?jY`1`NTJ2T17EK z?><_kYGk}`HSlQwihcI7p#gnp{rLE}$)W4pj2=IJ%o*Ny^y9d+t=-mgjf^usvyz+cH61Ku#r&)-hBc694tgG z$lp&s7)0%4I zy>Nln-Q9hZlaBnh6wl6`EAd#ezo+%4FKUXhd-~$VP+`{b)#M$mnwy#`s;hhM#~3A% zzfEhrckkX`m4fkQ``wy10oMOrxkz>m4i8rX-vSb;a}y+b#(+S+yeoYZ-s1MiZI*E{<-^{g}#q$h6@mRJ=d9PHAUgB{P z%McpDd$u)?zlA0&4<%o#M^N&!`azqa5Ow7`@=+<{_(9aw)ciyB|M7n~09qzoeuu-$xEp@? z9pqj#`r}K%ZX+8C_V$S=W?UXhCl4j?haV@u~sQ&Jeo1@b`5&k{t7qH+0oHb9$mXlTS zT_xPs?fWOmZKoeGZ`rbi^2N6328vUdpPMb|Tex16)|S_HQrRtsR&dL3tlM$A;`~S{ z{p!`L38s;glbimU9P#GUrAf7n#Kw&q4d3S)xeuhqsGT@b|Mv~m6MLWTx914* znH#^Ya_m_BElH<48X2RIiuffY*uaQ;rT5|~)R?O|o-w=?AD@Tz_~+M`+oKZ`3eZiU znp}Xqm6)G@jqIZUW`E^pQGUg!`%!UmpsY0=*PLGSZc0Gy-%O2waJnVuhKSY<~C7eGbL zlXK%04jg~4sy%$i-b1JbGb@)al{voKCa^b?OwyK8z7M7Kj<&g2z$r|>V z9_?Fx_F_cx!h0yNZ&F@N{k>jGZ0+su(vB%xS_(nIpz+)9O9VR_8k%X;_x)Ki zWGgL4QIP`L$^{g#oY{=qvsv9yL08QBy>O|`Jo_H&I}7itbmLCW7(L-%Yn-QG6$EVq zIQJ*2jdC{XwRcRFo76M?y@XNR&`kmIyr0_JD;+p&KRxn;jKaEZ-8g59UO~&0EL4%x zC9anj)?TD5?o|x>I0FMi1(3C&t!;Q+ozg97HzsHgbn1%5rNi^sB@9eVcNe}2U66_M zTyndmrKO>b&E>&COP!M^U4D&4t8M~k69 zY|OXt{6yHW9E@FEU4uhILkkPwmJaN+PC(!col`TZA`Hyi7}EnuC9eRqpP zMbXe=Z`?>sOr!w;lyh`^Ag391>(;G>H6JXOD5=c(rrjxA_uZ#YVG|S1y}i8@1z^X# zVZ$yR=EtVwrOpvuvsbkrxO*ZZb_yx1uMtucWg zWytuYxqKmulfC_Z$bt|X;WkK00&*&1)4Yx=Nf@3FFPfe&uONd zBWMU0@mcpEe9`(**Z{l#%+8LXB>3Nbipn-%i$b+!^XAP^4WZ-vi(Is7z0V7rtEHvo z8Pp;GBuHd$!nDDmG&ME5o^8He)BEhI!|+G3g%WREgRf|eiqhTPO$r>_0Y}VqQ3aOZ9&FvWVKYeLD*m7ZX&d1b!=NFPnMj zw;-aK-1&x)8GdR6#+F`IcT)4x@;8fByhcSzXzs~Vv!+tylTsz8Gm zixs*O7`SAQ?We!prFbzZv0eTB7r_gwYipMP_~C9!ocpm6Xzv$OcEk6ksT?)IX{HF= zAiq1fslqkdb5l=WzgAXNT?XBIF)a-p69lW=$RPNi_wQ^sDO(XKLgk= zRG)$B^yu5pI_2K8JXfD}efxG#OR*T8r!%qO1Jq9YuWxvdwE1XgT)8{k8JR$j^_igj zjOFbip=e|D28NX@gJ9KYX{~|yLY9R8jjx{`49|liH0RIKFIl2aYHXx7%Ki7J?d_bq zZZIA_dQ>e#Q{~X11E){llG8FkwaERCGTSI{jU%e$ifSt4b4shn| z+4YClEd9GhZ4NwRvDkg=*fA_x2`MT0R1ID+u{CK&<9DSWa?3|G29aJU+qAyu5l5l% zWACNh4b@av_apCRX2z~Lsn|Y+@9Y8pibK+IF#ry3!9KESRk)S2|8N?9!y0XV^M>eP z*P^3oC@6;5Uv#^7@5UPCCT|~E=BNAXV9{A!$UJtT|z^b{TaVKztDk29GkLRXz<6?#l;1iPyvOHs(AaB z#mLAgQZ;q2a*GU@@Xb9oRZz^KBQZb=EU;?*9Jo&3S7%%QLX)<6gLb8gdV3{#$|#Ii zTAF=eU;xVlMV)+eh|gGu?BuBgcPD9kH(CK9f<^NS2`wuuEHt;ZRYE5t@*etYS@l16 zrUwgjKm?VGZMP{Ss3{LKgPMlM1$?Pbj~x6*-h50w+cs~4{?$A+|FhHQZ`D1T)Ehh@ z7~a!0G<3?TDjt*D9bm>wEjkhYbNK+eyPSZV7KHri*d+fFhXN1ge!G$ZMN$w3E zt!ivs6%`dln#6vd{SFQy#>U3b-(I_!5ASgq)P%=q`Ed>i{2&2n#uY#`Pr4t+FC2AQ)^g$VdSR?@?1z^M(JEn424kdv~OLY98Ig>f;MQdw8^RDEzX~SFg91PI=iKQh?kq4T@Qng{BuHO>`t<3W{8=9$ z?uD9Iib|T&`T&|B97-rF^x-|YZxIAutywdMB?x_V;WOA6u!V92m>If8Mk1a*RngJa zP0Gkv3Jvg)`&g3PQTvwGsS*6uNZCv{K`bdTpw4WIp^k*Hi!;IUjT35w2dZ}Z^h-Px z*Yj;Y=&T7IgIgB%vi^lD{KmqlnpiNxrZ3sB@bBQ=0WJ@~T$9s|(zTkIhat!K;}Z}dnBNy#yell4G@E1+Fec2@_{EVLzOddt4eN=JB@p$IdZ z)P?yyek=)*%39j;lXd!PDFWndQyLU)?~ox{qj`G6u!YD z=K<9C`4|K5L;`OHXheBS7MaVtj+_7?LJ@3gZCyk?DJ^ATW~QQ4YMjn+I+2ipZaxDx z3FS+Sdnkw(a&mUQ1~m`vzO(d0%CUz<}Ns~F=H;td#0%D&xsjl+f7q+E9( zQH=*}4n2gTKv@YRXlc2Qr+2h{>bCmk!1kIn-Up6N-<%4r8#lb%e|UvspNPCXLpS{< z<@Vq!{5{7{pAPHoJtMkj&xIa|BWn)w2*2F@dgjWHv6rtF3%go*Uy1uYU*_(~ut^0N zG^n_l!Ghk*%P6H6u7%fKqaDsg7E3dxz(GGCBROi1H|v_F8BbwhgDm=!{M-z3%~KZALj-6AQ^2# zg>z(Iw{9IilNh;SBBx6|*3?IEE5lAfA2@h$IDR#MfzOWf31ojOtW*tWczfV4E@YE( z7TmaMDZmi{VW-PHAGe-}y0tCUX2XeFPSO8*`Qt^fp^3L_Yiny5)(vzf;qZAeDkEyC z7Nnal%k7F*D`q&Tbz4{O9Unczu$UHkgPnxPxpjKc55DFJV8*9D-WsBC0g@(xaC}jB z?%t)tRo0^zdYoNzP_Qf4)Rc$%x>05Dc58g{8 zx*P)yiXCrgfkc2fxz)wR7ep?;h}qc#h+X<1R8;Kg8iTzui+8kqOnZO&UaGFm9~KrC z;(hJ0eoq}gew;(vm64a1*Bvz(9b}pzo>|^qc?b(8+era#ebkjJw1M)3l|!ViYHdA{ z=8|CdFTJh9CGZ<;K>`cmLeo(J0RbeL!Y1SCe(d652Hy6IYyG|{;i*vOmX`94j=Nyd zQq-P3dl;54_fJg)Epg#6VZ?4(t$qJ=2#^oDGhi@L_5FZm$<|Rxd&3tQ;Jh2FfRTL< z|IJay#mXv?|`qD~JEczvVk@+i@$5#O;P?pL5r3|gCF5%(D z-n==^aq`RQNA@(A)^0}~p`{GrKuiy2<#c`dQtxyI8-EW`KNEl{!`MVWfz2$Irm5NH)_LGCkzucovL{CiV47p&Oi}vQPq`AhY0efNQV?* z)RZhHIVG%*#dLIh8#gW{&fGLp2-_B{fG`e8jpK@$>FN8et+&BiSy3F?1Pho$!k&`- zGr9<{71$*#ieK3n`L6Wr>^=g#2Cm7Rl!ODS0HA_M6M;wYD!YF4y|9_va}1ZpLM{!M z4#1PD0A9tV3Gd$BBl+OfmhIbXFYPP~=Jxr+LUna_`|T{7riIS-x9so92pLeLG(5rzxa}W&%q1HnCm`ZZAsA-o1OJ z){F7l^z`+SjK}%&=fSmMGs|Xr~-tOY42Tg?p`1x=S(> z4z#(s;2yhmd5p_rluqm1P=w#6?JXxi(nw9%e~`A-bexBalkL*1v|zP4S)Vz zm1OonJ8PP=PW&R&?qf>pYst?NOjx1p@K#el>@lukU|hdzG3?er(^C*L@zX$n*zzOw zHkOuE$Oq+`zvX**Q+y=c1mEcr?2rNLG)d)@tuPvSMMTt{EL#E5u{4Fn#QN|!Eo^P; zV6R>f@AO=aXD%!vQVF=ea9g&mCD1r!<*tvC%8rj0Bwj1U!8Bu*DUfsG<4oh`^2=OvoZtAw9wjI+B|! zjHj9})q$iUFsrV(eDSJ!hp(OGvIw&onwbS5mxEp04M!9$7(M>zAKdZ?0@HXXV?$|lCE7ttd@JZw-S7`}e;iThF5upr+yHPyJe$2UUrhOtFtlBz8I z56`*rk&q;YtZ&Qu$Q0UjGIh1^MMA;dZ)C)aB~BDxpzZt0pOaq%G0OS_sAPNX}c@CDJ&Vw@{{rTtY?S%_%+KBcw%CVleh7HQ^<#b1d>AM;k zqFNj4*nkW`cVLp$xtYh;u~j4AMm9AZXn$~$pP!uEzr?#4B34DTokgbkm-ercr~Wq0 zg&XtWLF2Y?55Qd{<;tpuSRsv~UL`vviq-|@)J4?VqQv1EVWmuO5ey@~dvU@V7lt%mQ^4I(g1D8>QqGDsy z9}7D|q+n)cU5p&2)0udnmmqLxV`Oj83Gu)xiSd}0;E~yIseF~ZO3^DxA=GCu_HekF zi}q&?U9sn&Pu2L)QT@{c6s-eNGC@60U%p&y!if6_@a~5UL1&^%$fb2Vjo&6If=uyI z`vtBV4Yg#Cj*s_9?t%jMwWaqf^!z#~kf8(b(oCQS+Pk;}s%&hO-XkQGSsFAl2}WD0 z2O)@ToYozPa4M^*MXKl50@(DBj?QlPBoTxg^>)2tbQp7^SFxBRnqn3581}tSqmC}I zNXJ`T#7gN#8+wGmA)5&(EiOvI)0XXaG%t4M8aqml2Gtd z7m@T$*!tZ>kZ-){Y}NObgo^{0uggWjV4_!H7?4CYBsLJ_36#jU@81VMG<(grWy>?c zOqJ|{S=k>?`2~f9KBG|w`wiG{^+6m8GJJPemv3YwbLo1vdWUK39VyN;cI2X-UR&@9 zBBP_1P$=q(XkF;1=fdwINnzCK51p810fH*Zs8_HV8?-KyH9O2P8!!D_gF8>euCgRs3(? zzJ;THQ6r-U>;b@2#N)5)(CE2z5N!-vJtVn>XW80nGW@Xs4}%~E1QAKMU`w=(kKE3j zpWUU9{U?A;2KgI!%>XB~6c!edE2X8SaY<5Osmh9qWJ`+{U18>bm%aLdmXOv>+qzgb zmd4iLT4N}}#!%6*;Yj}O^l1(_2x0wp08$(zc#A5t7|JIIAc3DySa$R}N|HwHp%Zi^ zxv0ZqIpNGF>VYh4p_CCNxx1T&BDS}U)zbMs$U11ib^q3F-x2}t!H@bHjY;E^s%UE1a`mY{%JF`>?> z8V+kG?5R4xuoUTrMc#K(f>|4iYVfXSg34 z=85_P$SwSvlc9F%)FlL^vHPmc|L)B{qtyV4!7M{q6TgRY_;EY%Y%7q50*G#Gc3_cb z3_wAJE%1ACh`x^PzJ6(FPkjV?b4(hs+zIKx)B+yr*zNg0^u#B{yS(+gXn-%Frd+qy2h9M<*EM9sLLO9QN$=uU`e{eyu#PunO+m>U~iW zxLR<|FCRkPKie=M|LBeweqbH^gT|RIwRi8Q^f_(3=@@d|pAlG6fQtq%VziFkju`hB z75KVQa#XgWFb@u)7WaOK$QXI$N>T>j3JwWoC}>JriqBuYs>7FX=&GQhNS=%IU66|( z=(`HHZS@vE$O-+ogvf*ks5W^NE-p(<7)j6yk6!)w@oG17y*g;uh(Q&>l10pQ47MGC zF?9FdMJFR`3E*b&0NE>HXaPA{pNFE+(N$*P7m+UoZI;$=|8d1#cZXu4qoccF3)ifa z{X60>W`mK2R75JiGbmt>Epa3M*P$i}-`Slb)3_izqe5 zGVWRU73Jm0p;*C;A{&kqu?GlvBm+Phf^bZKd-%HR^JjEnH1H$g!~!;|zL@yv5gQ&N z5pj`=g|-+5_)cUREI^_h_iHSwBNISU-;vu}kk~>+&EN1avI$qDXR4qO>0oXGwhu$9 z#!@Jio^Z`|at-~|_Z^dGJW?c}3}e0qv58CBr;3bGi=}4DfAzDI9KHKsMp~{-Tj3N} zSu513$X6ED)(V88U%Es?K^rF&9X$^eJv$&^F_GyI>7eYhhUcNuRY8bCdx|W}ICYgD zk#R5?!`BJ>vcqp6nDy$mYy?nc{@-dI;C%V-h>T^X%G}%>uywh*)(NtTKA7F^zAWAW zEw>!QF(fsG4%1j&TyQEItQqEWSx;=GKYen3em;ND=70%M94@jD933yHrplwIp*jZ^ z7fXS}?;@1}FqKpWP;(-@0u9lT$^bWmREDo_kq-fJ$;ilX0< zs!CyOeo$V675HFK@(2jf0U0W}SP<*J#-6l-U1yg;T3AT}z#cz+fP>NVWrbY^n~p9~ z##_L>N1iNRulZDnR^j#yu3$HTNbhiM{&CXHHAR5@4F%L3Z+fw-PHV> z11q~9zLxQaPD0*EqImKLFX#ZtvciwZkcnYKV_oUwu>BfBS8+>nB90Ep!iBM1oZgyD^M)9CusA8G(y+CzSAK5_wV2L z!$a_<0sl8Hs;9bfl~LdGXl$)oBzSmZUW@?+8zU4T;q_Y-V8r;1lqm1{s(_{GnVIq$ z8mj<-h@XU1LA$>bM>vOM<0ii;0a4zYksT2t7a=~Le(Xq$b##1WgTJ=saDxU(0#byo z*{Y_iyDYpX7ujjDnkKjuRhGAZ%x-BJnziRhbDQNq0FROPVltFpknHsP_vt7|w|xKN zrjxdIGb4CIgky>1_vcuVYgUG4=AQO25fQ!swzXTgZiO0d4703Df(ZnXNMUkvH01Mv zY^^)XID?otNgxM#;RvigZV&_|P;Bxog}+{T9l)HtnVCTtJd=1^kUk<`^gXigrbT4g z@Y?FwW}Er>->v2OALoi}+f{J^`80Z`vKN6_@?Jayf`usiBacZT5G<9AgkXKFWraYS z2j%FNtPd9@Z2L*Vy1jS}($k+o|2*`I610(`Ij{*IiRMK{wGi>fN0Q-+9r{nquW6;A zPv4U9ibFZGpwEQsx4(W+Eq#B zHm`Lh%f0?St?; zhU~BrMkMW>ofVKps|#E2Z+6>jEvD}-SVj)k+C5VaQ%62K{6dt%L{S~g)r>d&7)#orLIEQX*3QXh=EM@s=& z3%idR0(zs2aglW7=h{moK|xi>&W)|~M7mi?S}uOV2JTOCLab8_8_`kmzG^^%7-B-@ zHpY&oY_)FSZH?^gj;+}wUmUBE(WDtgktn>VXp{NQM|$WT1fFwyO|YOi!)VSks-{B< zG_u1WSiO6nWM)Q&99TPAU`fyDlfi+d9|jekVG{0P<_HpXFpQ}Xg(R23%RP^5BZ4B! zpsSGqI*^1#l)HxqpRn*7^DGr;8|A+Tw@%N@i~|SiC7%*FV}QN2+w4^UCh?Hz5_bE3 z8fK2x$a|Pxh%BUesjn`T2=^~7Ep6#z0-Aso4ZoP^Quw7VurCx(B{bG_9u<`@7{=p; zWQK$kJq2<@V)$T5kbw3yeVq{=X@z%)IRqN230a;8iABnMPWhB2y@QIja;^O0| zldG2n>O6Mbgh5ylc=}{)@>l;&GQ4^8woC&G+eGl5&B}xwAumo$_=YVHF}}J4QG{4= zJ1&TDL*K)$eno%~{~#}N0gQk&MOcc2Go)x_n4X6E2J)k(uYVNJQ%XfoQR75S4gz2B z!vO%`EdNw1nw|eM{u)e>%s`?MC%5)XDkIqLi?{?l(EP$LeSJp1Ck7z)DPx5KRWMQn zCXkujDTfj19&dPL)lu`w+ann#Zx6Em@2kfo6Z$~;m2K@NORg^++cTvNcv zm5615tv*7m|JRDJNY}w1lFGk+{UU>(uUvs48yqJf$X5fnyv++fSOn(-K9U5MInUZl z5EBD1W>ZsJ+YLJgSfN{5L7p)jGr(9+UBkmGx>}yUL!&0c)Z)3e#&&kZ#b<&kRfh}` zi36N`Xxg0^#3lc?1~&G6yv1?y)Ttt=h+5dEaF2=Vk>E9T64IgDN#7JF)C=%Xq7j2` zn2?DQg^avB^(@DYF)t98kxT1fuvi2c8eAneR!y7rUqe9vv}G zd}fCXn3n}klUY~5Hk8F^?TO7wsSRn(O|2b|-vvAYxjlSKf&~p6;SXgcrHkri2iY1e z3R(as5I?QL1+EHNT#;>8j>HcLwk|%K*R2JKNN{yH0h|;Gh@yVJ|L}o4--D{ERaJq@ z5gP-Ut8K|X1tV<2QqNm=IN|(vIRGX=t5cjNFr_dx$pik0$jtDi*a<s>_&R7uO&JK)6w5ZM>9nVR5V+EsWY}qD z{>YI4o`r4@!PQ`LE#T{CFF;q*m(SD>0S{u^3knKSCXCSg;=I6Ud@&{hPYAhL^5Ti3 z7)vIqs!31wVp1O2FrSUoE|%}Be})7;%|VnmG&Gl5bvbE#23zs^4!uwqD4pq@n?sCH z5PxHD5n-JS`m{T^>rY1ari;NDk?@?(?Kg9|q7br*_ahJ#ia=CM3_S&gx?e{pw(?cq z@k}_7T#LoF*ELW$mY0Y zFdzJJE<`T*#rCVNC{Q5U#1A8Zcwjl85_v&Z%N=hrTAG>^vIefmyhaWhUMhw{48ttS z4r^*@k&w-`rosfab{SjXI{-_FXvUEs$MJlu%W#nwF<4{#qxTu&k-m_+Ov``E$D`Gf z$fNqHQ!h;N14fl`wa|)?DDw{qS+L=M^d3YqhVAf!t}|OoXYQ=+!R*giB|Uc+);7tB zLg6N7IKb4z%bd;91Dr98k@5ln9L2~&PhKujZ9fll7bKrYj(y8e*96u1yyE^h(h4Rb z;zr(lw|o@Ac{)4Qd^qU<2X(?~%f%UUw|)%nf5r$PiLe3Tdym;+ z_=UtYPTwmymiJia&{%EH0b#MHli@#q{`j#lmk0f_8dTZB!Xh{A9@e?+I5_UpvhJ z+%I77^pK2|ug{+z?rH4M?m>?)Lae;Q&3(_JeqRV3M$>9eoS8I5d`N#;p!`)qh(dlM z#1X1^_G~Fd;)kgk2iI&e&i|W=NU|_T3QFR>SM>S=jYxSeU?d2u|AvurTeV*YcXyKl zr1VVvb<~SK$nqCTN=o`y?84w9bVW=Gw%qc7GTI}lAz$1wthEs$4oT)THeUJiwD<&A z2u>&nz+3@Etm8bFA??2)NrVo#j$bu=IvPb40!#LBJLD?toty|M+@xBQ4nJ-Y{Lx0u zD7qCJ0$>`lxvV(|Q;dky;CtZ^PE2D3*@RWYi&-`)|FJ?`LtUMP z++=25efq9FS`L_tvxC_Ymc$Hb4D z1>TrEz2(axWBZ|bQmDn4r3*qUM-18sll4w10nSIs^dz(#*sn;?;2S3+Kd&@<)>q^r z!9#?u10nrq@x9&QbFf4mrX3W8QZ|#gl?a+F16r_{0`#kFZeBw)MAbIe^^4^_ZvJO7 zqaxd<71yEEdFM{fun&CLSjbm74lLDJXt>;`!a!)OI5^>GX#xBrdLJL3_ib(Ctem{@ zR{=C3^c-2U1-w~CMn*{wp%7$XKZExo-YGlX<-f4X;21icsbMC5B2w;q(NW|2Xpq3f z*M_GQ-=1x?Urbcgw?`sVkj0;@Bcjt>zfSvAT5BWl67FbylHw(}r=RdxND8aE1#VAL z5}*6Ia~5T1(^68%!97z`Q*rmj!N3yuH>KE+ct4NJEM;ZrhDBo#~l@1hAi zhB;l+ee1hg2EPNX*VYEE+SEAI+}g^&Wy?}bNKh0=Bv#@2*~YbKn&2CCRyn0gHuJsr z|4KIN(=J}5Yse@9bZKyf&Tu$Z*suO6g9UqC?A_25O(DqB7^hCx+9amknqyG^;OvW0 z&?_1+pCs;o(qs^&at>Z2p=dTX1`1{*p`Vz`>i+y$zVq2g7 zXboF6In@NaH#sFm0XhgmbCo!>Wn@$ZamR%PW@qtu@BV`amm<=;ysITcoSvS(AwyjC znu=W&+W`zo;ou*zEyTj%37MQgoSVc5I=y~!Qe^TH@Dv3IU-2lWTJ-R+9y#>Tkxihd zm;K=tP&X%M_sGbxPPKzM8Zc7H<;C2^00%B_rY0#JziW$hG*(V9RV3rwG$e=tNs0nQ z6`dD==>z=;SoRm=Rc83w|9iz~7%v;0YP0J>nneLzvtafoM$BB2r#33%L13nWjH@Sp zW#C4g>2-YEQZ$A*E>e>BWdgk)FNtzD9-uTx4q-8|($6iFQ7PHnfkopFZ6yID+{C+ldjsQREB({#Ic( zp@N~N0G5M=w?S^lzfHRf4W~JDqR9YNoT;Wh>ZrwbPlDl7h+co0fmt(~+zM=kW}Zb4O3)jaD_X zupqo0UlHXX`kdfGUWIE zvD0S2`N-?nFCzGhVcEl3Z!b*f;m+=sa*=Op3b}%E-xAt3_*?_D1;}EHI^RC|*^b_S z!rYil{Tt>Ej)+h`eth+77T2BOJ;G+MIN($Ez=+vQbs^KVA6NiNf3=-`0A@E?s{F4X ze5fWE0KRMbU}`u->-;=f7nc8LB=o!&Fs(XwbUd4bx4 z!>G6x2tEEC!g2Rfq$5Xh>PHo9NWf=I@-`j3rmK;msj1G_j$9^u9V48jg{U%G9*IvU z=kU=~Ab~+njzTKv05q_l-`hzb5lP*p(W

?kPb;sNLDPOC5HJF9x*(MRkN`w~D5R*}fVb~9X9|M!xZ@-*d*m8OHV@5Lr`UN5mV0DWR3+Y)p5UYC z(>oUW4LT5_j;M5Gx&ab91Vz}9o9+4&6_VdCAW{b{x03DMEsx`yRhI=` zKw&xyE(Y^Bn561)yu&|9BKYAT{J)rT+5!Utx`5H9F>FSN5V6}ZG7qbAXtkConQ|~C z;zftoPdx|@Bt==2u33O0JIr|KQ2PxfbHN4AD;9g!sl^V7WUyiy^^(xU zi)qR62RnA`AW}R~BoT~(0C@!kpMG6;OcQK=rgOdynTM!^gmALF^@>?BWgw$>=&cI$ zY3D5H8R%O>rwuSU$B1QzWq(EL-~%k3-;o-EtN@zCv&CUp7$q6N^JazONn}9`7eHME zM$DjuL0oXW*YSqo-KcuKkRC8H8iC0PBsAv;r0&`_gIbgq@w8C*r6Sg{k}a`m&qo zvwM+fAifW%0ZCmr_Ng{#HMFhXNqNprwvs^)O`$-D7uwmnfpv9t+T7)Zg==wWStWEJ zK#FC+ROA>ZM2MT#Z@u~qix}=R64Z;Z$51oL*;OFbw#Tle&s=m4*w5pOLtS9m5{Lq& z3?mi%yu5O04;$G~|GID5l~p29h1e&z=V+yjh*9az#CMj%%$%H;5k1m%ySt((9c4$x zXe*$XUnZ$Me3(wYISUTXzzkO)7+`6H0XiSH6PbD<8o{+K$I;?0VqbYp_iy5kmKi7Z5zNDfi`e+1BO zSVuV12;drP? z&icOo{tgEpK#MGKq>uXZ8X6k(WD2tYdqhyvM=sYMl7gD2s0b{a)#=c6+UDeu)Vm{8-!v?b!f^5f&kWb0<@ z;Mgt`oQo5OgF$3IS(ap9zPtj*{AK<61;T(sncl1y$H^}xb)E4te?|{rZeg%J&F=h( z(=}qqLSl4~9EwIBaJE-Y>C}co+xAc#)w3QMq_;N$GIw^6|AU+F-W|BEtG74hgYDU~ zfxvQ{228b@{mAD##^~b$fpDrma-Gd+isYOIFf(m-TYQETI2uVt}`Y-hkw@h~x zZ7@0fFgNA&j-Ulz$b6|*^v2En5SU6>Fu*F{VVJ=u=_uX#@0X)7;f@b-A1-LMhaLb( z+TDD@pcW^xBo%Pbi*PIq4MD4xkkB!kLC=T3egRqxH#+Vk*e;O!LlD8-@ryd9qC!!n zPFF$VMKC%^Tt#9msHwft>p7yTN=s#*(LB;w2%d%P5lpQ?z4-ww3Cb78%*FSbV^T3& zHcuMW5Wa=3p56A(M#6EVC6HlW8< zPAn!EBpXOsLv%VZwh5-Ou5ZohIp< zJ?m?0+Z{O9cnz(0F-=f@p~_+scAcJw)6~o??c@=uty`C4W;;&k3V>{9w~I7NI+v(b zn6i9gZgbo4>fhf*3;X&u6$}3%?b{V0GUJ33934XJQW&TYk4IPmg;7MZ zD={e5+xOjeEVGm0kgBTgzUyW6K8VRzvRO?Qnl@~RpK;YnN;-NCy~q%HK{Q83*1_#1 z9gy%X&QVE}Ha|mt33*LYbLC)s^rNdOCrn^1Rb7(0nue40_Tk4 zrWQAdxTUp8zlZhQze;%UxShCNM@k<}us(gF=u!lS`_9O z7s+9Uij5?Zp<=oeW@QO)->i_ed2{C@hvbHOZbbCR$*D-Fu|th+`S=sTuH<|77%_8a zGagz&W&oitBEXTZwGUQGCC)HP-ra$LU4&7DNx%XBt7fOzj|4J0hTXW+peoYO$JKEN zU@L&1fX3j#jubn{cIbkQ7}m!MflBO^6uS72Pa8QL4P0jwcZ?T2ker0+mH4%m#Kg$J zMjeUai!uitP?FuK99NLi!mF;E1KKB0pei#ik$q_Oi2gjDs z(vU`lu!SL3E;$zo*(tK{Ire$5!b_;YQIK!?-MzfL7EIQLI5A@3g4Mu24-lZ9qoSbrtWTSY3NFBl!@1d2vhQM|M;siMb3I37UF_~W}9WIKP zo~eK>7pt?arR1)f6xmBKl>A>)1(y6>tlo)#|Aw1I!7g{rDq6Y;4)qzr1)IeZOjqD^g;F~ zQA-hO($H=a6CLq_Y#lIq_pIVG82nq+8;kZ*BKoMnqz;UZLal-ySPe;a2Giuz^M&&m zJivI^IVEy3ksC+79AA|ACROmUt&)}})d#n(2`p3kPaz)om7G`c=#d0?Snwq7F{T()QyXvnug1;(Loopm}cw|x>8Wv!nEHMVdF!O9ozgyd6GGdwh*1@sTVF@{1?UQKM@SBCwB-k zz2WO>xR)0!Cik(&;yKk-G&D9Q?lj()gnUe8vC5ny&M~@1%d?a+VWz;?&h+fi_`iGc zkk7p}mQ6gLFeG#ghjZG5ZVh!ffq9TNAmUz3Gb!4I{G5&CI?laxmoI7~DzMnq_cpY!o!pX~23(m$&Nudex| zC_b~7%Q~-LX2%@!Gc6L>`xm}VJ_bdRY_dl*r~HQU(EK;BD%9V+aIVPkwi~rpZh+{4 zdM~-#Wc2U;z1moMO^fXS-TT9Tcl9b5S7*G&NNQPor%SH45z*mqtJE4D znv&q8s15J8Tf3fC^`@AWv;g>;pXm>-=XW@qx>85A{7L=D z4berVETy|Vc<`VgoKyBo+ke=3y~6*&6PHOLGn*0PL@!H;h92f7kNloPhq6wc5wf1@ z8<$*(sKK_59%{7w+DK_$EG@f&7z{G~ee!9$ZvS2-#K}sOH@?*%ue3f;)&CkkmE%&lpBKTG2>!xWQFjRW1Hvbbq%_9p%|#h7ypVPa z1X@l(;ZX++Xw3mB5n9AB%-!`Gj9niCL*elAb0VA24U}w*MGY2C>sm#$DlfyHvhLL$ zgYvaUMt*&@kJ1v+p}Z6(pk_e1@@OTH9WdEzHWeQ`SXR+RX=5>;>XFGJ z`B)N_U+hH{-WPZRiB}wAc=VX@4ah_wCdL1Hx}V+#ml2Tda%+AD zU!bGIi!~NYasIK2o`*gCLZd{MEKqtIZrm(#E5RO3TK? z-vUzLQC($OekZ{LckblaJ!SO_UCc*{#Fwaf?*LM?P@putf#@P$^%n8PoDt_8-=hX5 zH&|UE`Is*~dmNDzB~(<(D1$n*;f0_PqymITL5}489dqi~Zml%|32Bg&R=#suOQ^eG zhUjb($O>0dz?@OK!6gqwJ*G@Qtcd!)OZF@D;rje?*xu4==f0RB1x+=(Fn<@_;PlJ2 z_(NrBQ@H_A9{YOH75AZJ8Qmh)_>H4+s0oIpaoAgA{|8A^WH;Vjy)=^g?S<$WNyPm? zd%JwE+O>J!9vaTdWHL7rG!C9z5!K6pf^rKb5?UV8L7&Iy8CY@@o*Ni7A$GBo{gnN4 zTzm%nn`4)p;XrnDAj_mROLPrGO}Zcsm5frbM`(hN*n$_+-pzSnp&CqOReAop$zq4J zPEDeWd*~}ZP)Unz9GovMqtQskLtF@QXksrpp7n)KQ~>WS?qtW-T^0Vfm9P4~SHz?5 z9pt&k`62!-yxtO;6KFZI`fafMz3GGia1uplJ-zvtYq^j5kPq|8ei1bYw|gBpo8yAr zWxm!(_KsSH3BEzCC&r!OfvbShmjCOFQI5&)KBY@pS{M`|HC=5M>}8d^+fg2RmQ z2PS)faB9+uk~uy2)lf;gUQyq3cfvDF3DOHBpzX}%Q|Xdr`}e2&z3GN-vZw);RgnCF z|Ano>Jon_=(18Pu-wb0r7HkT67zO)EWUPvKjHzWE$I;iyGTTtKGSId7nl>62n|$R%^L{a%3_PmoOl z{E=fp&TM1$SsL9`I*@ae%Hd9B^S`|rRd<832%V4QgwSs(IG*Chy=qP-Ws&X6MasaT zzp!Ry%7s|CMCDBd&ec)}E%YojF$08Bbk@Ev<;?t3f37x}J-cy&D&%`oF!f+X3$2BCM8(PEqz$aT3bRgT-Ze zu>$`Lp0us{Z~MuUo2Xib zso`m-Rb*?*XNuSc(25(6#zXdZv`D*;k#iK~;}10?Ypz4@KNiK)?g1 zWxCV4?-LP)pn()#ob0JRa=!mDyi0B&+b7{L8QCXi^c-~R2`fcOQO9ElJLj1>-=#U; z3Z2Bz3l}mU8rF0B`uWW*B;G((;dx)j_Pzf4bPpr4Z|4lO*Vnbeeu}^uN6~tUh5=d!Hr#-i1$`flU7!i|uMYMpBoQ>8Cr67-=NTCx zqBz5f@B>41aj|K4=rGDp37i&mo}>;qgKbuU!GvXOU;+(xnoS--!Gz}IZTLKuNeK-l z%mWjUz&;d9?GdJ5rE?4Li9}9{{q7~%8aI$BYLNr*NQazst19ZW?2(-&5Msaqc2Lh% zr}Wwz4gOA4k!Uptst|jVV~}WQfKh(TzCcRmQhKcsvI-zSeSQ6gQgjGZoR^nZvfgvz z2ySTL*X?lDUfKu|Y~pO$dC;5iU;p5=I|9RE>_=45mDYfz3&SL2jbrrKT^5+JU%q^K z5m@N&>f*EF*CR_Qe_74a6S*+69Ta937Osh_03&qU8kI%xZEGR$g%yy9W}uqi$bPHg z^ptSu4na_nm>uJ}OOh~+v8m6nVF#Jba9GA(EZxA4@51FeL@7c4`!GG-7lHAvTDYk> z1%NCyls2hWhmCwWAgUPbcvkS4~y%g@Ey*_vAOL){loz3cgGZmR%n$t=&k! z3QJ3q<~tYtE^v+FB5}<_u3Pa*1sM8o;n@eU$PHY}oStTi2T4gQD%KGOvy{ipl5Xx{ zQj+(hBf0^|i~}JJ!ZW^qb1KQcBe%%P#^wQ}4jPe8`|d4yPi-Z?-+$G=ruAfh`!+F= ztA81>1``}^(L+V9fMvLqSn{sS$jKP#9Z={EAx-38%Ew!@UsWx3w&&v7!pREQy`cp3`}OMAdQd zSy28&Nm9{hulw*^?`yPn!XUf4;*%}L+0>kx(7CbZZRK1&gQ735l{g}eV4)BuQDyi) z%$@X~o)ahye)mt$$+Z?1qn_*z&H0|n5GOW&g%H`ENNeBc#Rj~28{JMi{KZFC1O=t7 zZn2kVVAB4HiL+aL=LK5{B(Se*#qAO($9ZJL$Yd=@ zZshtCJ!k$#pZ*WEFYgLY{U^QuRd4o$nZ&XZ@hgx-L@l#KGCD8t(aQ=m!+J`n*0N|_ zvB*r=*)H(9bMXb0@>ZA`|j;QSOV$MJ`zA~yY`Ih22@mHiAE;~03?~9n(cNfJ&V>0njG@izI|=_4qM-z&<3>RY4av# zBhO(kn`0UiN5GLZYV-t$PnK*V$(GYUiAkgcV!l;UfEb9sPh-P}F=4CZZ4^y!zshCo@$(&AKUVGT zH)(e>=fGf?RZE5%Aj4H|HC-iuh76;s*966Iw+lt^O*|sS{7QX6mepy~SnLLHkMa$)jZEQ_|zkh^P zPUeq<`YqGWVym8GkQe;ZV8f|K;V)jk=NW$47;aJS{=!*?DC4NHeVs9S1w5~kI;MDD zTI$UQfAshHxbDbaDJyOgWvQfCRM~%gNZq-Zyt35ci>{73x8RgU&e7>-b3Iv_tXeJd zoWn}!gc3MwvmLyIQ{H-wN7O??`U4w&l zEFLZ?g+)XNa1PW*i9sPru61Xm_cT&BFuCxpkDRvu4IV0DqW>$V`lqw1031kc6_iGc zW*%(mZ!`JE!SCJtZ-t$YB&c7vUC<3OnWGl0^4HEdd3HwqsKTzlckfP;MPmkBAtowP zPYj@d1s`~)U)t}-h{6}-2d8C!N9=>z2^ZLh*5=d21;*Hi#$3*@Oo?@-3bCl31)|Y= zn0M;_(=OTn!s?P96zrROiYaxr7ni5^%3D0cd_S=>S_Hm|SCXFthgGs{hWgAuoDnzJ zOiT9pX@wJc4XjXIS=JI&;Un?V95{eweH1#liv}zGqBwI513fr>pGIy7wfnkX54_XE3|64S}>*o=jA{ZfYU zoIm+^v_=ACr4TwuyG>6BIswv3eb{w zz>Ql-FxEP}%TFLwiwihiskK{)ES=L$qLZ=LGUQO-<}eHN$fob>B{4qsr6Ia)@Js{B zA26z|JSt&3U!4pIvF89BF{|ZsThM=8SV*9RwDN2s3RLyI*oeX zGIX-bV7vA<)_##6qdZO1ObQCL9fdsS&b35DkI1JC~L>4r0x-n%z+uNFpSQiwaIGzdyM6imM*_)h{={sDIJ!vC-`X(SN!7 z=stVhdr6cjK}W(Paa2FTj*-O?xs?M#I5a}6#a|54;KHPM2h))G=_RZAPHjtv4Kh}0!6D2HRw^-yxp(AzSuWi zZ!;Y)a-fK5mnD5woNRJm(DJlX2vtu>LL3sZYt1fQ>cxi@r=nu;Koez~9nqL^0Xs4Z zB%X|d31LbjnG+Z0T79M|t=GXJb{*H2KJx@Ry_#;i{VN%*D!AK_gkf1bcJ!$1hoX0) z!#y@LG=i+(DTlPi7x$9sMHGcbR2Rqu;vYO{DrRR-=hul2QLUN`>$84kmjfe*?a3cK ztG~nhr`p@*8veE{mmvoiIgdG|Hn(aUhVB z<3@Wg?IHg$w9@8E=%zPPU}GKnAv~Q?)}Ffi-6dm$xSZDk_)On zY^B@pJ$ev=H6SeWt}3&oE`zFE0qc@~ZpaB8khNlY8yESL*zIeb%KvMF7)s zf`j2tGwlf1^yLW$Gu$;?dJi0US)(q2#%jI9VnE=|;ySc$-TLCv#~KpmM-0jC$#vI# zKO?S?-4)%vhZQ+epQcY%l7bl0sH72V?Z){RkQy}(kFk}~tX5FNVZDta~V z2TB5%&PSQ60HvTha-9};-3JST*@b;e7!k{~RY9*MNMO{aIIB;Q;E#DVr+X}qa@=58 zZ0rDikQTFv`zjJvHoj4^OsIIS*VXnvfwwc)#Fs z_%sVab7Tku&I#Gn0{K_=R&#{{}TT5Fh8GJ zoBaM{Tl=ron>c2}^iI)kS=)_T`P{8=W8dBeSSnNG_*&E|C4M^Tv(0$ zE@XWkLBoDlHEgRL6jwO+rv5VBm0?@SUw^Z{2Y#^L+4)M*`tk$!jcjDq_3LMa^Zh;4 zL57KmH2&a@Pm3x~=S?{_<04zZqtA?)yEyHRZ^VM1!+*$p&OX`o-V}NzBzSiA`jnq%xy8lD|DWkQozcygMwnJh0ICo$cz|Ps5C7u_4d#8O}ky_FL7euvMrW z>cns}<`v{!B{1O7 zjGnfsf1J)cq85p*Wi}${7*N3)$la&J$VVQrq?I^(`m|QGu{{{_drBJrZs>5_*}Fcg zAGk5#ZR+xbB#!CJUK{Me+A<1!hC}cN+Tos5tO$F3d@bg4#c1#|hP|C#6NItnQF=Jh z4NabBtZs?ESV^Y+Q8`$ga-XBSCpv^04&8@J5zIeL;B>b zp+a5Q@u6wX7_)UWZ&D#mTzb(>yU;eleMg^UF3iI|g{ z-g~mg)~(Aj_SAEe5FFl)7z1RdjJ2@lyazBsQjjv-&ibzStZBDt#eg1I(97qyPf!Jh z6~i4^Ft)aBNq7;VSO8V1nZo*+^XA#db;NS{f5cTfZpDgzlIb9*B|2)TRf?wDtby~6wId_*{$6_kp%|-!EXAb3 z{vu=eJZ!}q-qz;PN1g60c71HFn{@BqBzv{$(%|p4&+=Ami~*<$dv$u~^}OY5bP`Xv z8R79u^r`8#BK!OP?ep|<$WnA9Y26Is!;(0whobCNx?cBv1g|8?`$8Nefxb(keYHhs zT??q~V*@W(0TajtgIc^WRL?QxaiU+qkuECA5z!3-0s@F?@S&Y#M1x44IHbw!b0>y^ z$wJMYKh8&9N0b2;Si)Lc@#$qxryIw&bisGlzVo71oLD?Gny|t8&}diRRoTaVZ3x~c z36|3LTbI;8p?sjEeiysI{>X>W?^5%GL9{Haa?MaD=h!Iw{JitOv=HelUFi@X`?ExwLZ10_x^Ve1x0Jn?+?9xF zrj&=)RDk|TuCP0exZb9N)d9if(>pDWbDnFn)4yKe>y7;O;ntiQ=4U7!edM_X_6W*d za!FaOfyOq~^RMF)6`=wj$HC9Fr(J%WMH`p71T;ArN?d@_(ewwpr=y~NAMVM$;NUws zJ?OmY;6SZK_13q;P!&LX3QF3u8(B zjpTs2G*R^!6X@4in|`w%y40@UJY}9?(4M(A`Le;#!eYOv zkE+$b1n&>$q>>`=MEk3eZAm>mF=5OqY}h|6>s!rCLntsX0)TGbMz6H0{H+lb?b}-Q z3ag(AzMqOYTDNhww&yaQO-Wg0GQ3N}dmwv)709R~v_ph37#Sav9v$2pcF%|I3_ss9 zxEr50M_gK6liW^^Dbt-%6c#{fVtZ* zZh1(1wyXryMoq*=lF=eWwdA!ts)&a!z6ezVVs`$Kwy1`fk(of#%0=BFW)0H9L^!(& zI(_-H!J+?4jw!0vz4!#5Kxk*h_qzQiavwZb6Ie(e9;|(Jw^dj;O^7X)nui$~2c%%+8Cc7mjqN=$zev06 zWNg?suiQEJbl{B;=x80a>%C-P!c(A@Lle6#mNCdb65+(N0n+USR!l<9ws0yVv!%}h z=WmLHs{4W*b0CHh7%E1N-nKBmn{@Z?LPT)M=dA3ZaqzfiLFH$GF$^lu`Rp8JYSX*y ztbNq=!s3JFwfzX5%_3eaYq#b&gwP$Bn|^M=^`CQ8Sr$44K5D~n3L{kb>eX#Fgf**H zE%F)pY5O#rx)QF_y_3@d?^C@C|3#{v6(#6WEbe7hR4`PD7l&N$?pmumnc@#fwq2)A zkLBc&3lLTdE8^^T|Lrb{xTm0+gZkF-@LK{XZtzHd7ZK?jWO=ToMU!h-q5orgr_*~k zU`s?r1}*jVtNo#&S-=I*WOrm7U)(LtkviC3x)Yojj5X2ZW(M=MGUx_pag&&UEsfOdr z7f~uu&!Zh58!$?_NcW5*9S31R84gZoxVbms5jb(*MQ*ZZGm%XtGKZ6mdctbuzQ-}+ zsCccDT4^>cMN6=aQY)FgwgRnp(_637cfZQPXyaWFX zM!|@g$4&>G0Rcpc$7}M~BC5E-jv}F3w3sQ^_Y66D->A_H666*)d#iiKbx!xmnDhKA z>o`}CL(N`O95&Q+zobj<-TN2ULJ@UlBI1CCu|2Jp%1@D%Iu@HCrKG=M$IxHfwiQ7T zN(-t9H(y^<+|}|Pal2q0eM<~LGu6?iC(M2flHKYs7E})moO8uZ+U@c8klf>6&zRTH zvp-0<;C?^nH6WDQ%!{-7Dyt_IZumJ(wx|L27v^W-av9stFZlH7iTiqLCq%Rk=B*>s z{+lFLCKMgcvS1UJe{d&umRr;7a;}e-Z^;<93YyXPT9cPIPX2Wd<0HWf6<`kBl(u zKrx{zv=+v@Y3#ZODMNgrgu~!3=S3t(3zGrdG@xK?}HX7n{kTt%N= zv?uXD(P^KNgc*x67}iSkjt5CdMy+~o(tcITIXi~}5+37S8g|O(?Xe2ES`_+{wnrn6 z6FoD;Cw1-<>%rXp6Ycj+4w0CB5Op8eX&~uGdejwUAl7%ow%?}BzG%a%U!?eTNqh4dYiA7bHMKUH+;oEBoru+% zM$Q!`W_oWwpE#g2dy@a3{??NZ9lv9_{lU@W+nO{uA97V?#NF<5;}7R|iYxqXby%s< z?1tOQYAQZe4C0_qs!STV=FNj^Z;sSxbwqRBlucKIv&o%)CVCtR{VKRP5Hx&u#462F zi&q^7*2ab*>y^B{mq;U5d!`j!ii+}}HIk2d;6OVli{}jb!Ts=VWsdHzyVR`R+BF_m z9p=(A)Y|RZ9*qQ@gJXIoKl6$kO3kp9O(ptP=dH_zSz1~WvYlL zQw=8hXf=)oW$0eq zqB$D}^g2ufyG1OKu6gCx%d&;bUQ8L1L|<|IhH|6En=6 zohHGxipSQ*rYVvosqjYUtnn7d6HkjMhoYmSnKXK*Y5%|7fZ#qmlEWY#TN$kj;`vRz zFLd7okfL|?I!u=xprs{~7qwTt94{aeCTo_FXlWLp!xae+h@McH@4bgWRw%tv|`NS+4Jc*|^Tx(*!JhE7|b{hdFr zwq-C1>3)2?f-d3sM{LH-zOAAB+E*Vj$D`*b|<#n6^}^%Scs?Ml@p zAmJCZm=X-j0Fv}$Umo1Q&jJ4E!mU=stfM5y-L5S3DO&RQQB_aad%hIGs%F@md|`z}n*`M+P3tyA^ygY~ zs5I%F#$fw4oU5X}&)Zbznx5Zn?8G2ua&2St<>-|2i)Fas_!5HGaa;TLn7u~LL*kxz zNGT<+I+baNSAcKjpbyK@*pg!U_R^dYko2OkBc$IedfXUU;OCYva~x5B#;kO=^0}4K z{ZE@w2HgLhxhI6~D&i5$jd=E*%lf)q&!$j7SFQ{Z(?gp!T`R`}P%o{okXMQ8!`&0E zZQZ)7^I(e6Zp4$#lHjvi_dU19=xa9-PC*qvM?z^b>%*>sVR9f};EeU{mph6WMb+Pa z`|Q`8dv3aI+cs4q_E1BskikNE%3p2KZs<9_i=QgTF8S({LWW30f>5LNHoJ9_*81K3 zaYd^Ktxjt}f*NL-=YE~cPwn;3XcL?mn&}Gz^{bSc4VsIuVC$B9o;QkGv{kac}~ zX|9V_^Z{BQN*J@6_3J(k#O}K`?6VO%SMLIo=$(C0XWGYBw3BEGfP<$`pML0SLTmTF z{Fi(N%V*vSPLp{!0|tzxeFls$nl}$HNrBcx&uFju9zVOK^o!EUhzxe`-qyH&BrJ53 zO68AbWU2me;U6jEj7+zrtT#||8*{u-R7GXOl=S&PK)+rAcXqVSPkB|02qIo4V|?0J zVxgEC-EI=01p)6mfj<@qeZZw)u7f+w=%*Y)s4p%W7gPQF*4w=rTyw{D58Cp$EPr6A9it&1uS=2V&fX>02REF-}nz}5{H%tu!WJh#?5pSR@h zb6&cQPQsaq-wmWkbSZdX!#(6AdJQ=)FWdXh<<-LJ>VUI_0>;>CKNm zSRVTw}n zt=t$V8T-gVg0<3eUyjARc>>EyA}~P!!luGJ8=07FWd~i%Xhueuem7s%A;1z6C&3AGO=Vpi!Moxa$q~4R*Z~@D zF)0jrTK;@B z%PCnK=BG3)3j4;2nq0PjFFV|purRx!`ufV6XM)$XE$YQzFiE^_*SXM8M$Z-c&9Pp# zY&3-)7nG0nytp{K@8e*Tipl_eoRt*cp&Iiihu{hzlo}H2%^UL}g9l5EO+AA{*#n6M z5ysB9d$tyxB0b#9`i)xLW{MSqitjCOuCMXTZRU&_il1Xsfetvy?(1rqtBYTf`PWhC zQ#P~Mq%dYjl%)wSgY91ce>NJpob58JjlHGBx$}ypWKr=Vnz>0hO`(6k*>JH%iGX2U znx)eG>;0o*N3e=oz{`_vRJPw>zWl(eIuwVv=#Bw2po}(~ySPDe?s6aZQi7b9)zdN9 z6Q`e^Z$C#fK)!i=!=Il;a{wFbSS~m%cX11rNn4(WF`bVFDaP(UX0P9#a%mid2|Vj) z;}~U4U4%9iv_1Ru*~+#NeM|F%!SR^%=-t#ZRFYDH@cA7*#wlwm6jYG2=FTP1?97O^ zW1YK^Q+Kd@_<$c5Q2p(tj$HV$+nb5EmH%1tvqe|X!v{Vod( zPV8(^!o8XEXwQ$2Vw#d)VxSxO*WTbnW_Z_EP=lJT=_W5^Ew1V9Z-|Y8Qsq5{hPp=y4&CRViWSZy{ z++#y0m$A+4ii$q!PPX^C-nF_oJZkXbZWG#Eayf_z$T)n$T_x!)FV#hw9k`%)-p-xUICLodW# z_h0*LznZlcTchk`5>g4cMC;=Xu`+!Q_lkz1zPQ!)`Wi|*Rml0vLNqRlX&1g)KNRje0;nbC!Hu?GNksy@5o5~+^T!0#?qUzT^*F0m7cR}a92cR3=$wqi|d07pQChdh2FJQ)E& zk#Ys@=P@3M+&isoPT~gfn0ioaOh|OGUJtDV5u-n)>N!30Z`}!Oc)wA>8cv$Dh>lLA z^^@OOSJx0%DQWx_>UNS6c`#|MMx{$4oa(MHtNJum$U`0JS4i#C&XS^{fZw)Pw!qIa}gz~je{mn!IG{O5h8Q^>+Uq-yCE zsN=X}m0BV(##Nh8gPrN5R%E*0b=oqAaF@5m#RFK;SFjx)zzIR>^3C${a*_RH7V8xc z<9dp_;Wt9&Q!_J{q8RF$5l zq%HLvpYlR)Vp+G^on?T-zP(R9xs%jm&qT7kcQJBON#gBSGCRUH@9z$Le&c8#outMJ zDy(NO^rDT~S4Gnfd&=6K{sodEWK!o-v5F0aeolZVWGMSHGWQpl4UImRp(OF* z&<^6FqFRrb9Yt2443NU%LfD%J7L67ax(nlFP_+A`wR}EGjJGp##t#$~T1d!*7oTY> zZfsfs_$PzV*+axJPHQ@R`0!%oD03cYld5n4pyg|B_6lrcvIfNgl!oZi=FFLsoRTt@ zCtt7b+nQz8*7Ey`Hu``1Ui|j$Lh5o5zXbLv2}^#gS931C+JE)4VxC+dLTE*Ik75U7 z1p?(|mp8%!qFH_Z>ebeujn$3OezZraaO~{aovif49b*g(8nKoZZ2Qy7>Jp6vvXDD>!ug-I6xeY|gv}hH%f(9YTMQA+bJ^(z4jQxs;<)tN zmhQS58vgyko9-eZmb|Wcv78Tu@80D+dScBH>VOjU4yuUUS9qGu-lYsbW~nh)bFOh-K0fddu259FsU7~9 z>UCn!Yu<&oc0%;V4Lej)q#%>wNs4t^jfYVvV^{&F67@g2H5r0nn|0|I!0cP_u^$`xl_xp6ih9C%b#_K4%|TSX)a}sqB|)-h+DJe^F^UtHBQ7%yh0O?pLx@;+#gVTj}@Z<1mP*CTuF~xpB6X zO0-CFuQO(G^;?i*ewY}>5%Sb5#Mu&6Zk zCYrhWwfPJqHrhu|_X71(af@+DL_rB*LqkEOsBmMuH9I$JGh2nUg2gR7cCU~RapG}D z6tYM#&Xj7cWBH27AJs*D?}9Y-7WQa{S*7nv?&>UEqi>uc)9k zRKctXb?;m7W>^afsmW2gBR4JMnu#q>Ap?`151Y2AE~8M5no{uEl|bttPFH(K!IAtc%}~^OlQpE!e7tGA0L<0_)vZ4kQMf9PP=}r%+PyO?OnWiuE$)o^am1$L1 zuU?f4;4!|-4WT=6q|LQnz;Yc&GzzDNiJrz-=f3*NMFYimr$7$N-q$3D89VnC7vWbg zVf~8xJ3M?3_7g3|dn^K8+9FEpkhn_3jDR{1B2%9|dq{wERVFgWC5nSwa~mcTKluFV z(?OuF)Qk)(#5_$_=^Ge4i0o}Qt^<2Ks^b3U?KCu!6c<5mL;psw@wrn&w^(lalrDQd zyp|Y))b#W}xg6)n!td|n?QNxyO&*~{!Gm|j#ZRTqNL0g{m8C>35yrsynbMQzwC>P>eHK7t;YDzluu!$(p&8; zy~2YI{#`(V+1{9E{M-&OS18MRyakh|_EeAJO0S#4>H)oKXgOs?+pxiG(?9cm)-B0J z$SXsKz~9FJp2J|KILkawNiV)WWWa!%JfiSi9x~&Cs^{Afbge?#NGQdcHJ8WEO}UDk zugoC9V50r0inqPCs!ScY_WdmPODM^pg|N9ak$8}sUw8}@B5G17*jXeW^CLMXH(=E$ zy$=ps)dl_0tJa+r5?;bpPuUo<1`_v@&N>rydb?fEp3S3izrjJgCS~b-r=F#Ncc-P^tv502E!P?I=iSGT6A9&Bdx(c|%cCJ7?l@_`?)kXh#if8nOue}w3yIhkQF~dt7L7dgWlEP8En7+ew%Pf4q&+Oy z$gmq3HU_|!dUC4qHCl@Z1l(L&8{-4KE(DHc;s*kwdszPx2tpKUBdXs?SrM1Oy{+0>$Sd&_2`BE#e=Ca>RJRZ{jP z-=(zbVPQHPYLFjKpfTAYIdLaW2tGmbq_FL7jm{bpqR^4-fJ8x0u;EiXNY*HU4V*~_ zY@9f_Fw<|Qu8ymE*;L{t2rc5GEIiwFR7)535_#7=k1~@m5wWvO(q$m(gcmJ_D+FAC z0c=QtCgB2N@DX^ys*|&`vji(izB0>m1HVNqcMSC4dCZA)c650}t}%q`Qxx5kE8q&z zH#C8Zl)zR{|OsH#V+^W=Dl zy_ZvL^*i%c$3vG1g6X`QmbE>P8Ja|6r=AKoYUJmVoIttPQ@e}n4*)V{3j~%WhZfk| z>rN_Xr){dM9m;$)n5tu?Gn9*l2>eB*JHFw+psk<( literal 45708 zcmbrmd0dZQ+ctV7L((it^I!;R(xh3XB9#m!X_C^YL4~416B?0_q=C>RjY>0xOr<%g zR5WO!QPVya_w&5(v-i6{`;Yy6e$RdXZtDBJuC>l}p2v9{$GI-=)7NIC=c1=5igCA& z=6;G=qD@gW^>oYeClV&1!T7Ho9=i^E9B{GqIBn%_L+M$0xH`FbIN2W)@Un3~W$$uw zyO^Aqq^N+MhllGa1#xlb|Na3n7xxq5ZS-{HLsq)#n4F?0W-IdlG>=qM>?ta1X1C_f zgJW_o9%bLQ0GD_x%tZ*~~hMHvjAiko!} zm0Ye)Cm6u8gktF{r2hWvdG)Rk7W`X}7<)MREA`?3pMOv*$jr>_$}{!J&u8*yt>@z- zwav{TnI(VdS8v?we^@l=^5vDH+qQN6_z}!ueuk#Dwzg0#iluJZy1x(9IW)w$X;aE< z#%V9FlhqEb)x9+_W8qGFGd=%n#m|;(PcyBe=jUdmgoK3+O---XCTUx@7kbh$GyiCkla>9}oMnBu z$ZMN^yh?iIo|S@TtbxMnIz_oexOR2(GBATFz4~h{t>W{CTH56CWQRE#cHGA;M07ZnEVB9UYxNw%~4DoTSTF zC3X&uR3+QN#&kR@_NTz*yp4?wd9H*6!9BXVhZ-9i8^3i&sA|QloVT>GS+;%q_JDu@ z%kHW$J3G65ZB2J%9LwBz8q-WJefxIo#h-!H(a}-=goK2*kB(Ve3##0H*-!-PM(n2EPydaS#vr)R0PwYA4!BO{)rXyoIYS64a4 z-fjP=ww|s0^=mq;?wLRRDN&J;t{vShiEBbx61TS(ojyS|54RVYj>}UFlxeB2GHn^w z!K9<){OzP9;hHyZtaj|!!9WduEmriw!>*&g4-bF%Gu~h6|e@$e#b>vE2Wx zkI!~LKR;Vb%Q8v=|6*iek?GDuOG|5OXZIR=Bh&cIx|(UZe;B{2xy^G;AuC^BUj}Ns zgoLHJx%s=MrWn=6&hkqp&d%v?t&K9x-&0R#e*efm*5b9CD8^^C`U!4Gra_Z0+sesH>~r zx_WgLZJB(yb-drg^H%E{s__=*#+GfDmOhT%6jf5XYhr5p`{Pj0o9hk{j&^pI1}Zb_ zEmBfbw+^=Di?6<{Irp)xSk~q1yUgtDY-gVey;cea@_8dsj z@2z~-kZO!=7c)DC@IioN9SeMGeWUC3r%#bM%~`8B#6}07J$rWRu&6E7jBkGYOWua1 zCCfUv>PQeHcNO(?q+~YSVq;8A%l(~~IH^g5t z`zW+4n2zV4Z2#7@x$-Tr-eCOZC%9+^NiC`t`Im54AyLKxpOLDuUq?N__ zi8YJGwxL-Te)G{LCMIMTktPoDrOeGBBOb>Yw${_r3ze+m;^uA%UALW+YAp5h%WblG zb8Q0yb;{M%(r;lV<*mZwgHMb)uxoFyc}z@94uW<~8D{Uw%KS}^9C(9MRoOUH%$G*S2&D=gUpSQ#7 zCnsguRT1bgC?8 b{yrQP;k1DfezJhvBo;;ael#`QlskNr*L>jXUwHoa#r(BH{K;6Tv)M!_2S3K;2UI$*UTP#Yy`!IdpYkz~SYqSC7;AWgOVQeECv; zGbJTOV@q_!moGf1SUIiU;^N{A>!dWD&mvtfOm{QW1y;!0J@mL3tAr)8YVnz-VK&*sw1 z%n{lqMS1zM^75U|XV0G9(c%B|$B&e;tc8UI%KzujpGLfzj?A60**2AG zBEV8sUA=uDuy9!0(q+rK z>K`9)sD6@}85|YGS{1>cIp!9cHQPw5H0@dQ_U-<&?Ki^aXqPj5n!e?;*&-8Zgn>`# zaJTHl!~}C1dBd)dNT#1wN#awcep9h6^Bdz2>MNa{vNO;PPrsa-D}!j8z~Hb!>9&uiMbv+-!+rx24Va;6aMo zf8YTBq{lvER3EKN<AE-TsT#*((DlzkdBHAuX-x>#Ks7ZOz%JcO%96W4vld+E8?+v6j|y zcKX1^0iom#C>wRV!gv+eQlqapEFPs>hmO@mx8+8kh(Es;a7;u@tV&nBc)Uzbvx-w? zrf2oSprmz`f4BTymGeisU-}0#^9N}PT>$715)x`nv$dg9aP6d!$|CRGV=^{2X5g0H z>uhQh*_vftecs?eZct{XxTX;6lho8P;G`$%rM~kDz&;H(V{hNK0#@56!ahCs`(s<$ zz1elDzDf*?jBX2aGr{*Z(hdjY1I^xuj`Tagw!--bM%<{4_s&wwxDrDQ2 z!8fX^v_#Eca6W$g_-S5V)z-`xbq}j$WM!*x{H+`uuAq>R=n094FmBP>c)bgGvZCw) z?HEAeXif2p7t^DXr@rrvPEAGzp4s}-{QN|NtTM0sxoOve!a`tUT9=ZUjsA+lmUl85 zeSW>!BniNQl|T0F<>k4@@(+T7mh&ikv!D9$!P$^;$puRHrYH?1dL&23=D3znRHy3V zOhoI?AA11aR6@hT$hTU~|Cw+oR-Bocnf{sm@hKA6+J@xI$Vm9=L)7Vu--?YDzxVc9 z4YlTOsZm*&wsz$SUbbqjCF+nv_FaooUtWISX+ggY0)j`Hm_9Eoi!=^CIleIy)_WWgRG|`bUP$cLhzae}7^A_A+Di;_yA4U*V6HNU^k3?&}$iUdgxdACVWrCcccENGK_}r#$v@^_*|}moJA2 z{yZ#NGQ!Yent!tQhyAmhoa0?xU1m#pnh#Q?fY8oh9 z#i@Ayk94fkS*cSV9z7=_eNo_0!tYcKwQb>%^!lOaF*`Zrgur zaInVwOxH56KjZdc5vds&wy0?~G&D5t>g#XZj*ZPmb$tembP?h8#JOW_SKqf+SLX)5 zW}opp)zl;@%^f9&Un^c^bGowVPF zhTa3X-NLTN%DM1u6%{>R6C;i)4|w-vEDLzbY3$ppch%Lwx9{AMTD59bW6^1s$)4*5 z&x(Dfp6JQ8wz|G-j=QW`y13vY>)xYdd*VbjKJPFPQ!~!lhsmL~6gNL){hu>G;!h}h zk6WV2c-Pu`XJ&Rb+y3*zERA5M$M352LuB1TbIW^IqOtWjiHM$_pXwL`-mw1F^QKTN zIFjiqI>rExC8#=p&2|nB7jgKnUcH*3E0kQB804f**RMw2$2ZUEL0X!q5;hH;PZ`k9 z&hn5*ZB1eA0G7H7Xk51U|_KKp3=I^;z=$iXJ=0{Awga~G@)q2Ba8NZ#B;Hzzcyp=FjLb`-X{_n z6SG2Bw`TT(jxOob=PNih4jo#9qAvZxHB+a{b8n6>4|sL>xbUyfeSUw=#@d?*2Xz4>Kn6ej5m3;I6HqLF?}B9vKSC-xP zeWv8sa+Y~(Hyi+R3JDF>Ja$YdKwXd$P*mhDFE5vnl|=#C>Ec9rBct+SJYXjgR%41(1Y{82>Fm04O@(OlMCa#{1@ojY8pC}Y+)xvXseJWn8HJs5{n;#eD>haG>y0kGxw)X-VIv-EeAg-M za1owdb)R^nB4aBnL9{jCjOcJcu*^7Gc3<++(il#Byw`i8fAe@6iiwFm&C1$|I73K9-Mt%{m?*S##R@sc3}NsI=sNIENtCHj)FWwUIT;y2 z8JRUgQFK`6PI>)UgUcjxjHTV{>HjFXCwsVL=Yc;$u<=Opl)6+_I^G%4oNH9I<;{a`c z)oc?L6=1vE+10gkyE`7eRuQlLNH|cqQ$O45$PuoSUkYeC^%wWyz+SzRD`P0a4kQy(^x+kAeLDqb@pXU~mvLlK-`9Y}9G4g4njS z_WfaDVe)A6{9eUgxIlB%%&tC_Ix{-Z~9AZh@pmr_07zYA^Hu*}QL3lZz; z>G@#mSCV8lE2KZ%6}OpZbMKEI!kahm-uk}vebmjHL4!@1zb#fRT_V3@2ad6!upHUy zon1Ps{zqESUlUGE*g69+BhY>)DGs@vdwh@{iRtazJ-*M;duc`RwLeUm#W#5 zj)UyeJM8m|7d&HQV~-v`rlHRG_#_!;nh&XXS8{3`IIvpP`w;kF&>t%(OMCm3S%(>- zbY2=LkVO228NXkrC4a4^HzX{K1QYgC^7w0(xTzLt9EEM$!dx2jaCiv?X>d|yFtatlO{C8HIjO9}q3Hb?4m{j@ zwr5`T34(q7rcD|fuM>n1+=UX{WUnW6n*k6LKaV=o8xtHDh=uKcDuwg==SoG?+FToX zg>%yq_q+!-f0=dO_knL72xwrJ8hQ@E%q9Gq00AuQ?8zKS(w?54XnjRQM9R!D;8dJZnv(9j@%icXPG z--$^{QP;004gSLB^lhfsR=5!#-|?_5mxG&|8+na=czD<{wCip3R(YREQ9qm55>m-r zT_vhQ*IiOM|ECPK22Bc5al!f7ytjEc=L%1*;dGWH_@wqnrlvY(=bdtMt99S7?S!SJ zKgdqRG3w~Le3uTby?Ys!(9#O7z3ORpkKi*zFIptDpvdS3HuAYop3qB5N}|LB;pMMi zU#%yP5HjPOKS&P{BPzgL6ZlS4q)Yg;^z>~qNP9Y)b~3gcNe)YGpdKwrR3f6NYcV$L1J5Cf+KwRY`Vg|u*#8`JwiKbelby|HucrG`;$6qS?IJh{Jq{-*YNNEn7tmkw9`f;fYdaJYLU-%-?XmZ*Z=WWQ zq@|^GLZvZG(z5C-yMRayMz*_i=Z>YVZIz3$ksB$+G;2wI>u8*itash_xu*l-ne}J}BiZ7#&GZ5}p@w|2vOz)T3&2V0 zm*S!ES5GKITXQ9Pt1)7wyrE$YIYHZty!}MwJ7+o=tJ9+D156)y>IH$*(@kX^Et&p8^DtyQqOefC1y(}r?l)l ze96E=D_5?3l%Bo}H3O*aMRD=OWR9h-h+o)^j~kH9(O^*&pa3cEnnK0%Ul9?;YF#>4 z$2aiM@MfBrBB=3WSUa_9=SJUK;HgQG*tKidmhv~{R!5#I@UByO^98nq9Fl9QA3zIFs3 zH9Q^xLKaja4G_wKaVu7w=0t+txKL|d@<=Gex|Fn7-N}WR;MV&#TE8i+od>T{;M%oh zF&9)02%4A2B+&?a9$rU67EdFvWeYtTlE|9g*B?GK*0?6$0CY#=h~Ovvn9AI!MzJ%U z*w(EYfXM>`1CY7Ar-m8D#WVH&1};kmLGMC2XJ=;z!2LZ4fv9t2q{HfDIZ4_i#rFSx zUEA*9=r{&C#|rU4kzB(g?|Q#d+Oe&zO&GY~;lqat8ENNHW)hx0Wx+b4JFyHEl#;r( zxGoh8B{VqIQwk`rHwkj&-N3Wiv&YsM=@Q%`;GtfVO>QdK8TK`6%G=xHEN%j&%*~yP zUiW-r!bb;>5rW6?%yC^xp}5?VaWo3HZ_If7m|!)KX2Zkb{{H^@fHOSEX^$Q~>MWH- z=Yc$fmBQjpHyz#LZ6Ll2P&L>?;-4+Bbs zS(oV8@VW9A z6pgO;avjOGL;8=95pbAub94JTN>xjh(6?%7Y2oGKhQF5vu(*qpr&xz3Lffe^JR(B- z{CQ>6<%kGo@5!ILQ^M}-1aSyK9Bd^;q242BH~=kp93lI^J(z&-$-&Iha8N)r)Wu7e z{wDXCnf?N)_h7;ZI~pI)Ou4$cszJNjyLYe0(c8Y)akS8Yp*z^Udu8d-me)i5{g3lg zj8+OJ%Mhe=f?0b{>HH-3_YZMxIJ+)iUt9)M1!P{oeknSZLa|hhN12&Fei-c6m+l`K ziGZR?6`bO_-Yt)GhEjQo;Oz`Hf-{FMuio?x&)3+-8i?Eh$z#lA6nbxPX{oAN)z>A5 z4jlq48XX_+F-ZA}G>3i;OM`zvs5ZJqgCuY#FS-@o#)n7cKjQtuLD^k$Y5ZDrS{%S4 zX|PUuzX<#5w{I(ejc7=XMU$H(H?B>pDf`;Bzbs6T96{3t9^v!R0|Ajtm&nn}MN+wO zBcT6tim^M9tH!manek~7ADCjin@{7wTe+W73U|hQ^ua&rA1CAeG zfAr{4wA&!8!y(+=FDWr^df+J!RD@OS`<%i|pc44cg@><}Tu|G<8GZBSlYmXYuXbqe z*AxyQZxLM^DpW95Ajf{S`Ptqz#k|qa-MW{9I1z}cu{wS{03YDcc#;;-!_LX6$B!9p zlH1(08z&c61qvYAM;u-6*;6~|`j7PG+*xZVUzmhW3`AlTybO$Li?O+>si~~|{08s; z(j`HlZ-4^FY7jir3&|OCvp71+2wKgqzQr$XX+uAMhJqP921dP5Oi%cqdO&d0zqV#n zZ1KF4Zm1;K^qoRc0;;O~Qny1bIQaOOptkWEY18$;MaaS7;R&d4OV*jI!!PIf8@9jO zvFVton3xutV`Rm_!BrsT6y`ry4WH~0+q}6P`mmLq9bM;#*I9;kXqLmmy>f|`zhpD}5LqW8E6^hDyRxOMARCuG0!xtZS+UI)K_ zwH+dmRZ(%UFik&9G6?i7-uM2KCs!;ooO--~(cKv(z6qwSJjecXqtv~L9^A~4+0?$AQ@_&}-w51=nD+T`q{`npv}E z=gg0=Zh54nvYHy*>NRZY9+&?ujwxdOOm~$16=HFQC|gLI6CVh_NY2I|3C+S6yfD-^UIK19B5c@8Gb7ns;dw6BQ-o zFJQAE3mxzyjC0J`kvmunk}MIbs5Emu$)S?#*RL-_EF*zSNJxBgdk_$*P2dk$(|AO{ zi<$4SL}a8Pl!82j-e8ADrj;u%{Dn%w=};T4An`7tz%rvW2&}!@)!C_*qQA1=|KGuz zoI=Sv23&qFdlD}>HjtV*(EFIVit~RbMt>Rai3(3DDdCfI`5J)t)YD@u{Q|2MQa!pN zLUURL0;=8zRAmqMNCKjTWcZ5TZ;poIlsn1eHS+a+N3mJQzHF}_`_0XnKvVw%iBuUx zBr8Y9Rh{}Gjnz-}>1I^_W!GVeD$vYKJ-*`p8u&5!+E>4ZvZ5rZv-U9!cB8`<2?L3-gV5VeB?YFypuaM z(h*TsrbaRTotg@U*Z{s4Jju*cO$Aw$ zr8C~%$6iMYa!MXs>OC>A3@pp*UiWD8^?~XMq(%P|tQCIbf>!vrfuSKQ%F1VeE#v_f ze*OsL97^ruJ>{Ry%$N{*CX%SvUbV8deb~KaWk1)y zYF=j$DhWdzBTxvzk^r(_GfAHE=usZOB6>J;z|+Lu$t}8t40Zg(3HWwE zFN_hQ95lbeyp~SQzskUC08p>5r}qlpF`6&-jZDnUDYIuhchUzgLFpv4qKF7M3W}Qk z^QVHO^#<0utpEIKQ}O%ORcl4daO@{0SD!)h0+@ef5=&87^p4lOY-|;aUm(`+#Fs;q zXlc0(7O_X=KMRzObg>u62Ce`l0uc@rNF|Dk9WC(x zTOs~&<|L$CTry+kJp=ma!NHDW#=gY`ltf1uwX(I@;dI+&mmev%~}~ICru4L|u)V zb}PR%11c061k$d6M~{o+h4*9H8F|_M3SFJezvLaa%9?@nnC3o1GGIUuUghuudmtg9 z`0D8DcE{#AJ2`>k(R$yiiKC~6(}t+U$)NeNQVXT%$B7;dsa<`oxsKEDQ@g>Zyzo`c*4Eb5!$T8PEu12f9(|wmeaum%e@x*= ztWhyYXo1aF&CXNteP$a~pCDsfeH;51jIGUj=?49xK*)^A*lTglYNm@^rwt;O8 z+9#4+>Xv-K{acJE5{;wV${k9`}`-uUf*O;1)`i4eZq+oCz@MY5aiai0iDxj@X?60#vQR`%#Rj$1}8{JX}!7qRH){He=xp!!VVbgW=b)=`am^X5wPn^)c zuE$mkBOscsq>Ld@3SgercBB>i&U?eebXeTS2ie7r);g4VCnTC(Xcnk4m`_Fw;CiXI zxg21p`4`tjq4#+gE5}R-0?P8;ElB{dh_AAaV4k>%%=4G|oO{# zfoU%N0#t^siuC_j5p+X6i#fXTXyB0wamsPNxboL2c-T3f&di6cLG4ubq+G|jOW+U& z3Z=uqaFaR&d_q0V&epKAd-QcuFDiT&Fc~ln&@6~W&?Hb3Ibi_e4bjT1m zm?eGzk$KQ|LqOpHIajo`#aj5nrUR%;~V&&?MWep7z zC6x~$CZ6mco`8SEJrnif#YCfN^NX|F3Czq$Pv7;I^xL-fDk~?aPKs(sBT#rZ%-Me> zh<{W_J*{?Ym{DaEnkkr3++a(D)T*>Uffrq~yDD$G$t*+{aHY$w8ITlkwg;tuN3qTy zFf=UBKL!W#Xw7?6P&Cm0I_|=Y-%I@R2gH%`Wcd8%#>Ccz(w&Wj)lL8y`1(dYg1HF! z=hj^%sH$qw^x$!qJQFf<1-#QZGC@H3BpER42IR<~v6-HowT4T9v=PXPR#-FqKr}N8 zpcP1>CaVX@sg4b=YXK!|ff0Uolv;4h^ec1i5Mr+6tPJN>eCpI91%3f~B45Y3-;5je z?j#c^|K_;P?%h!Om#MvPO`aFwer^Vs=YQk&?aNrtihn#fK)1)+n4FdA_l#!IMN7b_ z`?9J^NL2J^6#Y;LHu(G7*RRXQ`#;0{K}V4(hb;Tg8;cov|1snBX?h-82){i*3H+Cv zLo@74Pn6yK{sq8eR#JZ9eP-d|IoEON-Q6A8GDS^}?BP1PSk%s;{En)Jda?giT_L-Z z>Iw;22+mt%(dBq>O*sb#2Q-Yp9Miop9smdiFl`m)y=_!!h^CpSCRmX85xs)b!}Pu9 z)_QDK6+F$O>dTXD`(7OtA5>=6$#w-%v%-Vos>`-lMS{shRo|B_Tl*Z0CYX5qgB+52 zbeT8wS^V1pW>NdK=|M^8lsmzUYGRZC>g3}U=5JLU8`EwMn~xaE3#A|0zyIRd-`~5y zI4q;kcyW<|pte*G_%}%=L~(-g59z!THoJ|#gcZvFinV@zW#l|#Jt+#BpRWT2=w)_u z_xLRYC#qlJ>5*`B7M^bsyy&%i3?-^Y8;mWi`s*I>J~L^<_5zAQp=@xoNlx7=F0O54 zv=*j^Lm!UGD<}ZLK(97TYS5h-92}H{1M*LTaPk`XwBZO1C@7GJ=tGL7fzsGgLV97D z2MU#Bkp25hWQx?>X)uP}*`g;4;i44NK_Pk2DUQgjb60akrV&>+3g?Y%*G zob?+wKKbHwgT;(^f>Kg8VGQOXtFYF6Bu=QF!NtYOrFvdovM`iwmWZ;oZ?MnB!dUgZ zj+*WY=9hq5+PVHg6a2vTj+^xFmVwYy)XdIPPw0%zy`X5!41S^8R0pb(kFS+&5G9k_!wUe;x4SKGFW*1dd=>gSTx0v7oY3e2f@;5 zasyV6q66*5SOg6f$|XbZGd0`*1Ol1@2!w=^PpBlOjWobk-M_EayTt(o9^M``(Y1q5 z)$7cWKkv-i-xTZJlx)Lk-TeI2rkkQB6>^R!*iC1?UEV}+7HS&c(YZxd8OPP5^xQ?BH>NE-a9(C`zENdsjvrAiB&Zi9xJEC%ssM0^~`7j}c z17elo1uj9jD{P5Kan;P(wtf3*H0M8NBh0<~H;^GOmlKU8o`Y*4T@oB`XBPrE3qg3G zP*k1o5qJR<18{qo>sm_T)w)r>(ex3i@zN#Quplp?Bk_Vo!ZoN$Fg4-iqB{*hD|Fn( zMrQWQ_CVn}J`g8Q9U8WrgX@;^$4Tf07{N3Eg7;2CeLc+e_gu~l!N-cuwEu}wT>@ZB z7xa|5CJY+Q&F!u?ZW;sMbphprSaczczi+))!D??nwmZ!vn^*=R7o8ig)jpi(v<}Ob zxPc%qG@{Tx0dHq|K`?Mi+p*|LwOq*USJKAti+A!d$X)P5)&3l^K6Z>u7l2Z_&@Ioi z!(G!k779Q)x>vX@F-WHkp(xd%2-uXJepmJ-R>2`>*l=|8^{Y@@1h2mt1%YA>%GP_w z$4FpT$PLRGxt5~QC1Vaa;RW6kTl!9M_KTq6pfZJEM2`R!BJ-#TP}X}vPT#okC}30F zS?R4jaN_c#a5%V8&>PE9@4Nc@FT>_3GwTKc1$Lz>P&P6ds%J(85IYy)803aV$|-f6 z4nG?VK=rFFAUIf9KtOG%quoUW=7n23PHil~(ArvMarM?%94Gj9L!n;-kg#J3;pw|y z#mw}l3|R#eDVN}7<+Mg5mErgr!gn)RM`te!n*g9AKA5~GVb@DB6w`_E1mX}9Gk+2C z;>8R7vo!mS8`3^Is4h$kfAt>RMu=;SB32$SP!3#}pA{ArWgzo_a&p12(g9~>9DRNf zMql_kM|O|?$HF=%)td#QxUTPQz_pDA_d-1f8)Y0zQ58u> zB356esV}FP+l_XpE_Xd~_ycr8{Z!Vt*YDoF1apZj%uJmbd&QA;_P1?uTxRe0_pv&9 zdcCo^@Rf9a`<6KPjPC4uqE21Cs6XQC` zDhUXJ#Fl{aakg&~|5*-OaZ29Z3f9_aB&oC`x%3z}SxWTd#f5pKhI=JUvX6@MH){ zNwEoQhdaT1JQb7#iw<~LV2|aOmNS0KpY2PEczSxo9H&;@az8#w~T#xr0+o#!^*xlmm908MPE) zzhcW@gF&SXeF1NKq#9A@wFzMUxoI>A+LJg{e2o;VSVv4m)55m*QaqK5qde@sY2kgUii zyheNZG#6FORRWOZ(BMMiICzr-WNk}}iS{AZL`4Wc#Argwa{0S=42yHWnL!GSR&l$O zRaVl%994~zPl_*bO^bg0s_cB7o6^CgrNrpI?DgqGbvc7 zy3tV-c=WUV`u*F}Q0Qz369<_OWI+w>u8v3*RrQ@Kcxx}KUdR^SJzB9!oXjVIpoIYk z2bdee+QGc*dmfw_c%0(F7l~gaPGzX+J4M(N3kon+(mj%iO`^zVj;4n1&F{%`Nq)u2AGaj;2(?EI%(mYUaVh2$MPbO;Jh z8Kz3(PZ&VLkd?hz>bbJGG^dzJSz{8OED5Y+ZZDu*#XV&o>dWn=?s!$QaUc*#GXm38 zj@balbH-cdef@?#fA#DpGc2>qWG-Q`vAeKfEGH&Qj7iMQGpm5?xG==_XKUF=7=4Yp@(bo*^GyJ`VgB)d6aH8jE^MO=k9IuP* z{AV)+3I@{H!ok!6*Moq7|JPG|QAU|YgHQhX1$SlNXDj|bFwhE**CxcKrk0ivI$J`_ zV%5N*pm16qJC+DTSnU=hTPRX+AgEy_R8`~S#XzC3L+GB%Lz4=W)L3J$sBpZmLXA2Lk-j5kZSk=dCFwHm$j`+#k&Sk?#_s_wK(}lNzxcw3a zdC-}c52YGG&hLNJI-C9UX&DaxC(I*Ir#F7S3Wp;`0N~lvJiT-Y4QSKL;9U|ibwUUy zgk^?#@fD213c~K$w8864xG-(~&f=b)dK>QNxm@8| zp_~_hOOqcygtMK97jX17Hg}9-*fHmYx3YFuh=(kqmK-UJpRA$?{(S1#B(`vK#h;C! zOYH6Kkq^+(v%iY|>(uH^)oyXm&R(lv?GHsA z$INpmVd{Wo(O3RemBM*b(Tc(itqz5%!>^~mf2FEwso3o0TahasKYHYkk_IdKh_wg{ zorLYXTVkGtVH+>hL#0aIuDf^dOYA*zCk+K=nw1EqG26$tu+x*72_Qvk_?-x%o(bic z*H|mW_gH8GNb~ps*pF^l#N4>C6n}w-_G~@-3Y2bn`*k%8+8mG3 zpJSMt5yo&Y;ZOl%U%u1?hl<>Y)dSxY$&~c)VUTM22RE3Pd5DdBV1SrFAUCow4*|C^kMMb^=MeOwCd z7SXxW(n5*OJzSf#f(*kEN_%dm&Uc5Q52?cy_~+@X+)WP|Ou0y!andYA_V5ssRb$WnN(Z{SklOC%RE zj3K!CY~m}D5~GCc8va+)`Y1K^5_m$O(5*@n^#>rHB;jn4GLXG&>fcm`_l|G7LR1dO zx}$`|#`AT8h9NU7|2%CNIQ2Mc^V{CD#|~+qpc}9S#G{Sb586!z_cT#Wq2ftcS25gG z_2c7~b+VgSEfCPhjQUuoD~=OGgHbbhWo7UA{;>ys(H3_x#wrF1Y6Oj**KW$R0LD)< zZ`>22lAV_q1d9NFD1{D=g43t~1J<~-;46uj^tWy6EE*}PM71Q^9x4Sn&-fAKCxjiIEde03 z|0yI@01e`dIePS6M?7H7QYs0#6}+ogiCx>FsX+D z41x5$DEUyTVkre{CGRus$A`Lrpcx>pFn@W7bO%(MzU2x_;v#9#latqjuB1h^87ZFP zgY7vKR1%0q8T4#@SPQ%r)Ls~TLXc^d$458?d&)65Ntzi9bXva*Vcjxi3-9gA=cj(S z!cl`X^$qYEOKWpvP55#KPXERE5kFL{&jl*56=4{B4O3{|{Fe_YcCb|@gc0Hc<}_dvK=0@=L0zN>}}eJ$ojLL=d(sTm+U(6|{F zqiD9H97$mE6i!fzdbs!YMRExYz@b9M-Z->52u#0$4kQKaG4H?=ZgqtvpiVIIggxH3 z@4mZk7N(26O5KVZH(i(`coOhEDKSyK_n~K~B&0NJ04r~w+b)i^3@rglw~mTcqmYk% zP#n3Bo&?_2X7nW*KJ1V#Zh=vRjQ#TL?*ao_PLkXcgJFdl>6O-Xw z!1bxgZMhG;ZVMIxC2NQC5X+>hh8GB+d#DBa=Us#m$UWuo&X9A3&VYtO>bgRFjN&6r z2O-9g{J8iU&=m*(*QwOG*#vl)yFb6nyfajiv?1i!<7vRP!c>|bt#}el1)*6m$XG@$ z3((o+br9DuoGk#$`C_uwvAN{DA|)Up!2)+9fSrH%)&td?8n7TpNl6ih%CTb?i7!Y& z0fDBjpbY&LK?ldpuVNi^HuZ0JnH^!ojF#5nLXXX`f^=(dd6txjY2)AXXku<*KobjZ zg+zz}e@=4QLru$*|EENp6qCzVw+?qga6=-&q)Lj>WM;$JaI&86)(Yx$&*MR%lTi{t zaWa@ky72!+?`9!90f1&ol3-dQh~?a2p#Zy>^z`(`i-x*DIC`=iE$t4K<2Sw(q)%mV;%fE|P<%fU(mSvZ{{bmxD09xka+O`n*0iISu$ANlIIt}#j z(%G_tW@@S4;lslcrL$7LxpAlNnt#>rG&IC?X2e4@v?DKmGs&L(LRZC!K7SNNlJFPE z;*8b|BzFQ{cck%hana*r&0@6Z>;= zPctv(0s|WW`2rOR6#E6RfeQ{?|E&qH>{n~GxtEl21>>2{Kad{I^hO)Ij&|;(PTRgG zdx~`OxchDWn_hR4KhW1x!_uzs4LRDVtVOQ)Jka0(OKBP8I#V+q!2l5qoqap|A?8A*b@+#H6G~YDFz{GJ)%FS_|>99IxDX z?)ePwSifO|ppX!8T@vx<-r4W#ieFsM2Yo`WNs)8hv|&RRxU%unQIYzINTsmjlv-IE zvtg!CRaKSv&6*yt)NMZN?fnh6HJn+TpTZdASS7P6?l{sB_u;_olmuqfp=Mv~mhbuT zqZ>bk$(aUEd)Q*?pWI6#MnPC3iGKmpT{(s^=1KR2)pJHL%zOq83s}V0fO&?oq#GMX zj9X-p|8L50qOG&~c&HKgGjL!}VMie0^&rNa1YFn8pJ`+Cc^HMx^O*9Ez!nUV*|wUS z+bTcADbi3Vs?!UzBNM+gOdv`nOooUjZ_U&$9RU{HI8b8pYiG^Xi4NcYY8z5o;nh zOatf#Ucuh?>C-0)4H+3*z9sLLh&6(g@;_#S5Z4%hLR{;_JU}9V$D0FXd+oL(0kDM# z<3cJ_$Xn9(Z05c*HYMwhSh(`6aB*?LjT)Q)f=GdtC~?zmkr)AO5)MoP37&R1YY8a6 z3d~MS_ri)-1}X)O|1)hd{Lt_k78ZXHO2=BB*--wV=>$NV!N*PNNmv+T5ECyw272Kf z!1hA5H5#el5rp#hC#MxJ^AwUy^g}0`-=47b;i!`T4L!18=t;wBk0W6$^KRc>u6%b# zW)o_zjZkvUrvdw1v~}IYp9v){5`&enls!qC z8XZ-`?EfglqCDFvn2Vh-2%!$O@1cI&72L{2v^iqd!AN7mU|l`%H$4ArXe1iWLpwb7>Hbb2SDI=!G`>K#4`r9CkL#pX7^3_F-N3?lY4&6{5-qxMBMvP( z!ENPkJhmV2g`t9xaU%Grpx7^eXb>O%SISnSe{%3dI6h>KI@jB1S!!ry%d3$%@!{X}RF`?FUF50vDR57jn;T z;=$Z@+jq#JE65Z!&O1ai!Z!nlgF(f$JhFgdv3VKB;-4PO1$Zc_&X2Dq1nn2pQ%T>!?8Xgt02(&p*~CkZzFck4*B&Sg3<&9n@cl45VUqnST!Fu=F0&+_ zcXoHlg&s$q1#FjD;>uX9aJv2iG&^}tEfrJzU>~f|esHS#zCa!jzW#XC^3W%9G(bphw@GQzy}UohqWa60B{Ul z|1hHA(6skLq z$~jJOK>;w5>c=S%XUsxF3Zm;{JQFUIQ=)fF3l|r4`J>_4&>Bb!S%g zE8rrYFv4Sk0%oCGdS);^e-!s1Rv;GfLg15Um54|>NjgV!ac;l9?CH&_`Y49iSf?3WJ zGi-NvrJN%)zP~CF`b@Ie>hGDf$=y#VyckSOc{Yc=ng%@&S$*uwsgGSa0z`R3$_xVN z&G9To>Kz8@o(1|Io56zlCIShNVxMBxuXq(7ni1G;0xC!(H?d`mTLWg0M0>cCmZYiq&G`wx-m$^{zK$zmz$sm_T{XKgKXyHq2 z>FNU_gu`f`$;Q78b0`{14S{{_(J?VGU(eq`z~l+A&=GO28#brWiHTh>loAI≈Wb zTi_-bVc@r~U%yiRpd*DgZd~5m+grFTltsDmr|WTh`^U2Dxj8siP}p1JwtL<3m{L$W zb{3mjiXbHin4EZO6eL{|A{sJ9qnpi>;ggNBl9FAqx!`G@oK`*%(09o;7(Qo#Ie;Yy zNnJ7XXpJzH1x@cz+X5;uxve2~%Y^2?TSXKmq?f^9jH67v*~B$TMt8SuW5fkRuYhnT zh8LW$qnxr%oM07UCO{5A6h>gm1`s90+%Mq?3~e` z9niHy-{Mp`S75GRy!glturXvKI3dy^k)ar>4AA81iI04ch@jXH{}aq`Z37wBL-Ohx z7VY_UK-U0AMMDvQ zVFENTeWyO{X($5!(X5;^pO17nd(B2>1OV${X&8sZ0}F~r^bMBSVDaQ3kU;}h6p>Htm;P$wvsy`SxA?QI+VT2q}batJAX)o zWOKe(0PUm5Mi6S&Qp78S?a&wzv&4lJWHWH%8}(oFW46vP|DBb;BV|KR!U08*^jyPV zo)g3cY}A(47k&{>*<*+{&Z+q+4epl|18G65oTSBp)Qv6{UYvlRq zQ7a%1m;zbhOC_W*ZU>NXZk2|?=LQiG= zci0u!jR+vXtYj~+CS<4wSEv9E!DNHEClgZkMTcZR8p1x##H$bhBu!u%Pc^x~tga@Z zE05fHmB{#aMDTrVFbsj5+}zzzKB|zk8iv*HJA!v5%Jkox@E|SX>K{O{ar2`0r;ek0 zLX3B2Rr3+q9ykW{7QbgcDnjqxyYj5?kq}Op_%AR97Rs;6ho%XO3I1^zNW6;p4Y;0V z7w!W-agf|@gqwh>F2f6!j)XPrx|mv6LMfpzOyhQBD4Yr1Zx7|3Jj*`ka>mKxc| zR;ef?qf!)QEm;dK3|S(3L?x{#iHwq6Dnx~Z(jr(>NOKp)rYN-9t-$^l`@CZ0OE z%*hIF>>i-WR!IvF#r^eZC zz4Tv64hsv5xS;)OW|%fLP?(O%>y)XBsXl}FvMVK=f!``(JUu*aqT&!o9%t>;Q?q(f zOM|L3R?IgwMStg6W8T)^_S-7WzaAYk>JX0N$Y|`$6H*X;YG&CY+qG)21kORA-bzVn zOFb#Qf0=z^J}PvnwGi+C__Js>bb3Ged5S$rIf$~!oH~dpNYVm1NPbiN{5Y=m7T0N6 z(ZQJ@+=5O+N-D&S!|}p;TXQ<|gSILQt7gpDu}ajHBGaHpANTW119qOf21(o-i96k8?OAB=hsXSKWUGXh|!2(J)A&Vb`r+4+SDcsn||X ziUxlEI`3cW2NEQG05_4?qjTWklnc*~BhUHUhZ)yvYpMj-;p7PHMLbN#1XiD5bL_LO ziZ*TAZWYA%`Ey);7yhTWiEo{XmQMZ2J?kjcQI_Aoyd`6LD;jyQ(57|RCtE8Fc_mqo z%?YQ~Lh)Y0%|9+9q9;#f5M8NQosrhG04(6a^X6LcHa2KkM$PK4&jW#{YRtJnA10Ls zi%c*)^9TcP0M<52N~cE7`~2?R*29N8OIRA)gtFjvx0Tnobh6kK_H5*9$lYH0Cn!)B z<8u_%jbafMlSNpRaVFO#p+LEYLlq2tW*Uz1-2ZG-}4Ye32cSO%N z!gb2s=HLLXq_%JYzX*DQO`A5Y0%GvZMF3&I!2)RLVRiqYb--j7JvlzT1vsUF9FM(R z@rF`XWS=yEXvgOn@?;P)fTt@d`|0d=NA+uo4`9!og;8hXL4qx5+=D8AIxHInL3~& z>-q(bK!l@e03Ar?{v}+fic0N^`lhUp`ji-fR)^zej)`kozM<%&=yDCSwUgeGk8L;i zLXT6FVv`LLJKrY2J8bz=RX*%Zs-aiWcDDg0sfC;Uhh426`aekz;+jD<{XAe_tBVhd zL(=`ujq=VZwYeJ|GK&*pO(S`~ltz1KNT4UcH)#4vtoamud84evAwybmy(rueRU{-O zx%)*;*~48OwLuwgJVjgmLm7&#Y}Xykq8=A&Q01 z>#whR3w?qRnhHhVnaA#C&+C5mS%dc7y6r)_jyzJlXnf#8^S_Xc{c^hgCvLrKQtE24 zz7U&()Tb^rMc8nDMVO>SUM;P$J3JgR(J!xd-&rDtxWmJ(rp7)WaOnO%KfimApD)Ew zhMv}|U8nt*yy~odyv;r^Jco5%bjYDBj%}Ibpm0N3J-c7MrUbiS1#W^M01z80kHq*H z{i51b>$iLl0Z4p?~nfN z-?8KyxK%ako!D7?Ff7Tu$cj(q{7mVDkwCIwsV0>=bvplhFYiyCetn_czcE4>w&xv; zjg8HL$PtY$!mz!5e!FSNm#}6D#Y&wWb&=F42^+>k48{#{-%Cfw&o63MSHF*e9i*Wm z&F5s7a28@1MM1&qYb29R>B4QPcm3zr`#MMdFx$_E7Y4lAfDTmBZ5{?^CO_@H*EIF- zRlN`0pAAx6Cx0U;jhz3TnyumlxFo;g*(Ri^fhF6|uHBiwpCR`Ge5{(hKZ`FA3O-6gzA<=a&@@ zkw6{>Gw_ioKDT0N9YyD}e0HpOuc=hrk*6bx6e0W2>tE`NmF<~*%u;$v@GIhNP%X>G zR#$)KWUb#A*+&8$C`jurm`=8rdo4hRM(fIbz=$V=f?zp}jemI6%z z7>YXzB=uyg1DAFr%x;d6o3zQ+G~{Vc19+5+yL9a;aJ_Usn;EJhGboUo5@`4PjjE#2 z{*9GTfpb>ja6tIp+x-(MeFL+Qnl&CZgldenIBxC_>(>qqZ&e*hbE7Ccreuvj< z*G^oso%yoJrJy&wT#A=M$Lm_Q=@KZd<8<*4pN-gA-+ub!K|5L4zPwFgwT_Y##9zZC z2Qm7tS=fwFQBm{C3=?;hb7a2Md^=IQh_>;k_3sP$o`6K}jz)2Rh$Rmn=;!MOyCV_Z zZ(p;R90c@`jS=chcs6s)KP=Im&~dz;`T%o$DHVZr+70B~m(eN8^~3sSw}@8SUBz`> zwtBSmzodTgwr?W+7bJxXugAJ$)X0O(qFb;5p?=BoKKI zRjyn@nD2zLdRFX&^SbNNO-SSdZTTP113tyU3oqt+Lzukd=uxwv3eacKJPIvNem>`p zkX^Ji)X7pUmQSjZ$R)7nt3+`lnA)Pf`h#qS)gK<*0*KV=(T`h`##Fp@;tdvThK@?G z86M!d=xYAbK0J1S{*jf}qO>^duUBMk^wxTdpix@C7PQDTG#AIgk|kk>>%FPnb!SOA@;&S77y9u~ zzDx=cN;th7+DzsVjIAn&$QWGEi;%KAipxa#`KOogHZ%qilt7S(^P8LX8o7m-$Hm;K zLLS_+!A;og)6w=M8@EurBDp5M01BINCNMDS5-dnokpw(5tAzV|VAX}-;H5lvOX{&f z2>uF3F;^jHIE-IMC8ehYA5(g4&N)CqQ-${8FfrD(SL_*pQ$UqcXtb>_l*S;WsMR1= z1^9ccz~sAHSe#bztE!oyW}y6nMIsYy=vV>=9r~rXT?I0sZ0IaE$)ci5Zbi@P#j9^0Qn-|qgPUYoYnm9Z*y~h!#qT)+ajXvd4Yzn(y4nT1s>pt*$m_*w?fQPGcRe| z25TDf5(20b$Yfy=OP3Fby9x`Q>ro&Qq2w71ApToX>%4ut1sA}Xz`)7h_GHz^ZY%L5 zBBAp7%uW?-OfocZhC}Pga2wgXQFg^eYe_N~nsEEJ*y+Dk!}H#TIsATNOG!?W?WTpr z<+!P&#p8z9fwdV+(U;_YCjz~^i(W=F=Fxyak0Km~Hl;onC<5xp?@rfFNIn1RRn@+J z&e$w<6p{vS^O@lWtrY{av%x3oVzwm4(UtQ_JIl4ouFiR6SRR1{_BDA(h@%W<&Dxvc zz*^QEBAMc}Qrl1Nx19@U6Y7A&!xqlq{l3#Ff!HSc|AQx`gc)E({Rkf(#wv4izf^av2Y>^#+P1hy=JKgma^r7BCp% zN;tW9f&!uI1V=Ti_N(USD8MNj1%@yb=O*Do3emboR#A zC|Vpq5T^TSHYH8uL-uw~BPAszNsEywQk>61Y48vB!QikPa&K`~k|>#oC_==>w9VIR zDnCmU9w^({^20PHlJq4H?Ze{}q2KP`)SCS$^BbG^EnT~@iDZFXbZJ5{*H8hxaItXO zwE|Q)h9<7aOCUnnzN+f#&EX*R(1_C3e!;TL%%#becAS|10B_46C;Y3iX^&?gB|jp} zwR(dv>*90dxtG8sEkQ-`BF0G_bN;I$Er(w&K>XMo3OIZ3Kz&I#?#QO8`i(OzI|oIq zPkm-XH#YG6tmUYpa=qvikAZB*kTZ`j`7?9Xl&V0!evQHMogo z-S8k5E? zwUKWw+wcA6Nky(>v9B^*2qje*Qasmd6HhJTA|Br0rt8lV1J@t?xGoC38>6v4Bz}wp z(d}3r%SG;D5o@3wa&A`t7t~E%1anSAq;^mKg`_z`LU$D~o7Rmc+%Rt~T@ zk9W9l`({Tt=Bmh>JUqVVtRHdcZQT=#l`PlZh*!LTR6JU?Amdx4U*BdswiaS6y%rPm z_3=thSe=ldu&{mHUM1&W%2>6M;s#H>umR5>m%E2@vq49a6x(o40+0{}m6SM(cmmlk-7h?Mf{1gaZe z%RHmhg9FZI(~hs)Z@GQv&WoD{V@bwG&S#d+*vcQUToJgy#F1l3LCnmS#QPCQJ$}?S zYS{3LgjS=#dr_b>)ckrA1rqA%5z*sFCU!1(v4vC()B|P#&-)UQ`7$Za7kDB7ROSy@!6UlXRU+GGg!S`o7*3!rzTJ^|yJsrcTJD1gs zno3)*vZLfZW$$-aw%C=Iyg)*Z*L^thmju@f$H2MRh-#W+>oT~|Po9a4{}2`So0*}{ zcs=fd$AI&S6ba&%XSo@T=&twrak0*JDS+g&K~f1;LaGrb@t0E)T6HOn_C zRZqPM(X+&ohr5x`PgKDm=?SekS_#g`tyvB%KO5r~G7%55wNJ zHSSSWlvnFJA`{%NJyMsrVh7CL^KLjhI@-)0CwW4kS`9L`t;mRMis;Inu%g=Pk>sE& zR=f=VjjHLPZc*HUDuiEi&5R}Jt)kH1(ZQjf!s}bf@(rGjqoEzWs-Be|Mpgh~_m(Mk z&%SpaU%j#9q+eMwV3w$u-r?%Kg+={)F6X8uC7Xu*>ruhTRe{`dpBGY^xKx(j;s)`O zq!knrV0U+@-XzkB@=EtpR`rO@6*az~MIQ28DX?wcPX6)jiIIwDh8gB@xjt8B>x%q~ zNhQ#)qXdv6DdDV1{93zYdVIWi2`COk+x0resW~X-FWte6*z}2}I}>I%?t=(G@-Q!) zKTZuYu-=p(%i_?3_%}XtuuT()j+d#aXQPW;+vP-cu)11F>Hme4tC3X>gY6ymMq#8%bz-ErQrb%EGch%@A<@+O*Re9b)2zf-m}lA zHbvaX?%?d_ZOMQ{?;O)HTT4zB8%*oBVnO)QT}D|q54_z7Xm0a)XKR&K-Z@0OpATQa z@J1>j2@|^a_4Vroq%nvG6m=!c&7eK{s$ddcW$K5lnk*D z$DJh_D0@zq*r>Xk`*|kBS~UvXB&3Jg`VJZ3k;HzQUPQ=>3}HC8&rUk1*K1?ps|gH| z^3-PKJkG!L?}_PiXryS|mnXTHARgwK7&E0Ku`>lXssUb0 z55*UT>IDu`K|LoT7s#+6|MZ!E{@GaJ^{l&)5##sF?C@C~R6)Y2B!CPIL$WLenS?hK zqbjJ-u;ouqI)@}ZdejOK0d2x0ywFhMgm!r%aB=N@=u>N;s%=N9c^y};mg-AR8R`vv ztOD-MC*Ab_Jr|jVNSNg}jlWJTICphweppr6XLlf{9Qu$^>%N@K*sx>bE%&VkLTAAp zFF_0gU$*m}POeAWRUg|Qpq|`payf-p8OB*=aH)ZUf(2u-59+NbC+tUCnY9MfXgH}4 zRiET*iR_*B?xD_1Zna@+KFva>{B7ljshC~%CEgg3+a#IXqaQVav7|mimYp$4GG4NF z$aw4Wcmsd63sBh<>OU*y4WL%P>3`X|!rF=;P3hX*Dg*!iI`3nqDD@ykm@J`8jSah4Mfkvgqh&NmD~;CA1~rvUA8C zh;o{xCZy3p_HZYVTc&bh^WNo+8AMW_^sIg@-7_SdL%_ezqgM5l9txO9j&ZOjpL8D= z+FY?%J|6R}1xt1-*srjw^wL6Lqgz*5fBG79vz&K(It~UB8Yf0*?p4Xqqg|HvD~zK6 z238L3qa*Dmg8Ir?1esm;bF$*9)G}{OwC2&i2)Am^P{30oV=g2A6j#0Qi1a?g7((Rg z-5;=@+UkPWhB1>-rBHW^!$CwjQi(&YisPGhtmnv)MpjixJvO!wfq2&DS{LZMx3JOT z1C=~J%3D8yvc2J%CX@5@eR$mB+tm3f31nYYvGkSq3+C4}J;f;pe@X)hpx?|$?(v)e^2WD;h zHbaQy7vwO(ZEppmlg44?^NtNsBi%5GWe@|ZTtc`$R3&-zc#`1j=zm-{Q&<1QfT7Qy zuetVZ+mKhaoWF+;9qJ_YnLrRCY#=dp)0=>weylxzr{b)IIXW|{mUxD7FGT|+b!^k- zCcyF$`1x<;<;!3q|6(+YeR)lnjbH3EWOQoZ>4$HEM<2$-*>2CB6|B+%O4Xku;T8(e zjx6N`38=!y^};~tNzrFL#(COD?IoTX%eOwpMJvd0+Hl=jL@e?N=cO0_dE@Xgw<`E-DE~9^5!i%v)J8AO@ zW-m<88sG$(wnaUY+;OP+1Iol6^S`mWP1tyZ^>-OTy@wZVW{U}~Y!}SuND3oKvgylEgbsSIM&g zl_4)ecJz1VeQM2gI5mT`1<8>5Z?V%P|2}zLu~LR8H>cL`!A3B99Xws<5astKKAlKc z{ggd*bf_U1IM;+CUpSu5mZFh*8e}{aChO>lp%=&E_p<~t!0L`5mI8wM@A{*nP6jX~ zDFG~9VZ3O|H^1LCWlEASPBEU)mzHlA%1pW&3>f=~nod-3u(Vb)asbAr?nL$zY{f*SN?|oG9`z0J8CtS+AWNj=S3n4k4L;h{Tj&)-E z`hMj`2)sjGIJ&xSgZlM#Qm;Ez?k0UfG*FAsK(UXN`LKuSG9X)v-E@R7XL1FU#`CmD~-_2klunaiH8q=hop#hEiOFY61%;SLBDh} zFlrY}VwLN-sfk6mLA2yV<4-lnm%*7Sp@=-YI{<#KtlB3D4iMo?18yk`>1D+2+mVly5>5j z7S!(vv^_hZ@j&BZX4p5W9>!+uJZ9anaVlx-$2N61%PP+MIr%VCBMB5w#0kcHIX%wJO8ULhb4Lq$3#89F#wg(IJScvwA&isMH_x^A+H zh{m^-9qdB?AraWZW^UhW^lnrusgBtZ!Y%*>zu`pPLM(@zs6;s}FZ+%G>l5RTX*!Xt zO1;5r#U^S~cGVwiy}Vc@CipEO%M6*xKd=qLhN6J!HeteqQ4a66r~721M3X5fG$I(< z{(&jo_xIluA7+H;zx7ydh;33_PeNIPU7@7nDGSmLIF1Hq&BD`|oZRhe zMTKCoUx=6z%o8p_GDB*rvTHGp*_6GIM~huBJJ_EcgOUhVTs|prV5bd59wZ|Y1Phbo z34BCi)hx>ahwD~10Nwa=5<1gZupN=oP-V>ffkqcSJ2!(jZaO+T++^dw zq9}rq{}n~?SnD-d*lE1Z6%tb#@Ige@!Iu1yO9>|B=C*CyS_LtArV#VaEfMWep)O*s zrgVM>hgx|5+mZ44J<`1kp!Vj ztV9Fv97z}E2H}eK?QyiKtowTAor6C^7}n|Xe&v|)Q56b60x8Uqao?tW0yD-Y{|}|> zk8k;))mRyax{Z#=$l6%o&REskLZzX7m3U@QAMvf_U=n#;WXiZje)UrZoRhx8I1+${`T36nBKdM!I=GDKd?>4Ld zif4AkS?9D4Ojyd5rf}oPRg90Q^g?nsY{@@O1h~i2!Is_%Yl(tLKutO!QfV5=tZ=#j zQOD5|G-iFdQO1a3iDxB>2uYRk@|qFh;{vKKUm@%Os4v6Aw#eXkexRV_K#`3RzAj5^ zu8@{aenBKCIXgO;yMe_}V6Dp@z3~w9SE1Qi6}L%B0Lqw8k0V?67SSZ)tU+{SQs@90 zI|0cpfMs)VbQG^Sryl#F`NDp~VcJUUu%Xw#HiJBs^yE6Zq_i5?)Iz#HM}DTavxRsP zWG)y)j!0zb_t?dz+RF@}qTLY``BYt6?dT|G@gJd-4#icCx*(XMSo5f@1j1v#r(O-A znz^^Wh{MR{*jHB`VC#_D31WR=eV3^9AKWDw8_Re)YHt1pyldAV!U^&d33qGH9c#3H zI=gIh4;S!ek+ShC@my018m*!cpufb%-$*ht`91H4ZiH~TOgC`7B(WQev~ff}1fv*KW669LqnrH#U6(h;~g?jVNupda~F|Y=z1Txbx2m(q+|0vNJ zt`jn;G=($SaBE_w=HM%&m5^Dbru>bnPn_xem-B5s6_TL06aG>xQ5x%PF0SSI;9nC2nlP+ z4p+GSU%I{%H9W3wl>XA+Cpd7zR+ z)!uQ?d^h^d>!juUybz7&Z{#?9ts2mL?rg`JsnOg~wD5tO4u8JLAkA6RDmRoaepUYN z{rj^T1@xmibe!|!HfXt475-imc^2Y(?I7^_$jlOL)Q#NOrw|;TUAz_Z)kr$ zueH}X%l_SNJceN6XnlR<5d{0PLZ3D%VWI-3a}Z)c_NwF+2(H92Si*MrR=wXYnbVU~ zT13ArW8{u9LX}JWLTQ|{*20P-?c&N#&hDA(yBz$0XLS<*4{MqHiqt}5QXG6|?5c=; zV0i6k)`vS)X|{HDEfugXF~v7Dm20D(n>Vm5xKV2qY93PlBD%A-+V)S!IVS+jcXI}g z`SxH0vNic|!g;W&Y^+obAmO-<4-}7Tt*T;M$1hnJ4DNzsOD zcq_1i!@Gv!L_+x?{skf*^ne48s{GlrO-}+WgnQ!U4OCZGYS-?WRq^ogWSPj|Xu%ya z9CQ0W$ZYxa-iM{T^iv9|4fOR-qCJuvjx@jH&slH)JS3!ty!6VZzj8;p%lU6O$%k(z zW1R!-ijiT~B8Zik3G|DBCGFlJH^4~ts^B9YtQQg}!~Pa+47Hbd5eXpRI=CguISpwMLZj03m4yPXJQc) zZbI-$^(hWKe${cSY7y`fD!GL=2)fg!%ZO!fLat;bz^lmzmfwyh#US4tyVpWb8Z_y9B~=AJmVsa1YML&?oTs`-~Q zJLH85?$5tLgU7A^+PlF`R$XRrVHe9mW5@arHoty&$IhL@IsKQdSn+^-u3N~1A5wFt zMQ5X2l?IqTvxQ=hr{@ZzhptmFx#-U_NK+{2g3%{5BE@wowijq)=%Z-0^OqCDs3H-B7LTv+j?-yxxf<( zz6Ko*SY~wkKwGPUfdQxM4Y}&J<*S`UVY<h)BVKHgljA879A*-)Y$cO!mcckK4^3J~&xXXF?p}M;I>)5reX$ZE_NFWpI8Wl=aa~lUR+}YnV z(kX}S_wSE-W#!nreKV`u-XI6~^Oye}rPi5W{2Y!K8R6&Y*+RrTq*sB?3Sf5`eF8Oz zrm00rat3GY7PuvuQ-V4>`dV_v;NjXj7#!)xpK|Skn_+I8f9erP@FH8==7R?h#yhD@0^8pkd7NYKTWcdLM69IZA9U&6nTa_?hjQ?v30`;*ANL9FW-qO+D@DAd4H_ImsmFil@A)w#LOo znW$mJ9sRG199uE;#LUfThdDAEgD!_U8!lDI+-AHHulHyqcDCIGI<$%!T|{ay3dHV- z{6JiAOdg4STKV+(b9XQ)A}tjPX;*pB?wlZ6rYiXPBr{LS)0;M}InQo9T|{nQ z+mE@qT2reYUP$X>FlY1TX4*;@elB(7l}aJWQ(&@n()IP<{8aT5J6i_dp%MhpOn|H} zNjF;w5G@n@C`c11mH5C}{aplGpj!%0&JZmq&AMsY(n6y{-H#kT3^Fs&v+dYT-)Ak> z+{I$B{7N?&M&LXEx*Xi=DoY25eXpqK($~@h8FcS5|7K*d7%A&u=H0uy22`I*&;MHu zSCu&ShYXo{X$&Iq+V4r@QrmanAV`nl&w+oM%P1B}i{;&m{r<4_>Sbqw&Dfx#17w-N z5ehU758}A?;FWpoHo(0u$i(XgA6;Gju@(P&>B>(C_)CzGECDi>qY2%jtojH3qm8|0 zkKlI**~7LK1D{-%lC+`d`AdY8z|w`)fBHr@W=)=;O6<&&p#MnTMLI*X*!8vZTQJgq0mlk5*#T$6)j&4lQ zAAJoNX53SC$8%Q8QYVQ;eI2*46KAtj?_A7&V=h6BoZ%cLHpyat#T}NI))i5`qqda| zR9aRB`mwd9y2_HFZZhGjYpQb^Q6w=j$ku7s=Il8ZoTw@WF{^C5F+`6^Y*MT-<9Rp834w@%xniVAeFP58)u27G1b zxlaP7b&^mF3QSqZOuTgG2?@x>*W1l{g(-qfBJ)xyhUwHe>0d1AX!7B`ht>s2vIB^g z5xP7SVN#(L0$R6&!b+mYLlxFYTzk+w~^2FRswEbY7AKz6KVcu!XsiL zDbcu%G3HCLMGL|)+&|SH>lyt&H;$O<`7VtK#KWDkH%%MTg)_rSvUuj_B)fZDhL;jX zp|Y~_$Pzo~YCuG)I*eODNHv5mTT2A9916PA_9s{8hOd$*c#Y}7t!CAyOKR_`S1EQQLUk8v%Js`tLh zlLvL&4_id3l>)tT_tvc~6e<0Jpvk8&8myt4M1Y%yFKtR(#3A^ACJMXE7{g1M;Ij3{ z^B2xuxG+UTu!Bv)w{zDSDBu+*G)ub1-CBP<%jq1*-BC6zrvKga)6*1pU%fI5s34Tw zMM5|BXLPWEG|)InLwkeP_`!H$Q(W1y1Rgw|YD5P`)_Urayn)uX*$4 z#6ull5J$P*PO+h?@QRfctxW1Ey~ZhNX#vbr`QsF${Z_n;$DL)>bSl>J1#$737XkE5*qcnd<+q~N{;2gEaZz9uja-SUp1`U@mGO#6iCt>r0XW7#)R;?u( zg5jZZ6}@aY_QI6?eHYVwe23>NL3dkUQC^g}hYu@Gs_@(%7fBM{TWz+%UB+GH9v^SF zc<~l;zfkgsqE5i>_Zf##Z<4azyUL`0dsQ6dfGUy`jKHOt;spbG2+;x9wlLONq%ZBu_nFr|e(+%4i*v!76X(C`w{No00;1W3 zq?QN?F#uvY6WW@?q|Efs8O&F@dGlmb)k}IUYzpK2i6I1#kBqq;XsBo41q{)Qb6or~ zAnlSC!}6DWLWyjHrd9Mb8!gsrt!U$U^CV*3#Bqtey}4}*+9Hus#Ds8;DBS3ycG%t1 zJ)H9Nz(utoDK0_YU8khwnru#);+0l7cW&^#@(Tlo43D5761);z$tCPp)j%e-caO3k&u%Yc;=B*rV_ET3RHu4X#O>MT;(;Hulq( zPfyjl>#S+Quj}FEbq5Rb(kVM{+|X6LVUqnvI&x1j4q!QLE3h!m=k+m+T<)^&7#pIXbZHSpT8kj45{S-wTSbW;WMc+O= zocgi3q4I2+fcgqfg?k}K8JV$-768*eYcC)xY1I0!O9^5kndrV(3gpu}Hr~8z%2XYq zG3MX~F^PK}4MO!6qryZS{_p$ldj1c<8;k5;RK1^iKUHxRPgUuAjWr;%53*}(uB6=m zYi-f%eX;G6I$NgY;X0t%6Bh;H)%E!gw@@9w{rK_uz(;ChK6+>!q1K@F7D5E}4?7YD z)$Saq)D`_>Z(rmCNY)B$-LYfGM0_DGDEdf>VLilyH6=cq-vFp8%paZE#poo8r>BT6 z#S|irF(O)5q56U9xZCl_G3LEWGB+Ac5)wrdE>dfW%EYk9oQE!43PBYb<-}OfUrE#s zr;yN_JhYz>*wVKG+gxxo6I2@l zIj5SIHu@)B!MZPR>hZm4i-%kR09`F@Z$V~kCz-B6$lMeD8vUH=h!O1=#30jT7h77+ zJ#@uGBYkGmhV@0qEp}R<=WUdJJuyq8O8Td<$|hk~0?WK3yJ}Uc)-O~Ht$)FfbMuL9 zXxt>8E&bFVI2trTvnu%zx#>^$-FxuhEHu`>s37b{Ubo0QaZ0-32vUN?CIh__}?@flI5~W5BsQ zn7I`Sk~I$}L<8t;mAZ_|VCmw;69l{i40>+uucj)&T~Fzjk~1>=P(ac{2tc04QP9jO z;B)uPYpc&^ltW|aVGKSS5;6nww0)a4<9YvTl;%5}yw4&fvWE6wx_tR${4^-dlc}f( zSf0+J{X1n4T8sH%!0Y}&y=`!YKZXPvN_BrQ@% zJqd|1aiQ#4kr~s0gU%}@9{lK7_p}Ce-NVOGMw^K?pMTA<5A@JZ1W3!#M7i_`D7AmV zPSKpY^=y8{?|ITu{b;eYr=?0qVxVA&+UpndUa|y{xG@Xg5!rj7a zhlI6K>IIH5^Zpsl_WyP&-3&WDaDmilXYh60vF_6LPPZeet80|G<(v0`t=8>WA)U1Y zL5nDZPiGh?vBR`X@vIi&@}%hsn}m)1+}E5Ho#*bL>ZX8k#}oDr$HTA@ zBepT)EyBC=6bHxD?UdSI-rQlw1zKy2Jq-&qiOr8a$HP9{OrG72kSW?O!?m63nm8qk zVe*VS2W{Q2#xLO+HgJ?Ef&Slra7O3NUqSM_+MRIon4CEwa^!!mj(RkqY09nSP0Ry| z{5^E_yw)k3DAQ&=+&7G`m~NO-T$qu#W0dwqcN`=274&>*-J8na&(Eidr*R8(9&KS1 z9Jx)Bcz55Z7bl z47EvfHeE4oJ_3@<$^GVhf6eE`!8*F#2EYCGqV1BvDbsvHM}~~8*}432l6kiq!HM&` zC9P#X#r(`Kca<#djxl+I!?SnWU3Krz!gZKF$M%9#>Y zleieYC7u^v#@Cb$B|K{hVGjAN0%q8R7OylLZ9mgHV;J29HT?a8r?B%-(%s5&18s&!7x2Z$l|>V8ezD@-?Jh3v3_g!javG?to{kPIpYCMQ~)+e?vvz zLRe4zj4bRVLMoMm_uqf{v?Pc3cuSi26)Rd$M1sEMK#_T+yIP)TXBxzdCNrYvs8LsE zzp~p^T#QtXEuxDlHA7&v#XZxO3G}z%5h$6Di9e&C?Z1wVJQzJjw8+Wf>32N zY|V1fSqPuwTQrv<5NS6m*pQ4-9F$5-jE8JegOdbxm6(R2?D(O^Lz$Egu$@E|p1 zQ&CX#Pd;@WL{5`*nrp-1+D)TZ94YunDgU*!^eiBbPIo2S!Cb=D+~W*jH7YDl2Bb?) zOx(#_7py2LfH4Q6eUBZRPDy3`&g0cN%$q*AZ0B7JDq3A%X*OL(97|x}omX?lxKtjX zo14OY@V>Bc29^07)6kIkc-WOvkhg&cu9W52+mCTf1$ktUJ6BM%S2F4jG)$Ez$&)jn z=N9TyKwlD+x>BZ_hU!}UPpG!7!p%d{qUA5r;RP9w3$_X>4Ao@~_MqVcE|P2*B#v)E z;kY&@S}$}u-p(`vm!Tzj;uI3|y~;B3tIP-k7C!^IA8_R-6Uw7EytX-0R1h2$6dc^! z!Xn$W*5>-!;vObuW^*C9Xw;t4n_eJEFhKWY@1u{ClY1%F;O@Bg6;W6UhI!<=6zzRz zib1WyjW=wF$#0xw+QBNm-+t}|6HQIK9v`!Qt6+5CNFt8_-t$dgsh$*KzsV$|p76*JP(%~l z=9){4p*@jUPiPX+a$#0%RA$GpUuKUDOD=-f{{$7WEn}G-6-*tR?y5qn;Gn8FZA40J zNe3Dl4ujj=-{_0x8e9$WWBVo?@{Nmg@oH!EIkNK!kmvTt9%igp~; zWDcNyk_bKNZ;3&t4^LAU0P+{W-Kw_odHm8`HU}GeL0raCyLLS1I8p=#4Ek6?#ovkI za_Z-A+7p~mEf`LIkxhR>nJy+0(#_|+I6pW(UVl5IHUeH{e0BJJblJkEwsCBD1A-C` z;JI(zW6a(-lv5q@nJKZC7)GxiWBu^NTIJ%sU%Fz&bUNLsRf5e2oI5uM^;Nx-v@%m? z%m|>lO3vTJn(5tfPae?5X(+B{ty=Y6xNu<%cw~->rftcWD?jRkCWE`@f%iFWy6>RXP$y0`I4 z`f_4RqzZmdK;Rnfdg;19$7-u?(C)PnZu&5b#?bVBmQMOpQIOV7ALWOowQ{!(Nb?${G9v*(Kn|5oYw+G+cNoqxdb=plw zVUR`?yp!9GQ38^7YW?h-j3N?!jdNkf4oNNpRSzxyyS1ff)LN>;^xb`W_g+Nh)hw6Q zFOx~!4&AR=Nsb(^3Bg?hQK;hWTIPL#I=J=B|G3-Den2(JCI71&3~u(H|Kv`!0(8|_hH6R%lv;R zjDsJ}F`oFe4GxnL;$Oj=k8`KC{aAAFh~s#N_ipyu{+?Rh$~J$cMw=ma84mNH(Jiz^ zZN|?;P7qo&OQB#ABC-rQDp&;Pi&k~E?I@LWMhKYC8CGF_E4>f=U~CT}X(hQdRoVLi z=L|_H5;E%Ko-jXxClls;Hh;l!OvH;2GJJXwcII^_4W+t7UBY-NAy1t-(*X56LS?g4 zaoly%m2w;lpG{PcZRC`6LGIQ14GK*olKnr_OeB9arx0~NP5Li|=jeW$1I zL{q_z@t`#g+O!FgMLh+&!ocf~$^B{J;^HFwEfng(D8`COAf|XAgcb$f5YI2>iI)q` z0VLmI2m+S_nPCZl*$2uc*PNPPJTmLb3y$0iqqqmsn?D-xyCIF8;wB{&W7V3xv09=@ zeB!^E+=6O~0;7G04ifp$?AuIKv`ej4UN32ci`j=(8A`JUl)=6LrpJbLBXj6{nAyD(=}78T3wkjTc!GH8!AtPsNqLaGbO;Xz z`lx^^J6tXjYzODdIaW7Zu`=)cI_waJ-aoi!h0bP;&bzo7Oh8f<1vC*iV(F(|SE@OL z#da+bQJAbwb7+b_USe?MeXoow>nM^rfk?x{+lnWG!sdKrD3sS}@+)`~#y{3|wv7S~ zu;&tL2U?FX-#4Kgut7=fF zLPdHhG)efn7#w`S*CK9%FGT`DFLPkpgU1ldr`XE{Z?$aMy?cB1?8);rHt$7of#hsz zRj9N)XHT8t3TZhfrn)#O9PK1+Vo$*R4PW2yTJVcpp_?a^X^{HRL4$fj2BbWFdJ4ko z31v^;SI$tNDXFQwsgVP!lm7Xqjl%Srqg6;!!hG8#H3QEmJpelx*K-qp&|nL{j)q`D z@3zVJNURX2oMJq_9k)Y=u6i_~L&Sh|Ntt9Q4vBlN2hE3MTy;>x25BX8F+t8JP(DSe zb$n=tSPsPqV8gc8xs&+OGT{Y!Vbb*Jo(Zl3=%N-{Iu1D45({xp)oZgO^QD=CQqbqn zwGD|StXyK>MJW#`Ni6xv$SmR|o}egs`It6C4mWN>(GQ3bH~Jx(QMhc0Faf9Nknwk7 z8F%&bt}2!x{(Q|wkb!&gF-tkZ5K#H#&aC0LRs3;>+ z>`>Marn!S^>Geoh9>sVrVXqR6^p2}!+n-&r<{iyIPgSPQEGo-SODi&)H?Mb91&h9p zS@`1*fjbrE-ZMq;CfNp@Fm8vYxO26qlncMDK=K+`QKESeZH3}qw#vy5^L@RtCL3fU zVX!I03fN6y<9}16!ji#RS($#gbG@PMUgR&iXnj;u>C)7$he~S)SlV^5FtQ9*p5uJ< zYU%xyhYwro&nr~d)isayxpZj=n%6yO05KT8rPz=oG}=CsA-`>|Ma%G8yt4q+q8s6cCD=5RQhHacuQlq}7QmWN5Kg2kQkXqq^AM7svsE`JLs; zl|Jai1J`u;9h<0l?~gRz_XI(>Yk{PN5wtjA|HU#|c8T|;WKL%ZA#UHRF7fc?PN zeU2h62K^!^Ycc-PsKk#&Mc&vfReSZCU{nsWErHV3E_q{)#7tLh#kdyWAgDR2_^&*d z4%AKDaBe3~oVZtbVJi36FQbTA$Qgr_hRN&zFvy!Y*t?EB9xtL+#odE$i-HH|XhXyw zLGk|vu1B=-lJblI`V=sj#J0j9b6jkP;JKZc*ocEiriP1{h9ykqni@SXwm_K)320|k zIFq?E=o0h^Wfa^O1vrUh3yu zl+;K`DFUq&1t1HCqHEatukUTDM4T@b2t;vL8EC;3Py!7CVRnY3AkJX$o}lb*3z{Qr zmI!%lMx^YtvkA5(y8}xuA`2P&BPK)!M;+}wL^XjUbTJEtc%*cE`dV%x1tG zxfRs{$p3Cgn?~XU+O_L~^ZSp*qV41dh{lFjr&VX0g2vD%RWd=@Q8COWwf_ogMQ$L8 zWS2lBZl$yr?#LWmYkvGCgdU)o*eikJn{eJ!<)vDe<;VDOBf{Ky@be}GRFaNH(T88K7k5^kWBq&?A}f|nPHo4} z#v18HCBRaiQ9TiRZsSSk diff --git a/docs/users_guide/ocean/tasks/inertial_gravity_wave.md b/docs/users_guide/ocean/tasks/inertial_gravity_wave.md index 2783f6426..7c27b5992 100644 --- a/docs/users_guide/ocean/tasks/inertial_gravity_wave.md +++ b/docs/users_guide/ocean/tasks/inertial_gravity_wave.md @@ -20,11 +20,11 @@ and the pressure gradient used is the gradient of the sea surface height. Horizontal mixing and bottom friction are also neglected. The nonlinear momentum terms are not included and the layer thickness equation is linearized. -The analysis step computes the root mean-square-error of the difference between +The analysis step computes the L2-norm of the difference between the simulated SSH field and the exact solution at the end of the simulation. It also computes the convergence rate with resolution. -The visualization step produces two plots: the convergence of the RMSE with +The visualization step produces two plots: the convergence of the L2-norm with resolution and a plan-view of the simulated, exact, and (simulated-analytical) SSH fields. @@ -114,9 +114,6 @@ n_wavelengths_y = 2 # Convergence threshold below which the test fails conv_thresh = 1.8 -# Convergence rate above which a warning is issued -conv_max = 2.2 - # time step per resolution (s/km), since dt is proportional to resolution dt_per_km = 3.0 diff --git a/docs/users_guide/ocean/tasks/manufactured_solution.md b/docs/users_guide/ocean/tasks/manufactured_solution.md index 0016e894f..ae97febe0 100644 --- a/docs/users_guide/ocean/tasks/manufactured_solution.md +++ b/docs/users_guide/ocean/tasks/manufactured_solution.md @@ -26,11 +26,11 @@ model is configured without vertical advection and mixing. No tracers are enable and the pressure gradient used is the gradient of the sea surface height. Horizontal mixing and bottom friction are also neglected. -The analysis step computes the root mean-square-error of the difference between +The analysis step computes the L2-norm of the difference between the simulated SSH field and the exact solution at the end of the simulation. It -also computes the convergence rate with resolution +also computes the convergence rate with resolution. -The visualization step produces two plots: the convergence of the RMSE with +The visualization step produces two plots: the convergence of the L2-norm with resolution and a plan-view of the simulated, exact, and (simulated-exact) SSH fields. @@ -111,9 +111,6 @@ dt_per_km = 3.0 # Convergence threshold below which the test fails conv_thresh = 1.8 - -# Convergence rate above which a warning is issued -conv_max = 2.2 ``` ### cores diff --git a/polaris/ocean/tasks/cosine_bell/forward.py b/polaris/ocean/tasks/cosine_bell/forward.py index 638b58134..3c2804a4d 100644 --- a/polaris/ocean/tasks/cosine_bell/forward.py +++ b/polaris/ocean/tasks/cosine_bell/forward.py @@ -25,7 +25,7 @@ def __init__(self, component, name, subdir, resolution, mesh, init): resolution : float The resolution of the (uniform) mesh in km - base_mesh : polaris.Step + mesh : polaris.Step The base mesh step init : polaris.Step diff --git a/polaris/ocean/tasks/inertial_gravity_wave/forward.py b/polaris/ocean/tasks/inertial_gravity_wave/forward.py index 599c22c56..9d50aba11 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/forward.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/forward.py @@ -26,8 +26,11 @@ def __init__(self, component, name, resolution, subdir, init): resolution : km The resolution of the test case in km - taskdir : str - The subdirectory that the task belongs to + subdir : str + The subdirectory that the step belongs to + + init : polaris.Step + The step which generates the mesh and initial condition """ super().__init__(component=component, name=name, subdir=subdir, diff --git a/polaris/ocean/tasks/manufactured_solution/forward.py b/polaris/ocean/tasks/manufactured_solution/forward.py index 3890f2895..cd79f3e52 100644 --- a/polaris/ocean/tasks/manufactured_solution/forward.py +++ b/polaris/ocean/tasks/manufactured_solution/forward.py @@ -29,8 +29,11 @@ def __init__(self, component, name, resolution, subdir, init): resolution : km The resolution of the test case in km - taskdir : str + subdir : str The subdirectory that the task belongs to + + init : polaris.Step + The step which generates the mesh and initial condition """ super().__init__(component=component, name=name, subdir=subdir, From e0c72e2af4f19c2f18434e95bd418a3209c308b4 Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 9 Oct 2023 10:01:22 -0500 Subject: [PATCH 13/17] Change cfg units for convergence tests to hours --- polaris/ocean/convergence/analysis.py | 6 +++--- polaris/ocean/convergence/convergence.cfg | 8 ++++---- polaris/ocean/convergence/forward.py | 7 +++++-- polaris/ocean/tasks/cosine_bell/cosine_bell.cfg | 10 +++++----- polaris/ocean/tasks/cosine_bell/init.py | 5 ++--- polaris/ocean/tasks/cosine_bell/viz.py | 2 +- .../inertial_gravity_wave/inertial_gravity_wave.cfg | 10 +++++----- .../manufactured_solution/manufactured_solution.cfg | 10 +++++----- 8 files changed, 30 insertions(+), 28 deletions(-) diff --git a/polaris/ocean/convergence/analysis.py b/polaris/ocean/convergence/analysis.py index 0821bb2ff..7c8d1f342 100644 --- a/polaris/ocean/convergence/analysis.py +++ b/polaris/ocean/convergence/analysis.py @@ -282,13 +282,13 @@ def compute_error(self, mesh_name, variable_name, zidx=None, config = self.config section = config['convergence'] eval_time = section.getfloat('convergence_eval_time') - s_per_day = 86400.0 + s_per_hour = 3600.0 field_exact = self.exact_solution(mesh_name, variable_name, - time=eval_time * s_per_day, + time=eval_time * s_per_hour, zidx=zidx) field_mpas = self.get_output_field(mesh_name, variable_name, - time=eval_time * s_per_day, + time=eval_time * s_per_hour, zidx=zidx) diff = field_exact - field_mpas diff --git a/polaris/ocean/convergence/convergence.cfg b/polaris/ocean/convergence/convergence.cfg index 3b9e45ad1..e4e202244 100644 --- a/polaris/ocean/convergence/convergence.cfg +++ b/polaris/ocean/convergence/convergence.cfg @@ -1,7 +1,7 @@ [convergence] -# Evaluation time for convergence analysis (in days) -convergence_eval_time = 1.0 +# Evaluation time for convergence analysis (in hours) +convergence_eval_time = 24.0 # Convergence threshold below which a test fails convergence_thresh = 1.0 @@ -25,8 +25,8 @@ split_dt_per_km = 30.0 # since btr_dt is proportional to resolution btr_dt_per_km = 1.5 -# Run duration in days +# Run duration in hours run_duration = ${convergence:convergence_eval_time} -# Output interval in days +# Output interval in hours output_interval = ${run_duration} diff --git a/polaris/ocean/convergence/forward.py b/polaris/ocean/convergence/forward.py index c482952cd..d1754cc85 100644 --- a/polaris/ocean/convergence/forward.py +++ b/polaris/ocean/convergence/forward.py @@ -127,11 +127,14 @@ def dynamic_model_config(self, at_setup): btr_dt_str = get_time_interval_string( seconds=btr_dt_per_km * self.resolution) + s_per_hour = 3600. run_duration = section.getfloat('run_duration') - run_duration_str = get_time_interval_string(days=run_duration) + run_duration_str = get_time_interval_string( + seconds=run_duration * s_per_hour) output_interval = section.getfloat('output_interval') - output_interval_str = get_time_interval_string(days=output_interval) + output_interval_str = get_time_interval_string( + seconds=output_interval * s_per_hour) replacements = dict( time_integrator=time_integrator, diff --git a/polaris/ocean/tasks/cosine_bell/cosine_bell.cfg b/polaris/ocean/tasks/cosine_bell/cosine_bell.cfg index 5486e33a5..fb17734cb 100644 --- a/polaris/ocean/tasks/cosine_bell/cosine_bell.cfg +++ b/polaris/ocean/tasks/cosine_bell/cosine_bell.cfg @@ -23,7 +23,7 @@ min_pc_fraction = 0.1 # config options for convergence tests [convergence] -# Evaluation time for convergence analysis (in days) +# Evaluation time for convergence analysis (in hours) convergence_eval_time = ${cosine_bell:vel_pd} # Convergence threshold below which a test fails @@ -42,10 +42,10 @@ time_integrator = RK4 # RK4 time step per resolution (s/km), since dt is proportional to resolution rk4_dt_per_km = 3.0 -# Run duration in days +# Run duration in hours run_duration = ${cosine_bell:vel_pd} -# Output interval in days +# Output interval in hours output_interval = ${cosine_bell:vel_pd} @@ -70,8 +70,8 @@ radius = 2123666.6667 # hill max of tracer psi0 = 1.0 -# time (days) for bell to transit equator once -vel_pd = 24.0 +# time (hours) for bell to transit equator once +vel_pd = 576. # convergence threshold below which the test fails convergence_thresh = 1.8 diff --git a/polaris/ocean/tasks/cosine_bell/init.py b/polaris/ocean/tasks/cosine_bell/init.py index 5eb41aa80..4c3b26ef9 100644 --- a/polaris/ocean/tasks/cosine_bell/init.py +++ b/polaris/ocean/tasks/cosine_bell/init.py @@ -1,6 +1,5 @@ import numpy as np import xarray as xr -from mpas_tools.cime.constants import constants from mpas_tools.io import write_netcdf from mpas_tools.transects import lon_lat_to_cartesian from mpas_tools.vector import Vector @@ -98,9 +97,9 @@ def run(self): ds['tracer3'] = ds.tracer1 # Initialize velocity - seconds_per_day = constants['SHR_CONST_CDAY'] + s_per_hour = 3600. velocity = (2.0 * np.pi * np.cos(angleEdge) * sphere_radius * - np.cos(latEdge) / (seconds_per_day * vel_pd)) + np.cos(latEdge) / (s_per_hour * vel_pd)) velocity_array, _ = xr.broadcast(velocity, ds.refZMid) ds['normalVelocity'] = velocity_array.expand_dims(dim='Time', axis=0) diff --git a/polaris/ocean/tasks/cosine_bell/viz.py b/polaris/ocean/tasks/cosine_bell/viz.py index 65b313692..c115629ba 100644 --- a/polaris/ocean/tasks/cosine_bell/viz.py +++ b/polaris/ocean/tasks/cosine_bell/viz.py @@ -148,5 +148,5 @@ def run(self): ds_init.lon.values, ds_init.lat.values, ds_out.tracer1.values, out_filename='final.png', config=config, colormap_section='cosine_bell_viz', - title=f'{mesh_name} tracer after {run_duration:g} days', + title=f'{mesh_name} tracer after {run_duration/24.:g} days', plot_land=False) diff --git a/polaris/ocean/tasks/inertial_gravity_wave/inertial_gravity_wave.cfg b/polaris/ocean/tasks/inertial_gravity_wave/inertial_gravity_wave.cfg index 84e8d4534..c33e82230 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/inertial_gravity_wave.cfg +++ b/polaris/ocean/tasks/inertial_gravity_wave/inertial_gravity_wave.cfg @@ -22,8 +22,8 @@ conv_thresh = 1.8 # time step per resolution (s/km), since dt is proportional to resolution dt_per_km = 3.0 -# Run duration in days -run_duration = 0.4166667 +# Run duration in hours +run_duration = 10.0 [vertical_grid] @@ -48,7 +48,7 @@ min_pc_fraction = 0.1 # config options for spherical convergence tests [convergence] -# Evaluation time for convergence analysis (in days) +# Evaluation time for convergence analysis (in hours) convergence_eval_time = ${inertial_gravity_wave:run_duration} # Convergence threshold below which a test fails @@ -67,10 +67,10 @@ time_integrator = RK4 # RK4 time step per resolution (s/km), since dt is proportional to resolution rk4_dt_per_km = ${inertial_gravity_wave:dt_per_km} -# Run duration in days +# Run duration in hours run_duration = ${inertial_gravity_wave:run_duration} -# Output interval in days +# Output interval in hours output_interval = ${inertial_gravity_wave:run_duration} diff --git a/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg b/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg index f241fcb8e..4dc4195b5 100644 --- a/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg +++ b/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg @@ -30,8 +30,8 @@ dt_per_km = 3.0 # Convergence threshold below which the test fails conv_thresh = 1.8 -# Run duration in days -run_duration = 0.4166667 +# Run duration in hours +run_duration = 10.0 [vertical_grid] @@ -56,7 +56,7 @@ min_pc_fraction = 0.1 # config options for spherical convergence tests [convergence] -# Evaluation time for convergence analysis (in days) +# Evaluation time for convergence analysis (in hours) convergence_eval_time = ${manufactured_solution:run_duration} # Convergence threshold below which a test fails @@ -75,8 +75,8 @@ time_integrator = RK4 # RK4 time step per resolution (s/km), since dt is proportional to resolution rk4_dt_per_km = ${manufactured_solution:dt_per_km} -# Run duration in days +# Run duration in hours run_duration = ${manufactured_solution:run_duration} -# Output interval in days +# Output interval in hours output_interval = ${manufactured_solution:run_duration} From fd0a9e53feb2bac47ab931bed672535f59b2775a Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 9 Oct 2023 10:33:07 -0700 Subject: [PATCH 14/17] fixup convergence analysis step --- polaris/ocean/convergence/analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/polaris/ocean/convergence/analysis.py b/polaris/ocean/convergence/analysis.py index 7c8d1f342..2e39fcc45 100644 --- a/polaris/ocean/convergence/analysis.py +++ b/polaris/ocean/convergence/analysis.py @@ -243,7 +243,7 @@ def plot_convergence(self, variable_name, title, zidx): logger.info(f'Order of convergence for {title}: ' f'{conv_round}') - if convergence < conv_thresh: + if conv_round < conv_thresh: logger.error(f'Error: order of convergence for {title}\n' f' {conv_round} < min tolerance ' f'{conv_thresh}') From 8e3ee428a6ce85c0900b30dc5a953c220f981d4d Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 9 Oct 2023 10:52:59 -0700 Subject: [PATCH 15/17] Add convergence suite --- polaris/ocean/suites/convergence.txt | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 polaris/ocean/suites/convergence.txt diff --git a/polaris/ocean/suites/convergence.txt b/polaris/ocean/suites/convergence.txt new file mode 100644 index 000000000..6c01e8ccf --- /dev/null +++ b/polaris/ocean/suites/convergence.txt @@ -0,0 +1,10 @@ +ocean/planar/inertial_gravity_wave +ocean/planar/manufactured_solution +ocean/spherical/icos/cosine_bell + cached: icos_base_mesh_60km icos_init_60km icos_base_mesh_120km icos_init_120km + cached: icos_base_mesh_240km icos_init_240km icos_base_mesh_480km icos_init_480km +ocean/spherical/qu/cosine_bell + cached: qu_base_mesh_60km qu_init_60km qu_base_mesh_90km qu_init_90km + cached: qu_base_mesh_120km qu_init_120km qu_base_mesh_150km qu_init_150km + cached: qu_base_mesh_180km qu_init_180km qu_base_mesh_210km qu_init_210km + cached: qu_base_mesh_240km qu_init_240km From c46720a64cba8486e8f326b8654f4973aa3a312e Mon Sep 17 00:00:00 2001 From: Carolyn Begeman Date: Mon, 9 Oct 2023 13:22:02 -0700 Subject: [PATCH 16/17] Remove convergence plot from viz steps --- .../ocean/tasks/inertial_gravity_wave/viz.py | 22 ------------------ .../ocean/tasks/manufactured_solution/viz.py | 23 +------------------ 2 files changed, 1 insertion(+), 44 deletions(-) diff --git a/polaris/ocean/tasks/inertial_gravity_wave/viz.py b/polaris/ocean/tasks/inertial_gravity_wave/viz.py index 93ea2d8fb..703b2bdc2 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/viz.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/viz.py @@ -53,7 +53,6 @@ def __init__(self, component, resolutions, taskdir): filename=f'output_{mesh_name}.nc', target=f'../forward/{mesh_name}/output.nc') - self.add_output_file('convergence.png') self.add_output_file('comparison.png') def run(self): @@ -117,24 +116,3 @@ def run(self): size='large', ha='right', va='center') fig.savefig('comparison.png', bbox_inches='tight', pad_inches=0.1) - - # Convergence plots - fig = plt.figure() - ax = fig.add_subplot(111) - p = np.polyfit(np.log10(resolutions), np.log10(rmse), 1) - conv = np.round(p[0], 3) - ax.loglog(resolutions, rmse, '-ok', label=f'numerical (order={conv})') - - c = rmse[0] * 1.5 / resolutions[0] - order1 = c * np.power(resolutions, 1) - c = rmse[0] * 1.5 / resolutions[0]**2 - order2 = c * np.power(resolutions, 2) - - ax.loglog(resolutions, order1, '--k', label='first order', alpha=0.3) - ax.loglog(resolutions, order2, 'k', label='second order', alpha=0.3) - ax.set_xlabel('resolution (km)') - ax.set_ylabel('RMS error (m)') - ax.invert_xaxis() - ax.set_title('Error Convergence') - ax.legend(loc='lower left') - fig.savefig('convergence.png', bbox_inches='tight', pad_inches=0.1) diff --git a/polaris/ocean/tasks/manufactured_solution/viz.py b/polaris/ocean/tasks/manufactured_solution/viz.py index 5ebbbf4e0..b485334e7 100644 --- a/polaris/ocean/tasks/manufactured_solution/viz.py +++ b/polaris/ocean/tasks/manufactured_solution/viz.py @@ -53,7 +53,7 @@ def __init__(self, component, resolutions, taskdir): filename=f'output_{mesh_name}.nc', target=f'../forward/{mesh_name}/output.nc') - self.add_output_file('convergence.png') + self.add_output_file('comparison.png') def run(self): """ @@ -114,24 +114,3 @@ def run(self): size='large', ha='right', va='center') fig.savefig('comparison.png', bbox_inches='tight', pad_inches=0.1) - - # Convergence polts - fig = plt.figure() - ax = fig.add_subplot(111) - p = np.polyfit(np.log10(resolutions), np.log10(rmse), 1) - conv = np.round(p[0], 3) - ax.loglog(resolutions, rmse, '-ok', label=f'numerical (order={conv})') - - c = rmse[0] * 1.5 / resolutions[0] - order1 = c * np.power(resolutions, 1) - c = rmse[0] * 1.5 / resolutions[0]**2 - order2 = c * np.power(resolutions, 2) - - ax.loglog(resolutions, order1, '--k', label='first order', alpha=0.3) - ax.loglog(resolutions, order2, 'k', label='second order', alpha=0.3) - ax.set_xlabel('resolution (km)') - ax.set_ylabel('RMS error (m)') - ax.set_title('Error Convergence') - ax.invert_xaxis() - ax.legend(loc='lower left') - fig.savefig('convergence.png', bbox_inches='tight', pad_inches=0.1) From 766432471a33d47482293eabcc6a14e2619b2b0d Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Tue, 10 Oct 2023 07:42:04 -0500 Subject: [PATCH 17/17] Fix days --> hours in cosine bell exact solution --- polaris/ocean/tasks/cosine_bell/analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/polaris/ocean/tasks/cosine_bell/analysis.py b/polaris/ocean/tasks/cosine_bell/analysis.py index 5200e2587..7ba98f8ae 100644 --- a/polaris/ocean/tasks/cosine_bell/analysis.py +++ b/polaris/ocean/tasks/cosine_bell/analysis.py @@ -85,7 +85,7 @@ def exact_solution(self, mesh_name, field_name, time, zidx=None): # distance that the cosine bell center traveled in radians # based on equatorial velocity - distance = 2.0 * np.pi * time / (86400.0 * vel_pd) + distance = 2.0 * np.pi * time / (3600.0 * vel_pd) # new location of blob center lon_new = lon_center + distance