From 5c7a5c99c6b0c632a5523462829a6f7eba575fc4 Mon Sep 17 00:00:00 2001 From: JHopeCollins Date: Fri, 13 Jan 2023 14:45:52 +0000 Subject: [PATCH 01/22] first attempt at allowing indexing in replace_{test,trial}_function --- gusto/labels.py | 155 +++++++++++++++++++++++++----------------------- 1 file changed, 82 insertions(+), 73 deletions(-) diff --git a/gusto/labels.py b/gusto/labels.py index 23fd2d1a7..1466960e2 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -7,7 +7,78 @@ from types import MethodType -def replace_test_function(new_test): +def _replace_dict(old, new, idx, replace_type): + """ + Build a dictionary to pass to the ufl.replace routine + The dictionary matches variables in the old term with those in the new + + Consider cases that old is normal Function or MixedFunction + vs cases of new being Function vs MixedFunction vs tuple + Ideally catch all cases or fail gracefully + """ + + replace_dict = {} + + if type(old.ufl_element()) is MixedElement: + if type(new) == tuple: + assert len(new) == len(old.function_space()) + for k, v in zip(split(old), new): + replace_dict[k] = v + + elif type(new) == ufl.algebra.Sum: + replace_dict[old] = new + + elif isinstance(new, ufl.indexed.Indexed): + if idx is None: + raise ValueError('idx must be specified to replace_{replace_type}' + + ' when {replace_type} is Mixed and new is a single component') + replace_dict[split(old)[idx]] = new + + # Otherwise fail if new is not a function + elif not isinstance(new, type(old)): + raise ValueError(f'new must be a tuple or {type(old)}, not type {type(new)}') + + # Now handle MixedElements separately as these need indexing + elif type(new.ufl_element()) is MixedElement: + assert len(new.function_space()) == len(old.function_space()) + # If idx specified, replace only that component + if idx is not None: + replace_dict[split(old)[idx]] = split(new)[idx] + # Otherwise replace all components + else: + for k, v in zip(split(old), split(new)): + replace_dict[k] = v + + # Otherwise 'new' is a normal Function + else: + if idx is None: + raise ValueError('idx must be specified to replace_{replace_type}' + + ' when {replace_type} is Mixed and new is a single component') + replace_dict[split(old)[idx]] = new + + # old is a normal Function + else: + if type(new) is tuple: + if idx is None: + raise ValueError('idx must be specified to replace_{replace_type}' + + ' when new is a tuple') + replace_dict[old] = new[idx] + elif isinstance(new, ufl.indexed.Indexed): + replace_dict[old] = new + elif not isinstance(new, type(old)): + raise ValueError(f'new must be a {type(old)}, not type {type(new)}') + elif type(new.ufl_element()) == MixedElement: + if idx is None: + raise ValueError('idx must be specified to replace_{replace_type}' + + ' when new is a tuple') + replace_dict[old] = split(new)[idx] + else: + replace_dict[old] = new + + return replace_dict + + +def replace_test_function(new_test, idx=None): """ A routine to replace the test function in a term with a new test function. @@ -30,14 +101,15 @@ def repl(t): Returns: :class:`Term`: the new term. """ - test = t.form.arguments()[0] - new_form = ufl.replace(t.form, {test: new_test}) + old_test = t.form.arguments()[0] + replace_dict = _replace_dict(old_test, new_test, idx, 'test') + new_form = ufl.replace(t.form, replace_dict) return Term(new_form, t.labels) return repl -def replace_trial_function(new): +def replace_trial_function(new_trial, idx=None): """ A routine to replace the trial function in a term with a new expression. @@ -65,14 +137,15 @@ def repl(t): """ if len(t.form.arguments()) != 2: raise TypeError('Trying to replace trial function of a form that is not linear') - trial = t.form.arguments()[1] - new_form = ufl.replace(t.form, {trial: new}) + old_trial = t.form.arguments()[1] + replace_dict = _replace_dict(old_trial, new_trial, idx, 'trial') + new_form = ufl.replace(t.form, replace_dict) return Term(new_form, t.labels) return repl -def replace_subject(new, idx=None): +def replace_subject(new_subj, idx=None): """ A routine to replace the subject in a term with a new variable. @@ -97,73 +170,9 @@ def repl(t): :class:`Term`: the new term. """ - subj = t.get(subject) - - # Build a dictionary to pass to the ufl.replace routine - # The dictionary matches variables in the old term with those in the new - replace_dict = {} - - # Consider cases that subj is normal Function or MixedFunction - # vs cases of new being Function vs MixedFunction vs tuple - # Ideally catch all cases or fail gracefully - if type(subj.ufl_element()) is MixedElement: - if type(new) == tuple: - assert len(new) == len(subj.function_space()) - for k, v in zip(split(subj), new): - replace_dict[k] = v - - elif type(new) == ufl.algebra.Sum: - replace_dict[subj] = new - - elif isinstance(new, ufl.indexed.Indexed): - if idx is None: - raise ValueError('idx must be specified to replace_subject' - + ' when subject is Mixed and new is a single component') - replace_dict[split(subj)[idx]] = new - - # Otherwise fail if new is not a function - elif not isinstance(new, Function): - raise ValueError(f'new must be a tuple or Function, not type {type(new)}') - - # Now handle MixedElements separately as these need indexing - elif type(new.ufl_element()) is MixedElement: - assert len(new.function_space()) == len(subj.function_space()) - # If idx specified, replace only that component - if idx is not None: - replace_dict[split(subj)[idx]] = split(new)[idx] - # Otherwise replace all components - else: - for k, v in zip(split(subj), split(new)): - replace_dict[k] = v - - # Otherwise 'new' is a normal Function - else: - if idx is None: - raise ValueError('idx must be specified to replace_subject' - + ' when subject is Mixed and new is a single component') - replace_dict[split(subj)[idx]] = new - - # subj is a normal Function - else: - if type(new) is tuple: - if idx is None: - raise ValueError('idx must be specified to replace_subject' - + ' when new is a tuple') - replace_dict[subj] = new[idx] - elif isinstance(new, ufl.indexed.Indexed): - replace_dict[subj] = new - elif not isinstance(new, Function): - raise ValueError(f'new must be a Function, not type {type(new)}') - elif type(new.ufl_element()) == MixedElement: - if idx is None: - raise ValueError('idx must be specified to replace_subject' - + ' when new is a tuple') - replace_dict[subj] = split(new)[idx] - else: - replace_dict[subj] = new - + old_subj = t.get(subject) + replace_dict = _replace_dict(old_subj, new_subj, idx, 'subject') new_form = ufl.replace(t.form, replace_dict) - return Term(new_form, t.labels) return repl From 4613dc003c4482bba5dc884e259f7f521ffa4d50 Mon Sep 17 00:00:00 2001 From: JHopeCollins Date: Mon, 16 Jan 2023 16:19:40 +0000 Subject: [PATCH 02/22] replace_*: allow ufl.algebra.Sum to be replacement for non-mixed functionspaces --- gusto/labels.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/gusto/labels.py b/gusto/labels.py index 1466960e2..5a854d1cd 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -42,8 +42,10 @@ def _replace_dict(old, new, idx, replace_type): elif type(new.ufl_element()) is MixedElement: assert len(new.function_space()) == len(old.function_space()) # If idx specified, replace only that component + if idx is not None: replace_dict[split(old)[idx]] = split(new)[idx] + # Otherwise replace all components else: for k, v in zip(split(old), split(new)): @@ -63,15 +65,22 @@ def _replace_dict(old, new, idx, replace_type): raise ValueError('idx must be specified to replace_{replace_type}' + ' when new is a tuple') replace_dict[old] = new[idx] + + elif type(new) == ufl.algebra.Sum: + replace_dict[old] = new + elif isinstance(new, ufl.indexed.Indexed): replace_dict[old] = new + elif not isinstance(new, type(old)): raise ValueError(f'new must be a {type(old)}, not type {type(new)}') + elif type(new.ufl_element()) == MixedElement: if idx is None: raise ValueError('idx must be specified to replace_{replace_type}' + ' when new is a tuple') replace_dict[old] = split(new)[idx] + else: replace_dict[old] = new From c32a7740a28a609b546250e6e786c8cd247b4a00 Mon Sep 17 00:00:00 2001 From: JHopeCollins Date: Mon, 16 Jan 2023 16:59:42 +0000 Subject: [PATCH 03/22] replace_* label maps: check new function is of acceptable type earlier --- gusto/labels.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/gusto/labels.py b/gusto/labels.py index 5a854d1cd..0d934b31d 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -19,12 +19,20 @@ def _replace_dict(old, new, idx, replace_type): replace_dict = {} + acceptable_types = (type(old), ufl.algebra.Sum, ufl.indexed.Indexed) + if replace_type == 'trial': + acceptable_types = (*acceptable_types, Function) + if type(old.ufl_element()) is MixedElement: if type(new) == tuple: assert len(new) == len(old.function_space()) for k, v in zip(split(old), new): replace_dict[k] = v + # Otherwise fail if new is not a function + elif not isinstance(new, acceptable_types): + raise ValueError(f'new must be a tuple or {type(old)}, not type {type(new)}') + elif type(new) == ufl.algebra.Sum: replace_dict[old] = new @@ -34,10 +42,6 @@ def _replace_dict(old, new, idx, replace_type): + ' when {replace_type} is Mixed and new is a single component') replace_dict[split(old)[idx]] = new - # Otherwise fail if new is not a function - elif not isinstance(new, type(old)): - raise ValueError(f'new must be a tuple or {type(old)}, not type {type(new)}') - # Now handle MixedElements separately as these need indexing elif type(new.ufl_element()) is MixedElement: assert len(new.function_space()) == len(old.function_space()) @@ -66,15 +70,15 @@ def _replace_dict(old, new, idx, replace_type): + ' when new is a tuple') replace_dict[old] = new[idx] + elif not isinstance(new, acceptable_types): + raise ValueError(f'new must be a {type(old)}, not type {type(new)}') + elif type(new) == ufl.algebra.Sum: replace_dict[old] = new elif isinstance(new, ufl.indexed.Indexed): replace_dict[old] = new - elif not isinstance(new, type(old)): - raise ValueError(f'new must be a {type(old)}, not type {type(new)}') - elif type(new.ufl_element()) == MixedElement: if idx is None: raise ValueError('idx must be specified to replace_{replace_type}' From e6f7ab04b42cd8b9b28c8ff6bf0ab42332adf3a3 Mon Sep 17 00:00:00 2001 From: JHopeCollins Date: Mon, 16 Jan 2023 17:23:25 +0000 Subject: [PATCH 04/22] replace_* label maps: better error messages --- gusto/labels.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/gusto/labels.py b/gusto/labels.py index 0d934b31d..7a16571a3 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -23,6 +23,8 @@ def _replace_dict(old, new, idx, replace_type): if replace_type == 'trial': acceptable_types = (*acceptable_types, Function) + type_error_message = f'new must be a {tuple} or '+' or '.join((f"{t}" for t in acceptable_types))+f', not {type(new)}' + if type(old.ufl_element()) is MixedElement: if type(new) == tuple: assert len(new) == len(old.function_space()) @@ -31,7 +33,7 @@ def _replace_dict(old, new, idx, replace_type): # Otherwise fail if new is not a function elif not isinstance(new, acceptable_types): - raise ValueError(f'new must be a tuple or {type(old)}, not type {type(new)}') + raise TypeError(type_error_message) elif type(new) == ufl.algebra.Sum: replace_dict[old] = new @@ -67,11 +69,11 @@ def _replace_dict(old, new, idx, replace_type): if type(new) is tuple: if idx is None: raise ValueError('idx must be specified to replace_{replace_type}' - + ' when new is a tuple') + + ' when new is a tuple and {replace_type} is not Mixed') replace_dict[old] = new[idx] elif not isinstance(new, acceptable_types): - raise ValueError(f'new must be a {type(old)}, not type {type(new)}') + raise TypeError(type_error_message) elif type(new) == ufl.algebra.Sum: replace_dict[old] = new @@ -82,7 +84,7 @@ def _replace_dict(old, new, idx, replace_type): elif type(new.ufl_element()) == MixedElement: if idx is None: raise ValueError('idx must be specified to replace_{replace_type}' - + ' when new is a tuple') + + ' when new is a tuple and {replace_type} is not Mixed') replace_dict[old] = split(new)[idx] else: From 17b1ec9019453d0ec4a1fddda19015c102bb389c Mon Sep 17 00:00:00 2001 From: JHopeCollins Date: Tue, 24 Jan 2023 10:04:43 +0000 Subject: [PATCH 05/22] replace_* label maps: leave type checking to ufl.replace when building replace dictionary --- gusto/labels.py | 114 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 113 insertions(+), 1 deletion(-) diff --git a/gusto/labels.py b/gusto/labels.py index 7a16571a3..2ce07ff76 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -7,7 +7,7 @@ from types import MethodType -def _replace_dict(old, new, idx, replace_type): +def _replace_dict_old(old, new, idx, replace_type): """ Build a dictionary to pass to the ufl.replace routine The dictionary matches variables in the old term with those in the new @@ -93,6 +93,118 @@ def _replace_dict(old, new, idx, replace_type): return replace_dict +def _replace_dict_dumb1(old, new, idx, replace_type): + """ + Build a dictionary to pass to the ufl.replace routine + The dictionary matches variables in the old term with those in the new + + Does not check types unless indexing is required (leave type-checking to ufl) + """ + + replace_dict = {} + + if type(old.ufl_element()) is MixedElement: + + if type(new) is tuple: + if len(new) != len(old.function_space()): + raise ValueError(f"tuple of new {replace_type} must be same length as replaced mixed {replace_type} of type {old}") + if idx is None: + for k, v in zip(split(old), new): + replace_dict[k] = v + else: + replace_dict[split(old)[idx]] = new[idx] + + elif type(new.ufl_element()) is MixedElement: + if len(new.function_space()) != len(old.function_space()): + raise ValueError(f"New mixed {replace_type} of type {new} must be same length as replaced mixed {replace_type} of type {old}") + if idx is None: + for k, v in zip(split(old), split(new)): + replace_dict[k] = v + else: + replace_dict[split(old)[idx]] = split(new)[idx] + + else: # new is not indexable + if idx is None: + raise ValueError(f"idx must be specified to replace_{replace_type}" + + f" when replaced {replace_type} of type {old} is mixed and new {replace_type} of type {new} is a single component") + replace_dict[split(old)[idx]] = new + + else: # old is not mixed + + if type(new) is tuple: + if idx is None: + raise ValueError(f"idx must be specified to replace_{replace_type}" + + f" when replaced {replace_type} of type {old} is not mixed and new {replace_type} is a tuple") + replace_dict[old] = new[idx] + + elif type(new.ufl_element()) is MixedElement: + if idx is None: + raise ValueError(f"idx must be specified to replace_{replace_type}" + + f" when replaced {replace_type} of type {old} is not mixed and new {replace_type} of type {new} is mixed") + replace_dict[old] = split(new)[idx] + + else: # new is not mixed + replace_dict[old] = new + + return replace_dict + + +def _replace_dict(old, new, idx, replace_type): + """ + Build a dictionary to pass to the ufl.replace routine + The dictionary matches variables in the old term with those in the new + + Does not check types unless indexing is required (leave type-checking to ufl) + """ + + replace_dict = {} + + if type(old.ufl_element()) is MixedElement: + + mixed_new = (hasattr(new, "ufl_element") and type(new.ufl_element()) is MixedElement) + indexable_new = type(new) is tuple or mixed_new + + if indexable_new: + new_len = len(new) if type(new) is tuple else len(new.function_space()) + split_new = new if type(new) is tuple else split(new) + + if new_len != len(old.function_space()): + raise ValueError(f"new {replace_type} of type {new} must be same length as replaced mixed {replace_type} of type {old}") + + if idx is None: + for k, v in zip(split(old), split_new): + replace_dict[k] = v + else: + replace_dict[split(old)[idx]] = split_new[idx] + + else: # new is not indexable + if idx is None: + raise ValueError(f"idx must be specified to replace_{replace_type}" + + f" when replaced {replace_type} of type {old} is mixed and new {replace_type} of type {new} is a single component") + + replace_dict[split(old)[idx]] = new + + else: # old is not mixed + + # or will short-circuit + mixed_new = (hasattr(new, "ufl_element") and type(new.ufl_element()) is MixedElement) + indexable_new = type(new) is tuple or mixed_new + + if indexable_new: + split_new = new if type(new) is tuple else split(new) + + if idx is None: + raise ValueError(f"idx must be specified to replace_{replace_type}" + + f" when replaced {replace_type} of type {old} is not mixed and new {replace_type} is of indexable type {new}") + + replace_dict[old] = split_new[idx] + + else: + replace_dict[old] = new + + return replace_dict + + def replace_test_function(new_test, idx=None): """ A routine to replace the test function in a term with a new test function. From ca1ed0c6f8a3aad488aecd16ef91e3153ad64f13 Mon Sep 17 00:00:00 2001 From: JHopeCollins Date: Tue, 24 Jan 2023 11:02:40 +0000 Subject: [PATCH 06/22] replace_* label maps: tidy up error messages and test skipping --- gusto/labels.py | 18 ++++++++++-------- unit-tests/fml_tests/test_replace_subject.py | 8 ++++---- 2 files changed, 14 insertions(+), 12 deletions(-) diff --git a/gusto/labels.py b/gusto/labels.py index 2ce07ff76..47ec5bc46 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -161,7 +161,7 @@ def _replace_dict(old, new, idx, replace_type): if type(old.ufl_element()) is MixedElement: - mixed_new = (hasattr(new, "ufl_element") and type(new.ufl_element()) is MixedElement) + mixed_new = hasattr(new, "ufl_element") and type(new.ufl_element()) is MixedElement indexable_new = type(new) is tuple or mixed_new if indexable_new: @@ -169,7 +169,8 @@ def _replace_dict(old, new, idx, replace_type): split_new = new if type(new) is tuple else split(new) if new_len != len(old.function_space()): - raise ValueError(f"new {replace_type} of type {new} must be same length as replaced mixed {replace_type} of type {old}") + raise ValueError(f"new {replace_type} of type {new} must be same length" + + f"as replaced mixed {replace_type} of type {old}") if idx is None: for k, v in zip(split(old), split_new): @@ -179,23 +180,24 @@ def _replace_dict(old, new, idx, replace_type): else: # new is not indexable if idx is None: - raise ValueError(f"idx must be specified to replace_{replace_type}" - + f" when replaced {replace_type} of type {old} is mixed and new {replace_type} of type {new} is a single component") + raise ValueError(f"idx must be specified to replace_{replace_type} when" + + f" replaced {replace_type} of type {old} is mixed and" + + f" new {replace_type} of type {new} is a single component") replace_dict[split(old)[idx]] = new else: # old is not mixed - # or will short-circuit - mixed_new = (hasattr(new, "ufl_element") and type(new.ufl_element()) is MixedElement) + mixed_new = hasattr(new, "ufl_element") and type(new.ufl_element()) is MixedElement indexable_new = type(new) is tuple or mixed_new if indexable_new: split_new = new if type(new) is tuple else split(new) if idx is None: - raise ValueError(f"idx must be specified to replace_{replace_type}" - + f" when replaced {replace_type} of type {old} is not mixed and new {replace_type} is of indexable type {new}") + raise ValueError(f"idx must be specified to replace_{replace_type} when" + + f" replaced {replace_type} of type {old} is not mixed" + + f" and new {replace_type} of type {new} is indexable") replace_dict[old] = split_new[idx] diff --git a/unit-tests/fml_tests/test_replace_subject.py b/unit-tests/fml_tests/test_replace_subject.py index 1756c4860..8da0ca28e 100644 --- a/unit-tests/fml_tests/test_replace_subject.py +++ b/unit-tests/fml_tests/test_replace_subject.py @@ -20,15 +20,15 @@ def test_replace_subject(subject_type, replacement_type, function_or_indexed): # ------------------------------------------------------------------------ # if subject_type == 'vector' and replacement_type != 'vector': - return + pytest.skip("invalid option combination") elif replacement_type == 'vector' and subject_type != 'vector': - return + pytest.skip("invalid option combination") if replacement_type == 'mixed-component': if subject_type != 'mixed': - return + pytest.skip("invalid option combination") elif function_or_indexed != 'indexed': - return + pytest.skip("invalid option combination") # ------------------------------------------------------------------------ # # Set up From b8a33be39e57cf37df64b613179869b04ffdf869 Mon Sep 17 00:00:00 2001 From: JHopeCollins Date: Tue, 24 Jan 2023 12:07:08 +0000 Subject: [PATCH 07/22] replace_* label maps: add more information to ufl.replace exception message --- gusto/labels.py | 29 +++++++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/gusto/labels.py b/gusto/labels.py index 47ec5bc46..4d5be32ab 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -154,7 +154,7 @@ def _replace_dict(old, new, idx, replace_type): Build a dictionary to pass to the ufl.replace routine The dictionary matches variables in the old term with those in the new - Does not check types unless indexing is required (leave type-checking to ufl) + Does not check types unless indexing is required (leave type-checking to ufl.replace) """ replace_dict = {} @@ -232,7 +232,14 @@ def repl(t): """ old_test = t.form.arguments()[0] replace_dict = _replace_dict(old_test, new_test, idx, 'test') - new_form = ufl.replace(t.form, replace_dict) + + try: + new_form = ufl.replace(t.form, replace_dict) + except Exception as err: + error_message = f"{err} raised by ufl.replace when trying to" \ + + f" replace_test_function with {new_test}" + raise type(err)(error_message) from err + return Term(new_form, t.labels) return repl @@ -268,7 +275,14 @@ def repl(t): raise TypeError('Trying to replace trial function of a form that is not linear') old_trial = t.form.arguments()[1] replace_dict = _replace_dict(old_trial, new_trial, idx, 'trial') - new_form = ufl.replace(t.form, replace_dict) + + try: + new_form = ufl.replace(t.form, replace_dict) + except Exception as err: + error_message = f"{err} raised by ufl.replace when trying to" \ + + f" replace_trial_function with {new_trial}" + raise type(err)(error_message) from err + return Term(new_form, t.labels) return repl @@ -301,7 +315,14 @@ def repl(t): old_subj = t.get(subject) replace_dict = _replace_dict(old_subj, new_subj, idx, 'subject') - new_form = ufl.replace(t.form, replace_dict) + + try: + new_form = ufl.replace(t.form, replace_dict) + except Exception as err: + error_message = f"{err} raised by ufl.replace when trying to" \ + + f" replace_subject with {new_subj}" + raise type(err)(error_message) from err + return Term(new_form, t.labels) return repl From 788dbb7f9fe18c67ea731c8514d35a79d040999c Mon Sep 17 00:00:00 2001 From: JHopeCollins Date: Tue, 24 Jan 2023 13:40:35 +0000 Subject: [PATCH 08/22] replace_* label maps: remove old _replace_dict impl --- gusto/labels.py | 62 +++---------------------------------------------- 1 file changed, 3 insertions(+), 59 deletions(-) diff --git a/gusto/labels.py b/gusto/labels.py index 4d5be32ab..2cd77ddd4 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -93,62 +93,6 @@ def _replace_dict_old(old, new, idx, replace_type): return replace_dict -def _replace_dict_dumb1(old, new, idx, replace_type): - """ - Build a dictionary to pass to the ufl.replace routine - The dictionary matches variables in the old term with those in the new - - Does not check types unless indexing is required (leave type-checking to ufl) - """ - - replace_dict = {} - - if type(old.ufl_element()) is MixedElement: - - if type(new) is tuple: - if len(new) != len(old.function_space()): - raise ValueError(f"tuple of new {replace_type} must be same length as replaced mixed {replace_type} of type {old}") - if idx is None: - for k, v in zip(split(old), new): - replace_dict[k] = v - else: - replace_dict[split(old)[idx]] = new[idx] - - elif type(new.ufl_element()) is MixedElement: - if len(new.function_space()) != len(old.function_space()): - raise ValueError(f"New mixed {replace_type} of type {new} must be same length as replaced mixed {replace_type} of type {old}") - if idx is None: - for k, v in zip(split(old), split(new)): - replace_dict[k] = v - else: - replace_dict[split(old)[idx]] = split(new)[idx] - - else: # new is not indexable - if idx is None: - raise ValueError(f"idx must be specified to replace_{replace_type}" - + f" when replaced {replace_type} of type {old} is mixed and new {replace_type} of type {new} is a single component") - replace_dict[split(old)[idx]] = new - - else: # old is not mixed - - if type(new) is tuple: - if idx is None: - raise ValueError(f"idx must be specified to replace_{replace_type}" - + f" when replaced {replace_type} of type {old} is not mixed and new {replace_type} is a tuple") - replace_dict[old] = new[idx] - - elif type(new.ufl_element()) is MixedElement: - if idx is None: - raise ValueError(f"idx must be specified to replace_{replace_type}" - + f" when replaced {replace_type} of type {old} is not mixed and new {replace_type} of type {new} is mixed") - replace_dict[old] = split(new)[idx] - - else: # new is not mixed - replace_dict[old] = new - - return replace_dict - - def _replace_dict(old, new, idx, replace_type): """ Build a dictionary to pass to the ufl.replace routine @@ -236,7 +180,7 @@ def repl(t): try: new_form = ufl.replace(t.form, replace_dict) except Exception as err: - error_message = f"{err} raised by ufl.replace when trying to" \ + error_message = f"{type(err)} raised by ufl.replace when trying to" \ + f" replace_test_function with {new_test}" raise type(err)(error_message) from err @@ -279,7 +223,7 @@ def repl(t): try: new_form = ufl.replace(t.form, replace_dict) except Exception as err: - error_message = f"{err} raised by ufl.replace when trying to" \ + error_message = f"{type(err)} raised by ufl.replace when trying to" \ + f" replace_trial_function with {new_trial}" raise type(err)(error_message) from err @@ -319,7 +263,7 @@ def repl(t): try: new_form = ufl.replace(t.form, replace_dict) except Exception as err: - error_message = f"{err} raised by ufl.replace when trying to" \ + error_message = f"{type(err)} raised by ufl.replace when trying to" \ + f" replace_subject with {new_subj}" raise type(err)(error_message) from err From 259f2d1e9286569725d618620cdf984e7a62300a Mon Sep 17 00:00:00 2001 From: JHopeCollins Date: Tue, 24 Jan 2023 18:17:04 +0000 Subject: [PATCH 09/22] replace_* label maps: parametrize tests for replace_{subject,trial,test} --- gusto/labels.py | 3 +- unit-tests/fml_tests/test_replace_subject.py | 34 ++++++++++++++++---- 2 files changed, 28 insertions(+), 9 deletions(-) diff --git a/gusto/labels.py b/gusto/labels.py index 2cd77ddd4..a4248f6b0 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -109,10 +109,9 @@ def _replace_dict(old, new, idx, replace_type): indexable_new = type(new) is tuple or mixed_new if indexable_new: - new_len = len(new) if type(new) is tuple else len(new.function_space()) split_new = new if type(new) is tuple else split(new) - if new_len != len(old.function_space()): + if len(split_new) != len(old.function_space()): raise ValueError(f"new {replace_type} of type {new} must be same length" + f"as replaced mixed {replace_type} of type {old}") diff --git a/unit-tests/fml_tests/test_replace_subject.py b/unit-tests/fml_tests/test_replace_subject.py index 8da0ca28e..90bb45a77 100644 --- a/unit-tests/fml_tests/test_replace_subject.py +++ b/unit-tests/fml_tests/test_replace_subject.py @@ -4,16 +4,23 @@ from firedrake import (UnitSquareMesh, FunctionSpace, Function, TestFunction, VectorFunctionSpace, MixedFunctionSpace, dx, inner, - TrialFunctions, split) + TrialFunctions, TrialFunction, split) from gusto.fml import Label -from gusto import subject, replace_subject +from gusto import subject, replace_subject, replace_test_function, replace_trial_function import pytest +replace_funcs = [ + pytest.param((Function, replace_subject), id="replace_subj"), + pytest.param((TestFunction, replace_test_function), id="replace_test"), + pytest.param((TrialFunction, replace_trial_function), id="replace_trial") +] + @pytest.mark.parametrize('subject_type', ['normal', 'mixed', 'vector']) @pytest.mark.parametrize('replacement_type', ['normal', 'mixed', 'mixed-component', 'vector', 'tuple']) @pytest.mark.parametrize('function_or_indexed', ['function', 'indexed']) -def test_replace_subject(subject_type, replacement_type, function_or_indexed): +@pytest.mark.parametrize('replace_func', replace_funcs) +def test_replace_subject(subject_type, replacement_type, function_or_indexed, replace_func): # ------------------------------------------------------------------------ # # Only certain combinations of options are valid @@ -52,6 +59,9 @@ def test_replace_subject(subject_type, replacement_type, function_or_indexed): # Choose subject # ------------------------------------------------------------------------ # + FunctionType = replace_func[0] + replace_map = replace_func[1] + if subject_type == 'normal': V = V0 elif subject_type == 'mixed': @@ -64,7 +74,12 @@ def test_replace_subject(subject_type, replacement_type, function_or_indexed): raise ValueError the_subject = Function(V) - not_subject = Function(V) + + if replace_map is replace_trial_function: + not_subject = TrialFunction(V) + else: + not_subject = Function(V) + test = TestFunction(V) form_1 = inner(the_subject, test)*dx @@ -94,7 +109,7 @@ def test_replace_subject(subject_type, replacement_type, function_or_indexed): else: raise ValueError - the_replacement = Function(V) + the_replacement = FunctionType(V) if function_or_indexed == 'indexed' and replacement_type != 'vector': the_replacement = split(the_replacement) @@ -111,7 +126,12 @@ def test_replace_subject(subject_type, replacement_type, function_or_indexed): # Test replace_subject # ------------------------------------------------------------------------ # + if replace_map is replace_trial_function: + match_label = bar_label + else: + match_label = subject + labelled_form = labelled_form.label_map( - lambda t: t.has_label(subject), - map_if_true=replace_subject(the_replacement, idx=idx) + lambda t: t.has_label(match_label), + map_if_true=replace_map(the_replacement, idx=idx) ) From 1d7852e02a98c44a9fe2486c3ea723ed63b97e7b Mon Sep 17 00:00:00 2001 From: JHopeCollins Date: Tue, 24 Jan 2023 18:33:10 +0000 Subject: [PATCH 10/22] replace_* label maps: remove mixed-component test parameter for in-test if-statement --- unit-tests/fml_tests/test_replace_subject.py | 42 +++++++++----------- 1 file changed, 19 insertions(+), 23 deletions(-) diff --git a/unit-tests/fml_tests/test_replace_subject.py b/unit-tests/fml_tests/test_replace_subject.py index 90bb45a77..374913648 100644 --- a/unit-tests/fml_tests/test_replace_subject.py +++ b/unit-tests/fml_tests/test_replace_subject.py @@ -17,7 +17,7 @@ @pytest.mark.parametrize('subject_type', ['normal', 'mixed', 'vector']) -@pytest.mark.parametrize('replacement_type', ['normal', 'mixed', 'mixed-component', 'vector', 'tuple']) +@pytest.mark.parametrize('replacement_type', ['normal', 'mixed', 'vector', 'tuple']) @pytest.mark.parametrize('function_or_indexed', ['function', 'indexed']) @pytest.mark.parametrize('replace_func', replace_funcs) def test_replace_subject(subject_type, replacement_type, function_or_indexed, replace_func): @@ -26,16 +26,9 @@ def test_replace_subject(subject_type, replacement_type, function_or_indexed, re # Only certain combinations of options are valid # ------------------------------------------------------------------------ # - if subject_type == 'vector' and replacement_type != 'vector': + # only makes sense to replace a vector with a vector + if (subject_type == 'vector') ^ (replacement_type == 'vector'): pytest.skip("invalid option combination") - elif replacement_type == 'vector' and subject_type != 'vector': - pytest.skip("invalid option combination") - - if replacement_type == 'mixed-component': - if subject_type != 'mixed': - pytest.skip("invalid option combination") - elif function_or_indexed != 'indexed': - pytest.skip("invalid option combination") # ------------------------------------------------------------------------ # # Set up @@ -59,9 +52,6 @@ def test_replace_subject(subject_type, replacement_type, function_or_indexed, re # Choose subject # ------------------------------------------------------------------------ # - FunctionType = replace_func[0] - replace_map = replace_func[1] - if subject_type == 'normal': V = V0 elif subject_type == 'mixed': @@ -74,12 +64,7 @@ def test_replace_subject(subject_type, replacement_type, function_or_indexed, re raise ValueError the_subject = Function(V) - - if replace_map is replace_trial_function: - not_subject = TrialFunction(V) - else: - not_subject = Function(V) - + not_subject = TrialFunction(V) test = TestFunction(V) form_1 = inner(the_subject, test)*dx @@ -99,9 +84,6 @@ def test_replace_subject(subject_type, replacement_type, function_or_indexed, re V = Vmixed if subject_type != 'mixed': idx = 0 - elif replacement_type == 'mixed-component': - V = Vmixed - idx = 0 elif replacement_type == 'vector': V = V2 elif replacement_type == 'tuple': @@ -109,12 +91,14 @@ def test_replace_subject(subject_type, replacement_type, function_or_indexed, re else: raise ValueError + FunctionType = replace_func[0] + the_replacement = FunctionType(V) if function_or_indexed == 'indexed' and replacement_type != 'vector': the_replacement = split(the_replacement) - if len(the_replacement) == 1 or replacement_type == 'mixed-component': + if len(the_replacement) == 1: the_replacement = the_replacement[0] if replacement_type == 'tuple': @@ -126,6 +110,8 @@ def test_replace_subject(subject_type, replacement_type, function_or_indexed, re # Test replace_subject # ------------------------------------------------------------------------ # + replace_map = replace_func[1] + if replace_map is replace_trial_function: match_label = bar_label else: @@ -135,3 +121,13 @@ def test_replace_subject(subject_type, replacement_type, function_or_indexed, re lambda t: t.has_label(match_label), map_if_true=replace_map(the_replacement, idx=idx) ) + + # also test indexed + if subject_type == 'mixed' and function_or_indexed == 'indexed': + idx = 0 + the_replacement = split(FunctionType(Vmixed))[idx] + + labelled_form = labelled_form.label_map( + lambda t: t.has_label(match_label), + map_if_true=replace_map(the_replacement, idx=idx) + ) From 648effa617633a043b58ca5a95f1e9b14d05e375 Mon Sep 17 00:00:00 2001 From: JHopeCollins Date: Wed, 25 Jan 2023 09:30:32 +0000 Subject: [PATCH 11/22] replace_* label maps: remove old dictionary builder --- gusto/labels.py | 86 ------------------------------------------------- 1 file changed, 86 deletions(-) diff --git a/gusto/labels.py b/gusto/labels.py index a4248f6b0..331fa32cf 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -7,92 +7,6 @@ from types import MethodType -def _replace_dict_old(old, new, idx, replace_type): - """ - Build a dictionary to pass to the ufl.replace routine - The dictionary matches variables in the old term with those in the new - - Consider cases that old is normal Function or MixedFunction - vs cases of new being Function vs MixedFunction vs tuple - Ideally catch all cases or fail gracefully - """ - - replace_dict = {} - - acceptable_types = (type(old), ufl.algebra.Sum, ufl.indexed.Indexed) - if replace_type == 'trial': - acceptable_types = (*acceptable_types, Function) - - type_error_message = f'new must be a {tuple} or '+' or '.join((f"{t}" for t in acceptable_types))+f', not {type(new)}' - - if type(old.ufl_element()) is MixedElement: - if type(new) == tuple: - assert len(new) == len(old.function_space()) - for k, v in zip(split(old), new): - replace_dict[k] = v - - # Otherwise fail if new is not a function - elif not isinstance(new, acceptable_types): - raise TypeError(type_error_message) - - elif type(new) == ufl.algebra.Sum: - replace_dict[old] = new - - elif isinstance(new, ufl.indexed.Indexed): - if idx is None: - raise ValueError('idx must be specified to replace_{replace_type}' - + ' when {replace_type} is Mixed and new is a single component') - replace_dict[split(old)[idx]] = new - - # Now handle MixedElements separately as these need indexing - elif type(new.ufl_element()) is MixedElement: - assert len(new.function_space()) == len(old.function_space()) - # If idx specified, replace only that component - - if idx is not None: - replace_dict[split(old)[idx]] = split(new)[idx] - - # Otherwise replace all components - else: - for k, v in zip(split(old), split(new)): - replace_dict[k] = v - - # Otherwise 'new' is a normal Function - else: - if idx is None: - raise ValueError('idx must be specified to replace_{replace_type}' - + ' when {replace_type} is Mixed and new is a single component') - replace_dict[split(old)[idx]] = new - - # old is a normal Function - else: - if type(new) is tuple: - if idx is None: - raise ValueError('idx must be specified to replace_{replace_type}' - + ' when new is a tuple and {replace_type} is not Mixed') - replace_dict[old] = new[idx] - - elif not isinstance(new, acceptable_types): - raise TypeError(type_error_message) - - elif type(new) == ufl.algebra.Sum: - replace_dict[old] = new - - elif isinstance(new, ufl.indexed.Indexed): - replace_dict[old] = new - - elif type(new.ufl_element()) == MixedElement: - if idx is None: - raise ValueError('idx must be specified to replace_{replace_type}' - + ' when new is a tuple and {replace_type} is not Mixed') - replace_dict[old] = split(new)[idx] - - else: - replace_dict[old] = new - - return replace_dict - - def _replace_dict(old, new, idx, replace_type): """ Build a dictionary to pass to the ufl.replace routine From c6a882207abdc04e78652996de7cc510a92b41d6 Mon Sep 17 00:00:00 2001 From: jshipton Date: Thu, 26 Jan 2023 22:53:05 +0000 Subject: [PATCH 12/22] hacky fix --- gusto/labels.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/gusto/labels.py b/gusto/labels.py index 68b769460..535dd600f 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -185,7 +185,14 @@ def repl(t): # the subject if t.has_label(perp): perp_function = t.get(perp) - new_form = ufl.replace(new_form, {split(new_subj)[0]: perp_function(split(new_subj)[0])}) + mixed_new = hasattr(new_subj, "ufl_element") and type(new_subj.ufl_element()) is MixedElement + if idx == 0: + new_form = ufl.replace(new_form, + {new_subj: perp_function(new_subj)}) + elif mixed_new: + new_form = ufl.replace( + new_form,{split(new_subj)[0]: + perp_function(split(new_subj)[0])}) return Term(new_form, t.labels) From 829488c04f646ec6a19d5df7adb3c25ab1e177f7 Mon Sep 17 00:00:00 2001 From: jshipton Date: Fri, 27 Jan 2023 12:30:13 +0000 Subject: [PATCH 13/22] hack for replacing subject with trial functions --- gusto/labels.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/gusto/labels.py b/gusto/labels.py index 535dd600f..9bb493d76 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -186,13 +186,20 @@ def repl(t): if t.has_label(perp): perp_function = t.get(perp) mixed_new = hasattr(new_subj, "ufl_element") and type(new_subj.ufl_element()) is MixedElement + indexable_new = type(new_subj) is tuple or mixed_new if idx == 0: - new_form = ufl.replace(new_form, - {new_subj: perp_function(new_subj)}) - elif mixed_new: + if indexable_new: + new_form = ufl.replace( + new_form, {new_subj[0]: + perp_function(new_subj[0])}) + else: + new_form = ufl.replace(new_form, + {new_subj: perp_function(new_subj)}) + elif indexable_new: + split_new = new_subj if type(new_subj) is tuple else split(new_subj) new_form = ufl.replace( - new_form,{split(new_subj)[0]: - perp_function(split(new_subj)[0])}) + new_form, {split_new[0]: + perp_function(split_new[0])}) return Term(new_form, t.labels) From f5f653b4e98c9e06e86d275a61322b3e40a38166 Mon Sep 17 00:00:00 2001 From: jshipton Date: Thu, 2 Feb 2023 15:24:24 +0000 Subject: [PATCH 14/22] work towards checking for the perped subject in a different way --- gusto/equations.py | 18 +++++++++--------- gusto/labels.py | 38 ++++++++++++++++++++++---------------- 2 files changed, 31 insertions(+), 25 deletions(-) diff --git a/gusto/equations.py b/gusto/equations.py index e8d51a9aa..b61872c5a 100644 --- a/gusto/equations.py +++ b/gusto/equations.py @@ -665,13 +665,13 @@ def __init__(self, domain, parameters, fexpr=None, bexpr=None, f = self.prescribed_fields("coriolis", V).interpolate(fexpr) coriolis_form = perp( coriolis( - subject(prognostic(f*inner(u, w)*dx, "u"), self.X) + subject(prognostic(f*inner(domain.perp(u), w)*dx, "u"), self.X) ), domain.perp) # Add linearisation if self.linearisation_map(coriolis_form.terms[0]): linear_coriolis = perp( coriolis( - subject(prognostic(f*inner(u_trial, w)*dx, "u"), self.X) + subject(prognostic(f*inner(domain.perp(u_trial), w)*dx, "u"), self.X) ), domain.perp) coriolis_form = linearisation(coriolis_form, linear_coriolis) residual += coriolis_form @@ -749,13 +749,13 @@ def __init__(self, domain, parameters, fexpr=None, bexpr=None, self.linearise_equation_set() # D transport term is a special case -- add facet term - _, D = split(self.X) - _, phi = self.tests - D_adv = prognostic(linear_continuity_form(domain, phi, D, facet_term=True), "D") - self.residual = self.residual.label_map( - lambda t: t.has_label(transport) and t.get(prognostic) == "D", - map_if_true=lambda t: Term(D_adv.form, t.labels) - ) + # _, D = split(self.X) + # _, phi = self.tests + # D_adv = prognostic(linear_continuity_form(domain, phi, D, facet_term=True), "D") + # self.residual = self.residual.label_map( + # lambda t: t.has_label(transport) and t.get(prognostic) == "D", + # map_if_true=lambda t: Term(D_adv.form, t.labels) + # ) class CompressibleEulerEquations(PrognosticEquationSet): diff --git a/gusto/labels.py b/gusto/labels.py index 9bb493d76..583df1159 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -140,6 +140,18 @@ def repl(t): + f" replace_trial_function with {new_trial}" raise type(err)(error_message) from err + if t.has_label(perp): + perp_op = t.get(perp) + perp_old = perp_op(old_trial) + perp_new = perp_op(new_trial) + try: + new_form = ufl.replace(t.form, {perp_old: perp_new}) + + except Exception as err: + error_message = f"{type(err)} raised by ufl.replace when trying to" \ + + f" replace_subject with {new_trial}" + raise type(err)(error_message) from err + return Term(new_form, t.labels) return repl @@ -184,22 +196,16 @@ def repl(t): # subject is replaced because otherwise replace cannot find # the subject if t.has_label(perp): - perp_function = t.get(perp) - mixed_new = hasattr(new_subj, "ufl_element") and type(new_subj.ufl_element()) is MixedElement - indexable_new = type(new_subj) is tuple or mixed_new - if idx == 0: - if indexable_new: - new_form = ufl.replace( - new_form, {new_subj[0]: - perp_function(new_subj[0])}) - else: - new_form = ufl.replace(new_form, - {new_subj: perp_function(new_subj)}) - elif indexable_new: - split_new = new_subj if type(new_subj) is tuple else split(new_subj) - new_form = ufl.replace( - new_form, {split_new[0]: - perp_function(split_new[0])}) + perp_op = t.get(perp) + perp_old = perp_op(t.get(subject)) + perp_new = perp_op(new_subj) + try: + new_form = ufl.replace(t.form, {perp_old: perp_new}) + + except Exception as err: + error_message = f"{type(err)} raised by ufl.replace when trying to" \ + + f" replace_subject with {new_subj}" + raise type(err)(error_message) from err return Term(new_form, t.labels) From a5a63016f8c22ff90757194b7474eafb29181c49 Mon Sep 17 00:00:00 2001 From: jshipton Date: Mon, 6 Feb 2023 10:23:15 +0000 Subject: [PATCH 15/22] only label coriolis term with perp if we are not on sphere --- gusto/equations.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/gusto/equations.py b/gusto/equations.py index b61872c5a..40ada6ca6 100644 --- a/gusto/equations.py +++ b/gusto/equations.py @@ -663,10 +663,10 @@ def __init__(self, domain, parameters, fexpr=None, bexpr=None, if fexpr is not None: V = FunctionSpace(domain.mesh, "CG", 1) f = self.prescribed_fields("coriolis", V).interpolate(fexpr) - coriolis_form = perp( - coriolis( - subject(prognostic(f*inner(domain.perp(u), w)*dx, "u"), self.X) - ), domain.perp) + coriolis_form = coriolis(subject( + prognostic(f*inner(domain.perp(u), w)*dx, "u"), self.X)) + if not domain.on_sphere: + coriolis_form = perp(coriolis_form, domain.perp) # Add linearisation if self.linearisation_map(coriolis_form.terms[0]): linear_coriolis = perp( From cc1e7cb9b79379424b0a5882da61d70e3be053ec Mon Sep 17 00:00:00 2001 From: jshipton Date: Fri, 24 Feb 2023 15:51:21 +0000 Subject: [PATCH 16/22] adding simple test for replacing a perped subject on the plane --- unit-tests/fml_tests/test_replace_perp.py | 48 +++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 unit-tests/fml_tests/test_replace_perp.py diff --git a/unit-tests/fml_tests/test_replace_perp.py b/unit-tests/fml_tests/test_replace_perp.py new file mode 100644 index 000000000..e3a4da2a2 --- /dev/null +++ b/unit-tests/fml_tests/test_replace_perp.py @@ -0,0 +1,48 @@ +from gusto import * +from firedrake import (UnitSquareMesh, + FiniteElement, SpatialCoordinate, + as_vector, FunctionSpace, TestFunctions, + TrialFunctions, solve, + inner, dx, errornorm) + + +def test_replace_perp(): + + # The test checks that if the perp operator is applied to the + # subject of a labelled form, the perp of the subject is found and + # replaced by the replace_subject function. On the plane this + # relies on a fix in the replace_subject function that we hope can + # be removed in future... + + # set up mesh and function spaces - the subject is defined on a + # mixed function space because the bug didn't occur otherwise + Nx = 5 + mesh = UnitSquareMesh(Nx, Nx) + domain = Domain(mesh, 0.1, "BDM", 1) + spaces = [space for space in domain.compatible_spaces] + W = MixedFunctionSpace(spaces) + + # set up labelled form with subject u + w, p = TestFunctions(W) + U0 = Function(W) + u0, _ = split(U0) + form = perp(subject(inner(domain.perp(u0), w)*dx, U0), domain.perp) + + # make a function to replace the subject with and give it some values + U1 = Function(W) + u1, _ = U1.split() + x, y = SpatialCoordinate(mesh) + u1.interpolate(as_vector([1, 2])) + + u, D = TrialFunctions(W) + a = inner(u, w)*dx + D*p*dx + L = form.label_map(all_terms, replace_subject(U1, 0)) + U2 = Function(W) + solve(a == L.form, U2) + + u2, _ = U2.split() + U3 = Function(W) + u3, _ = U3.split() + u3.interpolate(as_vector([-2, 1])) + + assert errornorm(u2, u3) < 1e-14 From 90e813637a5eddbba90b81021c536b904b6c98fc Mon Sep 17 00:00:00 2001 From: jshipton Date: Sat, 25 Feb 2023 23:50:29 +0100 Subject: [PATCH 17/22] start of SW fplane test - needs checkpointed data for KGO --- integration-tests/equations/test_sw_fplane.py | 125 ++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 integration-tests/equations/test_sw_fplane.py diff --git a/integration-tests/equations/test_sw_fplane.py b/integration-tests/equations/test_sw_fplane.py new file mode 100644 index 000000000..2b3ccb38a --- /dev/null +++ b/integration-tests/equations/test_sw_fplane.py @@ -0,0 +1,125 @@ +""" +This runs a shallow water simulation on the fplane with 3 waves that interact. +""" + +from gusto import * +from firedrake import PeriodicSquareMesh, SpatialCoordinate, Function, cos, pi +from os.path import join, abspath, dirname +import pytest + + +def run_sw_fplane(tmpdir): + # Domain + Nx = 32 + Ny = Nx + Lx = 10 + mesh = PeriodicSquareMesh(Nx, Ny, Lx, quadrilateral=True) + dt = 0.01 + domain = Domain(mesh, dt, 'RTCF', 1) + + # Equation + H = 2 + g = 50 + parameters = ShallowWaterParameters(H=H, g=g) + f0 = 10 + fexpr = Constant(f0) + eqns = ShallowWaterEquations(domain, parameters, fexpr=fexpr) + + # I/O + output = OutputParameters(dirname=str(tmpdir)+"/sw_fplane", + dumpfreq=1, + log_level='INFO') + + io = IO(domain, output, diagnostic_fields=[CourantNumber()]) + + # Transport schemes + transported_fields = [] + transported_fields.append((ImplicitMidpoint(domain, "u"))) + transported_fields.append((SSPRK3(domain, "D"))) + + # Time stepper + stepper = SemiImplicitQuasiNewton(eqns, io, transported_fields) + + # ------------------------------------------------------------------------ # + # Initial conditions + # ------------------------------------------------------------------------ # + + u0 = stepper.fields("u") + D0 = stepper.fields("D") + x, y = SpatialCoordinate(mesh) + Ly = Lx + N0 = 0.1 + gamma = sqrt(g*H) + ############################### + # Fast wave: + k1 = 5*(2*pi/Lx) + + K1sq = k1**2 + psi1 = sqrt(f0**2 + g*H*K1sq) + xi1 = sqrt(2*K1sq)*psi1 + + c1 = cos(k1*x) + s1 = sin(k1*x) + ################################ + # Slow wave: + k2 = -k1 + l2 = k1 + + K2sq = k2**2 + l2**2 + psi2 = sqrt(f0**2 + g*H*K2sq) + + c2 = cos(k2*x + l2*y) + s2 = sin(k2*x + l2*y) + ################################ + # Construct the initial condition: + A1 = N0/xi1 + u1 = A1*(k1*psi1*c1) + v1 = A1*(f0*k1*s1) + phi1 = A1*(K1sq*gamma*c1) + + A2 = N0/psi2 + u2 = A2*(l2*gamma*s2) + v2 = A2*(-k2*gamma*s2) + phi2 = A2*(f0*c2) + + u_expr = as_vector([u1+u2, v1+v2]) + D_expr = H + sqrt(H/g)*(phi1+phi2) + + u0.project(u_expr) + D0.interpolate(D_expr) + + Dbar = Function(D0.function_space()).assign(H) + stepper.set_reference_profiles([('D', Dbar)]) + + # ------------------------------------------------------------------------ # + # Run + # ------------------------------------------------------------------------ # + + stepper.run(t=0, tmax=20*dt) + + # State for checking checkpoints + checkpoint_name = 'sw_fplane_chkpt' + new_path = join(abspath(dirname(__file__)), '..', f'data/{checkpoint_name}') + check_eqns = ShallowWaterEquations(domain, parameters, fexpr=fexpr) + check_output = OutputParameters(dirname=tmpdir+"/sw_fplane", + checkpoint_pickup_filename=new_path) + check_io = IO(domain, check_output) + check_stepper = SemiImplicitQuasiNewton(check_eqn, check_io, []) + check_stepper.run(t=0, tmax=0, pickup=True) + + return stepper, check_stepper + + +def test_sw_fplane(tmpdir): + + stepper, check_stepper = run_sw_fplane(tmpdir) + + for variable in ['u', 'D']: + new_variable = stepper.fields(variable) + check_variable = check_stepper.fields(variable) + error = norm(new_variable - check_variable) / norm(check_variable) + + # Slack values chosen to be robust to different platforms + assert error < 1e-10, f'Values for {variable} in ' + \ + 'shallow water fplane test do not match KGO values' + From 219e31dbe55b1e1d3094892223de2c938eddd3e1 Mon Sep 17 00:00:00 2001 From: jshipton Date: Wed, 12 Apr 2023 13:34:55 +0100 Subject: [PATCH 18/22] add test for shallow water fplane checking against known good solution provided in data directory --- integration-tests/data/sw_fplane_chkpt.h5 | Bin 0 -> 151160 bytes integration-tests/equations/test_sw_fplane.py | 21 ++++++++++-------- 2 files changed, 12 insertions(+), 9 deletions(-) create mode 100644 integration-tests/data/sw_fplane_chkpt.h5 diff --git a/integration-tests/data/sw_fplane_chkpt.h5 b/integration-tests/data/sw_fplane_chkpt.h5 new file mode 100644 index 0000000000000000000000000000000000000000..d4107f56a0a9c671f4088a039595e7c9fa8f8b86 GIT binary patch literal 151160 zcmeFZXHXT}8ZAnYpkxJ66piu{>-aWhOoJn`jp8b7ej4#YyYZ;=Zq{6U;V+$GS zKU!L{O=Qgf`SaiZBmVu*TbtrPr~WlTqV;WMn)3o@f5M z{$9#|PU!qSudk=9Ooo~NuK(}&_cGZ5>VN*C{rALw$N!@cP*XnsAE_As_4;3x*8SDl zzhKxq*}2+YASvme|J44y@c(lBSFdGS|2!z~?>h8k2mbRHIT;NZJsB<8IS)^F8$6d8 z0{us$3V*NvbHl&(ki)Y7dpxlGjfa@O@sLOL&;0*gJmiJ{Gf(KB;lKa?e;5M)Ydrj; zMl!OVzxw$X4yC_;C|R97fTNnRLROEjvBsuwCB>zote@9AOCtYFLOZXr8v>B zsc@Ebor0X=KMTa%Gl zk)U0rJ9DLhNLa5OmeK!2B>W-|DHjb>iBrr0mmD1-iW@bH~0qu)%XP~PU|H?pwV&wy}~a7(D%^a+;Y;VfX&($MiEqLWz!Me*y5j`ZUm!61d+cKjggtzdp|k z%+le$q{c_&30NnQi9eJE{%@SwyL~&JD^wRLRfy-#B-P5_qa;(&g~4z+Jomf?DF%$9+>xSl#tBCt|RB>3Lm#Fq#=#UJfrWWeZs zGxH&R0g$aOr;ovP5}3WtO4N& z*uCky-A~_bWC{+Qx~l@PCnarFMS+P-!7_>L<8Pc7t37%l70+GEd`dP5JLW~#3-fS% zv01M0ZOsMuIm&0G)+2;xzRQ zW7g$mBH>=zxNh7ZBH@`=7{@R0lGkXV(ExjQX2?w@*AWRpmU~jBT8RW^)``|D@Q1ol zBFP%pXFb>KNxaFc>bl(rk71IDGP|G;SH=+5-WDQ(KTGD7J9rb?R{V|lVTXRa1j#(=k*g_{p)w9=B~8mKTQKALRM-*O2F;MU4~M z^N}A9UcJ=_|EqU+zx)9{rjlLl`w)N6)4Dq@<9_Y_>Q)oPPpn(_d)8fK3ZW;<8=u0T zsENak)zHD8M}9fe9ApZeS8Jc^!=G=4oQVT?{!|B#t1awWzA*T}9d_!a)>3gnFYCgo zZ-#(VQs%=u#n3_LwR57B@ZUv)yg(FwN1tY}><13FRjk8$vF_ap=PqtMZ))|!`Y!ml zbyCK#4EFhQE-TUQ1a9xcBNREw6zqM)PYi>H1)UB$kqPiIvv+H}1^gsS=#$pKgX6{( z9VOskf5wN_2L4NoalGW@CR5;(fBbbG=WiZ4*c%O8H(&GR>;$fi2dhuNH5=H?mT?)2jXZ4P*2 zV2C?L%2%0256(=(?!N84Vs8+KRYzEk)kF6pSw561zmW$QLfDMK152x~SOzWlh~aTQ zfqI}l-|(q*fk+^aIBBs1b&WFj-Ep5`;K8Kop^Ciu=lD3?F2t+BiM%r9Qq+OzaIR#; zX`EDi;dbB>$9*?M2fE;Ge$(ZQ^}dO^-1v$6&p5Fb`S6e_+$*{Jd<=emJ>j~&3wpGh z>u89BUvgsIF=t^neNIn+p-B-B(@+qDd z{3V^@H27-j=sr1u@#&qX9^t^<>Pe3GBiI*Dr@(##d_8x%n#2#gFY$@6+QBX_Be5t0 z@aBReO7Lrd@uts}rinubhK0J@sH zS?61a{K!k;ef$dg5@&Jh=jO1JOD}TM6!89egvlTp*C!YvPcY*;i&%IJ7viG+!JCoN zS=1R~ptd*suRWO2|Gf`+QOex75B0$L7IUHs;=n_RXX2oZ`o*Fw5^Kh`i#zb0ApSI)Tt{))`{;d0y^ptR2r%c#x12|ot?mB13hx|pbI&uWZBMO8| zRKU@%o1Y_+fJt`xr`vMck8A- zt4BW2?OdQkzsPaG*iQoZgnyb5{0=<*mHLc+Abx$9yqj21?~fS;C((5g2~s+0c^`o5 zak=0c62A4nbLA_(5D9AiDWXfbuBOhWcoq3WDOz|j1^L)ZcR64Kc_vqelW`F<(C@nu`bi?iiwZCjfsl%D5PNp+{xwenUg}q4CzcWdyo< zp;wn?iFmI&uKk`J^YUb8cgx|vq3MVoGUUZR*Xdj8$YXh|+7TlZ!0oho+{z|01@eby znohvpsq4m*J@EHFmr_kLbkD4@rXCDjR#iJve!`Et7xGJ5ppVtI2b}LPCTB5B+hD$R zleml&>@PH4qOybiDrEeWrKpb)F_AmZqR;J~3{`Ug-yTie3Y3VSV)vwzj>tD+4b{)9 zk)P9kCi0SaXP`RR%nW;7<}A}oVSmeU)yPBOUF7T+ds1C}G@ZMU1AlkPsG5I4e%oRF zbDtU3)fc!EaCQZK@r_%5+`vnP^2N$p=x!*(dB?&EPXNl2ik~)ZGQWnilHu)=}6z0#&sRa{HyVZuTU!EjXL;~AGQDC4Dxl@ zPa44<#GCR%V^dNeKj4zC=Z`q4vFMi})t|Uz)WY zYN*He5^&#KH7f6p>l8tSifO15cd6*a-H>N$vPI=BFwarHoS_MIRaJQFlnUa*PUXX< zAlw&h5jndHqiTR$q$utW4~4S=*u?xPcVPM z{!M7)aqCZ#PG$EUxsLs)>rLuY^T-2tF0uzEATJ!-ucGCP_?wD#sPo44o6gFWw!lTs zutm-S`P{77-)RkYTaWKepoJefkE`paFfXg@jz<}uGvKcHD;6UsIkE3mIrb}b_vLRQ zPc9#*u~$SLVlf_*@WVRnlCFsk=!dL?g{FDow}@G3%RTUZRrh*yJ9K*R&u3n@Kj^2- zx%VH3?(TJIjh8@Q)JHBRX(A7Om~sA?34Np`xI7z!j%%ZbuDihBk@Dc!8(5EHzS3s_ z>+$56T+Bj!DwP{>&Vjuv&g{)Lhzs4_(s3=YYhA@Igp`M6#F&le(D&c$l7B_o=LM-$ z+PR>wSfbGjwkjhM_Q(=kuOM%`s#Hs!f&F?DD-C+U#dKj#G9UhzEqSLefOqXeLe4(; zB~X^zFNJy4GxYp&E2x9a{u4~l(;3;@71*-Iwzce-9E&&0h}KOjQ&iAJ{3;SNQZ&f z38v`~f8aS8u|MSz&hO)@`@RFbD19AuXofC_DFtt;BY%w#?K`dX zasP|oM;v%@|29jiN2hWBhxFEMU%=Cg@USzz;ME~fcv=nm{<2kTWCrUt=Y{VphpyU< zx%|ozm)!<~MdQF*K|4x?A9(wm>gl>f4SFz)tQ|KKW|Iq02UgZ@kFgG4yA68u0vu^*NeQ)CtPw-eA&ukF~e$ z=WZ{-Zw8|RO4u<(#cO2_+@+)5bg+ZBPp9788-+eXZ<_a!cpqcPlBh@g%ZI;dX@~CF z1m?bIARcn}u0HJrAA1bScHdk_oE+Wf`4+lvea~8a8}_lEK2;HodZs*V*O3q0=qy8; zB`}`;$ld1xz53K#F|0%$2#TAo7=gWB(vPPNVDHNu=M%f2-zsrF-wNalqb9NYx3KQ@ z;}6FK;fMBGW3C$FpLWw~5hHl-S$)otiG0X&@m515{K~X2zt9LBhtx7{QA7PTuT+aj zZ9?9e;U&BUPyW<RPgM4!Y?L(5r&bw2OB%dm!h*9#i4tZ2jm_YMK4G2nGNAp6+kRU$$1 zd(-B4*pu)!{!$xsZ{(5lf)_YbZFMUsMf`Y@>%J2OP7y|<@^@hG(V4~KpQt~#JF{D9 zael1(^>ZImolFQ*-vmAQ*GRQ&pwDtSBTUbQKEN_^|0ZeZ#`EKM<3q@Eh1)F44;By! zggndM&T8Q7<-NQI_T^7REu75-KB}jOZsWS=k?VA4VE^lstpZPgi`o?;MHceq!*{n0 zKH*qX?+5Eo`2XVOC!15KbF2nY=j?&^!7#EJCp>qTe(dIc*n2YI_eCxEyOljiM*uqJ zdn#x{3!d`kX~RDZ15bLMdrXtq53TrnszO(DQpPj_@TcI7>C92+s6yAP=MM02z7)yZ z4?MWq&WA`tFG>|wl?S0~YGPJ$J?6(rtn@JB{A0Gpya!n4go0zAHt?Yq)Vgbp~pPneVx*58(sLe_Hh`eGc`>TCtFOkrdKfIL& z`7EZ6=Ilw>bAT~XpB8!T?ylS@`f1=k(K0E3`TL7$XK%xf61A22S@0xsj>$LhKYkB* zMF5BBq8;BDfq!6x=50pcCMcoYdmeZnxI*=|96X-6B7FWV;*CbED*6!S9Z1kRb`Q7@ zvJ|-;!SSj-Yq8c{ioOf)?v;fBin5 z20Z(?myoMtecG&n%0`R~3Teb`ScmES_i+uZ>-qZ7SzYk5vxq$5BA%C~R~R6O`m#`S zo9Ki4O=ki*CSd>7E!`vBm_OMmpy33*eqWQ_lm>g5&o7_4iu2_{udgUT2a#X*Od`<| zatsr5R#9K2m8GOe@5kfXuT4jQ_xN`19v#^G%HW&{5Afzy)I2T@{g=9Sy=g$)XuHam zKSBN_&r0s6!x(wOzvw#hE#;}jrVGdmWq%4SG+-A~)Usk9>e}HP|7ru^93vQJu?SqL zbq!8XqrUwTvhVgLjbU;%1lYx~%-_Zizcg>Wa(@dx(~D-SiWl(S_62cj61ppYp6pQx z-WF@aw@x9BT|O68Ax(i6WxkLdzpPe0NIt9D8 zx3ouUBhIszon6X@!0(u7S{Cq{JE!-8^gi;4kwU@^jPi-S0;%Xrn)Hktx=y`>CKea7_GU zJ?00z#BSwSklvG7NoqgYz;($ZG~zlKFQpP~=doVk*Tk={5k~^gieHfSSF{Q@$w>RZ zSn*hHo;lRb`gVOw+`eZ?pGO?E3re)a^>m#ex$WqSsy}LH zETE1?#VbZ}Bi}wFPiHVf9XF%MuOz(};~H?AWWe)`8Hi2W(Fg3jXYfo1@i71R^D`;n zKxMaAD)0OM3+UVfOKAth318 zctaoe9~s#8OAU5)s4{oRtzqB3>V4)3?6|?sIA{PnGh%dKlXO<7o~Bj@dmn{cu8V{B zY`y!JR#4~aGc;#BVW0Z&P*)bltZ(G^yO94}zl=|{g5S&Gl5-(=j;vmBJ16w?b7=Z) zbt{oD;-VuS-HX2Fz5nbe@>j53^NpTG=zB8rvmX4A zjZsXTuDGnkM#Pujme%=+)Oi-kLM+u`5Iirdw`MFZLLzPc%PB@n27*fBlDX$ zt$>%TNXXAEcuwWjT}I8&-KN4=>VxPni7r70ccV|Hwp^>Vga5RfQzuF9YtufpM>nIN z7u+MRtAjdW6Yn$l26^SO-yprfD3OpMFPB$<`dVPs)3+IM@zjP|IS_f}aMIjE63%@} zIlGMDx7~?Qx&w&2>8-ohR*+xv#?2nd!0rRj2N|)%9 zd3)O#do`f%S!Sz}JDBe(z^O~2D{beay1U){XvG(WoCUwmP-bGHui%Wu<_A0x#n2-Hyl|A=$mgN7InSx*wkZMM9dEYYH-+wluGMyu@b+<$ zrS(BR@R%8sQ$id&$2m6h0eA80I4fi1x9hKow0wxiqTF|cROl?6IW@W!eS)L9)0sV} z7oT_728&}K=TyVaCWpR{&+;(@a9_y9aa@Zn-nB7V}Hr z^UQrJL!J#(c60^Ls$rGeR?vr*^Ng<=!yk+3l@pbie@8lt};4mn-8*!;U)O@oV&+XyeDmVq5ET3p# z8A3b@?_pR~4ud7U#j#WG1kv`k0H3ej~jU^ zazWIEZ{R!fQ_%jt^dmzh&n3|}GYAu-lF|Pg_%5EZhP|uBpRYW`^V$zzV?6^s{&})yn1}r1XqZ`T z4!hs*Sp5ir-6FY1dam(-Z#ku=NsRX7>%LZeNR2IpJBf2ulX`^=&JeM z$_5|$8q;L&1S#~r>sN(7kUkfv$8#Qf3%O0BZ;-4mc|k-Tp^ja#{0N=s>&d(n zLw@1)zi`kE{_kyTYMTum3vl@O znrq@ba5;QSV<^@AOkeYhvjA`Vpd?FqdNy-XBJ%19$mi@(&& zRD=#r?7Dq-XbJn=j;Xl2uy;_yl${EA6xGh(RR^yP4;u115g#&DJGq`q8?&@Ba134?c!9;m?8pd5iOmzn~|lKaVcaAn!I^RcDaH`7Zm^ zy%%x+TNS;bXSjc=wWY-wI6u+xS!M^XU#ts1c7YEIfg8jq)Zv<{?5)S3JDcD~rB@3hDbaWz^1%+YtxuBHBid&`ZtHmmL()`#?A+*E;z6Fy*A}zk)o`Ho5!>ad@)0 zI{7R73*5gvx&RyvU7v*JqMumYZ)cPP9{B`%RKgIy)9xnQ7IEJC%<1nxq1%zkrvqQW z!z0Fk-?zcbYzUK4FAw^+ldr>*q1Ox+37Ip%pJ$VYn-=i@D3ik2gR$RAZS^F0Esu69 zF@@eAyfif`h8_1u3zo+4yweWGn#r)IF@!qV2l>E9v?<90aoqI$^etP&pR1hx=RoMm z=GWEelrHdKb~f)l@cbwh?aB%LKK5JWx(I%lN`5ywARlQ(N$pQXUhOs|iqwFo?MAJB ztKd;2syoUH_g@UFlI{ln%UKb;UC4*NQd#aUn5T7d%c0lEqkRr!L>BaWh8xVsdx7g! z$-`^=$njo3H~HXM;CgOXQ_6MVPZ+&5I|@ALB3d7Z0f*8@F8d8&fAaHelXBQ!v$S4) z8G33~Y5UHFTclK)ojka132Z~Iwf%*eW4ZQwPS|Ri)>r?g>SGw>s3{o6!d66 zki+bY=VY-_Y*vGPGfVTk<>Bw`#Y|lj^he(D2iXi@zkm!eNd?z)8EqOe&|eC?6%2ih z;{nUL^Ni4g%uBi{W8lL{sCH08-fv9H@v6doZUgUHOkj`MkDc25$X{K*MZ+1u`@!E^ zOw^z&@-C)|T<~`zkG9(i_R4qfkT!##`Yj?IQ}BECvGI!@;GB>(Yn6+Du_3=YTGEDnAVXeDw+Uk)$NQ%rHEk~uzf;>x*m+S$zVF_5sttP4 z`p)X0hWzCn?soDc{4_j0^EeUqBwdfSssqlB)ibZuVYhEi&`l=Ttw{aO=?~&D*?+C{ zHtLV?PUYV}#)6^_t1y-H2<`jWAaiN9zo$ZzEIhT59IfqLuDjl0+{)mGJ}C_o3A>yG34!0+b4 z)#s%5L@rt_Di(W$K%9gd zc#@mQc_d3*HAg<}0au#r0;h zgH?MlKWU|fuMc&k#hSU;6ZQPF>mlV5jAWAL{O^%Ji~_!1BGr*iqXUGa=ws>&`Q8tK z&+v`1syN`0Yo6|?g8YBd-TI|0aQ?wWmgxuHsuSiVBEiS#Qp?^?;NyMK?xAAfzOR2Y z{tCu#6;9=wab8!b_Gmih9jTQ&ID)tqc7Ac@B%ZUN6*^agy3x*WJ24Ag{8{ff+y{G% z-#*caz<7ErYA6-`2nDY;rwROMpjhMi-Hm>xK&+Y%IJ%}>laPac<7xLNn8IHPHG9rE z;88Ou;C2V&(Y7;dxnvJ1=q-p!;rBY$7yWfC`WNCbc2ez0KJ*>D zc<=$ML*owozHm~Oowc?5Mf_tcS>4ZuC2ql%^)c$|@Ww(SCNP2I$hyqgdE+`JH- z8>l;{>c*QGVOMK$5K$XCo_ib-V1fE9SE0hw2fsg2^S^%xz3%XH>5YS3<)M#c&EfZz zh~j%W(6xjtp@8&#!6TofpOe0KAiCp7tp@z~rI~g8De^r3tqMH}{Z+y2m2_8{nCP%V3_HtcsP+_RhXo@St=+0YkpBx15@FVJe$g3T3QxQfz&jwx7+w;Y6fnV9U z=Lc?rU)@=eHeno#Mr0^oL)|bJQ?EV-o&xoEFFpcqKVrR(Jx1Q?QT^kv0{$KK@61QT zKAlM!T@mPSa*+IHC3tY^e@h*JysY?KM_w9!NLI~$dJUc(^Plea!kEwaz;7@3rhD4L z-h}w+uX*y|DB?$R%dX7xSZ5$Q`pqxkaRhxz9`xnomLy7rdYMPqt)~SZJdWyFhJhCo z;cpAI;7xjW<{lx~t$M`wVln)@r4V1g8T{BwyM?;J-XVup|0eJsdy`*v41V2KkD*e) zdR9tazdmE1Yt_V5dks2Lx^}*z7I7)4)qVMv-{s}=k7BdoW7UZJN;D7e|izi=jJUi7DNeB0`!{ZK3H!2$SLsbo+S41Stlihj`rem6VpnzkT+ zHd{uEOTzChT9#Y?;5?nr!Rv82-Zg8QOsZ2wyw-tiI6v=@{WJr4bbQ{xT>$uH{Ol-d zM&A>0?Qmfb)-mgizO4n^MOMG5llqE>YpqX6`}R|k4x@^&lb=(2pB8j8+ZjDf2VP4H zjIZ&+-VgS!M!#W~f$*$xI`T;3$jRCd(8=Yey@S@kCoKKt>u=D-ejbji&ZuL@XvGd| zV4qU&@-``C-3$0%tL@zgmBz-i-d6 z``SVi;?mIm&DjIEzx#N-MiO{3A!i{BgVz9gq065UC;aNy%}D#6_M1ni*O3qD^0F$- zePUb63 z_#slTap^1a*3&Q-)qLb(nxB^qqM_5|_J~Jq;78t}sn`Sfu>Mw7vcUY~XY;(cDpBWm zOvzpZ?&)q6Yn$@1|CFA)&RmMVZsXL8Zp25YMqHIJaBqCx6qlaBTK_!rmM zasJ~uCu%|D^JK#)H!9RA@+Ny9H{`KzkB^ZD;km6Qg0^Hh?s?9_q5|I91h~_qP!|qV zk2O93?h@a+SY#1Lv-^en$>As8{jjMs@cU=lFmVp|oe^te2}1pSE~_cXfakVvn)eFA z^R$j8efCEl`Pt8XstWnjf9@HB8|DcbefhKs-5STsyUPQw9JM;xyNHAMw(wwO>=z!t zYWtK6eq%OxElA%xYjWQ5@(}V)xLILAA$YkzYCkKDe8jvEb1)q`ZL?RjCVgI6vBxM> z4|$hK%~RYTe6yy@eH5!ezh9Qfx2pvD*)&|HhWwzPZFYpa1K-07Oj?#0M82FV+cgM# z1G)G&o*^$zQtEg6HQA=C@FUUJTK+Yjx9#w8 zI&qBgY&~K4k-35o-e84dzNo2=qSju-FNOG+FzgpQmXz9zQI1Z-ZyR(+&zY6F4m>U8 zF1g)B+(w+crRfLV{MeEvLk+)mSgx&Ud<5?X>!v!0w=q4l8(rWpCUU)R@dxVD_I$Yy zI3K6RvZ7v1BoMOw>bNlf!LyZ*!|?l7gXpht;PZTuQmzqrD_1p}KS7=EvfcOV8R~b? z`}E#R@WZ{R=yg2ePT{BZU>4%=eIxfb4e0-n*S9x9(OPldC^U>G>Z>8A_L=i$9bQT<9ia1gM2Lv{{04%jJL4SzxZ zY}wjW37uY6xqFB7IoAiK&qOK2@wLmHMmpeqgKR-a74I9eYZ^b6?%#LoH>fL$TYxwPB2u*V;T@fd1m=p;vO~r={X( z)OiuNw^om>L|`3#Hgm-g;QcO0_5eNd?2Z^y)63xh%Jyi9DBxD=mh#2`cvp{QGRT2X z#Tc@0S+LtRbHGm-e)1-7cHIemy3D_OJ_&nsayMcxpw4`0cP+h*^&Sqb+u9*-aj8gs z%7(pLN1L4nP!}b42y?EjlHM!*;T%Uj3*Q~yH&ctg?(vO?#kw2}}8|W@E{zx4g=6U9GJK8|E)JF%{Nx%QNQ+pzm2L0BgTHHhg`W5+vwbD1h z)#>+-{9))}Z_<621K`Oh=a?=Bczgc+*2*gI=XKh{a|L709VeGAj6M%j-k(C8G*NVK z{)zF-JaaB7ex>G{6I!vJpf=6^S0UNKRx5hn)QmQRwtzadv%OJxQAZNEelP5z*M zteoYj!~W-Rh*%rLJob0oM^{wA>vRA5NEhT`!=>snPT+nj(2Q<3^p&;nQ;Zw9Wkjs6 zWTQTwAGST-1H6?M8Xphhxb*6`C1vzEu?Lt8J|TVvz9f>H0q>>&g5rOE=XCw8dOP^| z(Z+vHYPsgC^yA~g9+pzxcV{5!xz@L)8d!!V{3m$IYuHyLS+ecqW z-^rCihwjCQAl13h4zOmoA9_Qm?Hp`|$uk|VK8YywTAc-5+Ff&z2NYQDGJubAsR=<#;9l$b`*|+rwNM39Wg)M&|7qqg!ug@g z_a?-_hhACBZ|+`v-$}Lou^aXg`&laoUt{0VX6)>;4&F^Oj+1Gk9!r?)sdIonhSr8z z@Nab}{3s1kw#N82{m=(}_*a~t_U9OIzt?xOCIEF*KtRs(5%4ywYTC#{ygA-jy!qu9 z-seR0_Id#?*C_Jy8_?s$&av2O;MS8^A4>W?W0POMpx6}bRaV{n5OusJfXCq5%5Le!OA%ok&8{y!+&!pd5y!54B zX#{U$m2TR$m~U9@PwfJI-!!ZV-G@HLQt?X9Mbz1?4jpAls6TeH_X`YQucf%ufAUy- z#T)r9)EU`|MXgxyLWbA}%GeqP&@qSKOY~nQEf%2!5?Hnvc49=l=IC z&amITmCU;r^R@PST#bVtrJZXIq~G72IHj9PB-vH5mp2J^#jx?q$zs03X)OmD;4AxO zb!QOp=^QnSut9y^ICGHo9{jd56?SGqJk^vAP%$B%dXBZ(y+!;TY;2us1)n8bl7BG5 zKRdz25&JxRKa$NV%^$~NJ3W8L7vueo?$VD?=x{`kDYX*3I+>TvllCv>6;H$~p;vZA z;r9m6w|C&46?PJz#ZSm5!3&`}>RvMZICRb@R04UHYGi*$590U8L}f=8>kKn`c>%89PC9Kz2|KbYpG0|hsJ^?>jSPnnGjCgpUKIq7gcuKU$_Nak= zipi{lNbhH5*5kjCe)pSk`{Jkr{3{z=Im!b6x}Vu7oy9o&-JDUV2%lex(Kw%1nGObruNtFmtx%-OOyF0(8xvRIcjPGzuT_FCcTC3CGvX>Nbfm4KQYXE!Othwl9exkhqNl0 z|3&2Qe4}?Ao1mjRC%wxhkf%9ApPc%DdXk@bc=#>sNB`MPP zcts*}+P8zJ9u7m16!0b|AT1P)`Bcp>3vA$D%hFz&<0I(%`@4OwV!!mtFaO>I?D?i{ z68Hf4E$x-k{RW*Dj0zNU6Y;r%MO=>$a7&>p(G*7hdi2Yf*oS!dRX^_j0{$s}brY#V z-+AW9-Bi+hfes^=WhMB3*=K$D3;aFJf9hT~aN^*h;vjvm^wTcQcWTH3kK#`DAICh( zJ0;npz{g?i@(&q|wx*Y#?nWLlJKZQk`krozrTAR}p0i#htm%UJDItFfZ{s>+9Ctkj z&UZB$i>v_0hg(wKMZx}8mn0OfqrTm9$sJlnU*W8(uAYy4Gkko)e;)IkPg=$EAYM)t z%3cV8p9B~F?wj9$=Qw@+CGft$Sp7|~48I?g`#{!$cr<<*xt$VmXEf5bI|({gn!nVp z1RW@8pY7X&xIN&-vFRX2`fHzbTY$3^YnZzu@_+B*rMvXFKa?lF&=NcqE9Nf=AwE2F z9|@d8|3y}#RMn4q!FK(qVHNEBb;#iUX7Ie-PSraS`T3AgVLdtej8DFw`g+h;b+-ty!EZ4~Mt)PwFV8N&o(dk_RqMz#k?+`TzX;yM@mN4c%YERw)rI}gKJeFLd-BB; z#!bJjza;h3C5dlr#KDtj1-&>ae{Qz$8W}_U2_EcsU_f4~{#7gD3EhY*R&8``;C)k6 zQ_@lJNS;riBmLe*-gJC(9C!-7eecvmoM-9a={tnDnRKdjIk1M$yHYJryoG)@mh7G? z!|!_~GgnK%&-=~2Lt4OtbM$r(BjU19vr)4f@j=LIIF^gL6&OGf{{gtKiYBG`fCttM z9daSy{>HJ##|%7VH#e~C0+lBIi4SClgEx~n4o zh}NkAE#Q6socLTT?sq3npNqr&8p@-zC5XeT6JatP$eZaUy25kt%jIpu-D8Mf_3kal zNZHX*wO|Sv!`+CwP4$L;-vHwJP3+ZH=Qg9`T@{TvcbdM;Nz_wH`+zUPw9Q6<>p=ENMEihRu!!8CCU{>ff0 z>;8rFRR!Lfm*Ib~Qb^x7?2j(?@$Web{A{y(4&MKc&o5Z!8kiZeKEdU5`b4u(F+Dx&hukXF{`f3#;W#w*Bmj0hJz*6feScD9lrK{aBjxfA8Wz|`H|HP51;6)({gBau zU%odS%S8}Jlx-0-3-J4<^4Loi%nu*xp)W=s{`2W9QaL|K10e z%Gpa2u=nNpS#}e|vv{~;@Eho2y*ImI8gXsKV!bhj``lO5Sd5V`@7xQ~D}%n~)I&ua zkf#-my+u7>$0I>ovsC!8pklI{RKHA`wy%AGZnS2Dws`?>g8B{l^SIvOe!Iv3cJ&!; zl=i^Rhv%lQJwyCF%06b9_1>-4WeOp#j|(+#xFBDzn!d3N7{mLNn|q%W zf`7_K8RtZ>eoVs2^()Yc;<5o>F?bkoROw5Bj`ol7wmpWfE;&Eo8Nfd9QP)(s6znc| z@TGJNaWyx_Ho1g2W7?)8Ncugzyv2SRA;h`Y!zKk^oT0#h2|^ttViy~WVehg$By zQuraFzcjB5oCiLO@wWo^Cvno3lcAFU_RqGYeH>*Y_jy0~`7rKa=PLNLS4w8#M&Hn_ z62~gEj`ssUB=7EnpQq9U?PGDDHZM^$26awIL-%$A?Ca&xd-WE4%{obLe}MSnv%2tM z1-wLlX)Dx(e|uVIUNm5|rFJJaq8`c~iFz~$eBStFSdl)5&Jgt%_P~C2_jt$IaOC$5 zeuqW|Hit+6Q^S>)xe}T}JS;bt!Uh9p-nBa;VqfImbgSV!i=yrm{`5`_`~8 zKffh{3OuRh(+OyUUj~+&_ww;PQH!hMr=b68=WJd6ujmtNmCup#f{OFjp&>liP1o6+ zYcKlmN$>AHc)#`aGg0Loa13IgKT6uyN?o$v+yFo5?`Q3P3cnP+O7hb2ev;Yu{=f$G zQg2kg7J_`l!e?Pas>k+6^X;9m?)tFR=_%kXa;~J*1naP}@)StmdcO5e3kC4`Sj2iQ z3G0QpUHX%c{2kx3c#Z;f#xgd3M;nfJZZa5N!2Y;GNdGkH_h?p!pAr(_Pw^^0mmsb; zI?RQgA>ut`g}0{#bSly^s73mo`S}VOK`YcbW>!~~An<%jfIRi#XVh^mji(97_jNQC zn;Vgzzi7rJpTNBIq&uP8F~+jQ$5leNUF9x2n@~Rh30A$zkLd z&XL1mIXGYK?08H7%$M{=3v)gNEqcc*m+Af zUGfBYG<(^-UjTJ#h9b&<6Fi;w+_y4H%HK|fWk;vbN4vCmlHQBbPMQUfbU9T>r^N!l zn>5M>J!kPgbIISx8T`nTy>i$NJ-WuWF|DAF3}?LlY^WP?V23!}@ltkPSpDBBf-WeQ4 z9cJL#(gZuGU&R}nLI>=7?S##chemGhmb{1OGG&{!4#4h5-c<#p&u^oB`81=U%To{d zlf0mJTjB4{q~C4diSms$$37uJfNfq7{qEi3LDwz7*E4-8;|9i+BS$Vo^5Ff_jW4f* z!9!d4vkB7s)Prf8UoQgp`0UOT*TLh$JLfcA)EW17>QkG6`@;qUcMag#MlbL~6ME(; z@sVu=PRoOg_d3C&2AzZ!EjRj8s!?h|Jcmw{Y4;th%Wf?3pMMu7XD!;1^zXzZHe~A# zfloQnOOYSI_czmxR$J&K>iVhe?Z^xDX%5Z_$bU~MlW0lrWA1VDoN_`u?4$qUOgDtQ z_Rw@p9{E^`ZvWv!*mvE%^7j5a)YV%V3@x{@KYQI+4o{!JFI_8p`L2lvV<9W)KVuF+p?_JB- ziSs^MdlkMw9}^G!bv{E!e4Cheu>$WKHyo&#kq2$dCf$nRCv$H-+dkO6^`l&@JL;;| z&Rzx!==F=u1=dIC1HXM(YAC^ZuMd4;!pJ-FcQS}|h?@?-aRyrW=R3?T(gPf-^eDLF z5XauuI#uz&i*|piI4z#%7-~~`5qO^-%*;9fygier!<)d%_p=LAkD#9;D(UPR@IO;m zCv+Kpaxk1eHw!;UtzNf%fX=CnO$LjA!ww^r_h#T}yU2oFk`=#(3#Xet%Oh-XGi` zuH=JPcIOU3nHBguVOD#wA9y`dj%-9+m-(?|k={ozNo3eZ!(M_^ecXBIx}G4PiU9Ph3N-Whz$sjlFu*C?*w1N>}oRFAua z`KKa(9#n##e)rB)Ji&1tV`I8H@SvsDZMzSgA z{b7jDQGv@!(-?RE8fw!z1#nh8wH30eX zgWigSJI3dlU*DfYJQVF;q1MIy8C}j73vnz$n<91{@k747LY*7?*0+f?VKdO%<_3d3 zec-X_z)4~l){Ww$I8SS9H3+3FlG5C4o zn%csF-*_%R^XY8lu@6~41B`+1>pd^-1|VLT{*;M!g2z-3pVD$%FIUnOA$|W$@6h2H zQ{cTl^^pzfcLI~1*@e;IqmwmO+nE=7$c%k;9=xb4hd6%1I*LcPR9^%B3ktXPH)1^w zp(VvI@UU50#IqUtQTjEcHVj<7Xv=lkp(A1cA;*iz(;Yq>hud*JG4#w_JaiopO!sUJ zd8@lcyd)oa&qA;1k2CV?`hg3V3*eV+)I zm_%ye_Zf5gN;c%(cq@IucJwWVS9-}w-=|a9eJJx6#+Mm8oe!ZNW`500Oaks%BU^8g zzTZo4zLLVxh|kYbdQue7$9Nj4vx$Dh=dSN;HM7B=k=~W|YQzykI%|^$?Avu_^a3f* zu#IV%i$Mn`{PrjnNsTk&#_3^8IqBb> z$=Q1$;2d~r@gcgb0r#_;r6LuvZdOLDvKp>0zTJGd9%JQor~9dhmw=!T))AOTX{vfY z2>9z(3-Tx7Il>Np=SiQl3*@&HGJwajj~)5De`24+v}Yd;aOT@q;Uo*aX>d}>7gM0$ z6ykQ*f*xn*ikprNU>`MG(R~bk$q@#dQWoGbbm2jsIr3}vv5OZwfY){nyW587A532F z&4~rhoWY+{GLe7gWCf;MQ~ghUX^CQ+(oe!4`I`UfQE>4ZD+Cn}zjZ zZ-)1sXAS5d$hX_<)q%d^L+7m9F^}R8d3G0Yc=|2sHWBz3QJb~RK=)g|4|{WA{5(5W zaRzzI-?hcY9eJVnbU-Zu%^3?(sgw*2#zJMPWo}L=b3%s9GGvI%GJBcl zc`76nMY4z_vq)tuk}^e#h@5ACKJU->oa=ji&vpJdf1K<0J1>8<+Iz3P_Ve7&{oHHq zz1Lo^n`3&QC%fv)!bHK&gv_%uj*wTDOM^A%z<&3Pxj{-OwGTJD5Z}4m@3|#OJg4X0 z`okd@>IeSHHzag(YE5!`Jr!Gu;J%lVpZ;?f++X#EW$k4Ie$O>~@z6tle4sn-ItMt5 zEw>Ht0sCu$j?5ecUh-!r4Jm^k`YVO3VPGd~%?+h4(1B|dy~@P=-f%&q+;6~Fjg?4= zQxKQ&umwwh;BEEUITEp6Qo)m+SY9Y+g~ua_=SyOP44+U|+{nlEM_&C#M_u<`Lw#UsJ||kALC&?Z7Qn$>5(%mw(Aa@Q<~KxuX-x)KJ+g zJ>bU_qtsLn#NEv8RX#K1RSmCP!4T-V!(acvFVJ4R{%ei}xK6J`k2*0w@8x!BA7Wh_ z%lhyqPEg8r>=F3_e#z|$8#@U0PN*H-M+dmSJG-2YVZVKwci|HBccovVlk7qNoPTm_ zZ-RPkCAHVr;P{Xomlg57-?`3(8?S%|*F{c6vQR-@o;E*D0-h;{Xf6=zUKmV<{|$yBJPxF=md%7ZxG!k)N$cZ+cqZMwdEOApY|&0*I?$yuVG;KO za31ydnJW*Ww2eA{oA@50srDslAMm;L?bv21xGvwlc%?qT*_wm9&lcKEKC|!B0mvI9 zUZZyy#yvH&dZQ6gMqL`WO8*Y;t7$jdc|b27s(7sigC2yRPdM5P`IIqOreXuSKPE+0 z<_WyIk$3r8EbycB$i}zLz@xirr!0u~!il^#FIWMujpcO%Vjfz@p~tg#0p|+QQ>Na) zqtQ)miv^%(zuJ1QIzc?_*p5*zLdi^G946*tf33~wr~&-j>y=+9LwUflz(EDB^R)Dl z$#ob%y<54~girxYGyBmyrLwqo96JEdc@FkSn59)n1hIono*s^07 z&M!|M>23tROZ2#|RzkZuFI1L;Aa9-+-e-ve{)x?Odwv1(^rzC5bC~E~dhN`QL44Bd z&=RpOUgXYAYWa|t&OJ0ScY|fumRsc7dkwMg_0Db(X9b- zRB2vKJ^}ea@vHD{2DH0&z}BQ5@H}{z!p9Zj;j|H4ISg?+ctOnKAn+)+Ih|4k^uPRl zTh(=H=-;$|G9%#AWgp*iTfog;MNZfVa10JAu_5MVeq_?Iw1WJu=XT~Kc)fVy)}J=$ zU&m8l3$K5G?v9qbN9q~kqM}0x%}u=9{5wn z+Q2^!yi3k$c-jSi>xNb9Z-V&g{8~MC5A0=YmgC+6_6i?7RQw%EYCf0Ddf?k@2fmG6 zpu46NQiFA1x8u_L|E~@Cfxyh-n#cP>w+I!Is3j|gKe?%xVR7Exc|ksF*eXE%8QegR?yB#uIuw& za6D&dFy|KN=Q54!en+s|aPW>~BeeerYks^3JaKFupuItBK6`thiaL)R<}yCDzlK`=U7_0DPnT7CY+$cCT&WS$+t- zJaHsBE&_OOoVT?`67rxzDkzc}+LI`E5qJuEP;hV@5%Vy*Q}$X}gCAwToD1Fn{>NB5 zKBWM!s+`*!Uw}R`eq>q90GxYy4QB6xz1qKe?6e`jU&Mdi$_#$aUD~Sm2lywpf1|n# z_ch=QgZxq`DFANs>;t#-c0UpziNa-fHU)7A#kc9E++)(OO1+bH0 zK1MZ|&}XH_O5!`{f@S_q`{B8Qu1$@MS{vkwD z$S3BLLfjglYjamwDj$OV*OmpliRVW-jR8S#!H!O;A_qUv4-S_nv&4Jyeu^ZrFqAa2nc!VgoSABZcW$T&+7fKcb$1KQGs7_p(Hz1jaM8ch9Mo*2BDuC%=R{ z0EaiECNnzVW!ffDG2(kz0Zxjhuh0)RiR;VPLjDzLOr}_a4oV2khVcT=-`)H3p&WRr zDIc1;1obmM8W@`ZF5It}t%Jco<3H{G#C#Vshq1)#kgvm^Xp(6u0Jp5dK=Mdyi z8}D1)Y1r>}=(E>(xUSoN*KA7Q)o2g?LX2N34ZFX>Z^X$TcdM*qh5nK+^en?8h|}{L zhD(aT`$I1Kg=hhPBhLGly`U3yU9wLvLA{}pL;+%b{H<`+Aw;lA?l@;+t2qpzh)q7urr!0!{3px&+`wCRTQo_(*&=jdFEP)=W_Bp7~Ltf9y2PY(h zzTydr2x7fGifb>2Ord@yzmDZrz)^?V{mgghUk>*qM1P0;b$!hkqyqVH;lg#+&%g`* z^b_?vfHz(29~lW<_P1FZ)PgwJ=f&p}^DcFjzpqh4DY(L$+X%X>{7hk?3+l^0;$EBs z9trXHDzXEwM(?wqWP$kDid^ixo?syg;KN4K%4KJ`U*EZUvQ!uRxUYZwHT43#X9#z- zoPa#z&FwslARa+#nJuX>9@=rw{M}LDb4zA{m=EZ|bJ`7kTZm^(LIvMxz?pkPZlv%D zu@1TpNA6SDo}gXS0X;sUE7vywyp7hpDJ}(Zw8%-x3w;aY@7f)gnqdE>08%M8#DjLn zRpl+9%e-&LW6r}k=Q&UQj|I4{{qWSvUEq`Pg~|CbIR8)YuI*Z&+h1ObCKBs%$2Pc% z34lGr?(n%A;C}C5f~hO;^S7-+eInq!v;FhWx8SE}fLxKr3VgTmMKi|%aQU%iO`Z78 zNabw(b~n&X`kLqZhY7CqP^%`uoL*c_T}O zt2qnc-@^x!1qj2qq`z&GH|UYdY^Wx2A9!QE<`*U8v!n&4?uPqGv0^jZbkL^)}SjkfvtB}&NSE~f!B3gc+1=DIQBJH($U_h(Llec}RkcL<%m=dkzv1A0F2$Ldv{ zDI@UYkJgU7_i*0P0Z+y;u*>0aQw(vxRq3K6Jqq#vsM$Qv3**roDQOF09mB^x`3Z+1 zug(irG>XD_=d)T^W-ItT#n2SB40$oeq4l8}@|a=hi@`~V-?8I=nGS>hGnI5*RDkQ} zX^!F)@Z-Lj-1$#yF#i7bNrLd3>$Rq42-tn@sq{e~uvb7q<2X0i%ecM6+6C$>$2pyg zfw*~@%GMIkF=7MG2E;&q(Vy5fHwpd-dtC8vgFIZi|D=u;blpO~{K02vzdF=sv zXS+;;WwoFT622H%3VcG$>v5$JKY7VpdInH`v;WIu#P>x>B8~lzfsfk4+D4NwF5h&Z z+H@J>`hZF8nKqQCD^Xh&*mpdhsd)?JN%bA;-<`w2v$p1AEkn>hNAB3Y0X%Z&$tY6? z{|@XNYz>5x-Aa`wzZ&jyUnj?iLYz`$o>B#a-bUY07x)GKUA)4$@CUZlY1r(tA)l5R zEDDHuxNc|e@x}uG@UUt|AjI|8p968?Fz))gM7ePd^p-yPG^=?BJm>lq{plu*AIjQ= zwh{A+B^JV!Ep-Z8VH=BNK(ICF-Q1sqyQ~>4B;^FFVkgoyO3B`jD z|IDo)PZ95jSx)K2RYF|-MmneTV0%2E=+Sq;C&V`C-g(&HI*I?*LK)?w{3ZhSr)`T0 zx&ier1|u3`KD+piM8(sv-F5sv#e1+Lp*4s(3GlW!RI#099r}@4RoAVc4;rGZLZiT^ z(rrFZD}iSQHN4{+P=9X3>D49BV@L0rZ#EE*R<kB$T1Ml;I&WuV@EIu|->h@-v1CC?eqK?O+*SV$b=dC+m?H{i@` z-*fK-#7AW?-o+g1nQZZ7xCMMTMxW068|=L!R`7^sIB3g8fl>DL?bzx@YyCJ~e|lck$8xA?{mlKeZYm)-Qd?`)HV$556!NSV{D! z6lZ#S`N96i*wqc<_f&d1#%Z^Vzdv~d{-0y&V$rLN3mIs|%y`DKwL)@pHd|~4NdC-}~J<$Z?OYY?`KP?!aTs$H) zOT0(BE;3bp0`LlUbkSJ^KfksIpc|0CH$2!iiTNpfN4{D!gFS*)N*qIgAHRWzgB|RT zy#8|g2dHP^7)m7#I{9(!Ca)^&FWB2Z{TkXkx+Sjs3FP}>OL65w(A7ET3l@>UW1E$` zOx@5Ayej(3k_UO-%<|&Q3(#j$rYJ4|9x4ZN-XDa1sC~_`iw64TWcMfjf-vq$s+Adj z1@Y-_HTc8>{#|lRb2$p*+)4`ZVa84facbulUOp4PxGv z=V2+cmk`f0rk{m~b$(u$l}8ful8<~+D?q@9*vvRH;``%4t-{rKh-dlHaDGp~Z+{!9 zK@0TrN5P+-8qo3RykD0P#7ifvFs}pbJH@wBKM3~yJmabO7_M8n{c-IFh?m-_dzt2h zE=p^Pa072by;LO&puh6a{37-N@^W#S)$9%6y?KXv&^5r1`O2nKL2zAr@lpFtzytgA zkQi&glk4y{y{BO3)}tj9brAoV2&r6Y(4lL(M(V_Kj6HU$W*?xw6w}ioPuOnLyVB`z^Q}@m4-~G_eUX@33#QS000%eYQz(d`W$%_Pjs~2{UazM#X{pF?&*dsy` zs3z9^d7V|~9s@epI?Cg)4*97U*r(kBev8lz9!0?Ox87aqS#Vr>!+l*0aQWGD|5z2| zrK*g6Ll?w@>x7O5u?~5Xk#>g!_5CeWh|wJASJh~&z@vl- z;L&e|p|w)zZ*HynA1A)k61g0cKn?BI*Dl@Zf_`s@qMuS2oKHEacbk|8+Vrij#uMz> zGV!dDSie5;q>6?)_`y$;z!eRC8cOn|6W=@jIZ}Fo_}#NZgRjot0Nu3ps!(->{-m)% zvLzjGZ?&LLQ~xYFs!$Ob2XXx&zv}1;*ZE9;^GzgN_e#AecM13tKiTR^4|);N zQ`6rGagLN>^kIkVJEd?rr~)s)Qik=GLcjU;==C}dh<~dcmu?%x)x`Vkz%1bYPUIF#C=r~{Bm* zeZAb=31tH=UEmSF+Go2J_`t?z_b~;=9jl)iE{6f{_*E?LutL9>@R%k}6?huS7cB7`ct`TS zI7fURui9WpqFjghDFT%;d%)iKj|Oq&lW>2QN}V<|4|p+9t_cF(yLz6NRsr5|rTa=X zAij}zuSG0Cd?r_R&oe<>%;&yy218!cHT6I606$$G{BULE0i9Q9Rajqz_e`B{iPwiBJzpFHoaC&AMaN0}voX3V}htEllSU<5{ujctJXwTc~@*i5T zt7*GT#0pNeHxq|EJ5CA6;BpEgZ)Vtj;9mr<37F;G(8S> zL`F|npf1QOlVC5gW*B!cesMkl`e4Kx{`TD%thb-U8EOK$I@YV=;0XAhmn=L)tUJyX zXRW^i`7)U|N^7|W@68m-$0R}5LsRbC5bK$=(q(NB>jC~b==_0#7y28n`22fN+Ghrs z5cANe&Yz>8hUXZvG0)_J0Qc@ZTdyg|+eUL)b76R{^t_}ewFC4bUHEj%JBa_a^ZUO@ z!}a>MA9(T+>ZkWiiAF-))HEH1W1yW^Z{8HWhjtTkqf1|cy&wC}VL9OGiJvCQ))23% zfzHwjC~J0KP5cV^9mjj4{RQMn(Y^Q~;yvd_2J<^|Fz!;wu+h{7y!h_^k-ZLJc4&ySrO|1ewALh;|BZVv^zVA_4f+Y zf6%T04{b#Xv=0JrrFdGBE1|yN6XZhy{9@2;abo-n{rZzr*NFARy~b+N?!fg(94jh% zV0X38{hZ4Ka33JFN$Mbs4+ST6l}ACRcOGfd;RHRcSMG5W1D6|= z@IWh#SZ_zRUoYweoM&gwnaKj}Ec$jE?uGbo+k1^gAL2K=r%sL+wqxxMFV=$JW5L}5 zaZoS!XITaDJVn-d@}LTgOIS^$75#p~b6mxw!6?u%8y>ZDw19{2li0NTpzCVyiad(| z7vJ(H96_LCL6_1~=OIqZxos(@pWNt%O~l<&Ohp`LGF-m5*87r+|XW)l)8y5LDI)HwB0b_gB7CN1>j5 z{>?`tP`Wk+XJ|lPRfv~7{t9_E`)Gf9GhDC!x$J>H$ln;-S*L4&+rn$HiM@cw)$6of zo-qDNZ)WWk1O5B^`IciQ@Xe|0(g-o%YW)uV!VKsm9?f}m5d2M(oHYmp9NH{j)a-?J z4Bqb(bc6n9oNH8Y4*DO{n2nxt;DKgYwVVg=wD$0YS2>_(m*;;J{D$`Z`m;F0!M{TK z&T~9a!fMKC#5|W%-#ul101m1*Iz_#qzpQ@5wn6-^pvCqB{@cLMlAxh0ir{x^-|kUI z@QbSQTUaCT!JpG&xCeNaWm8+X2-j^>t83T_@y=o75Nm=sf9ATTs0Vs7r`Z_n2zqkA zZHwC~lxe0i=EV^A>6V4HJmAf}Ly5d8z@yfUJyA!Y|915A*S`<#ahhDo`ULqbI-yrI z1UM-cT>df%dEw%-Z<%=SD^52XR}JxIEUjN6p8t6JKfK}r=kK~zMk2n4AGKTIC4To# zhiaSDGL#o3EsqfM>=T?!Jo6x5&ogVLhyl*uOnnW~fsb2-NzaJ>nz=G9NuSUI=_Tw$ z1NnJTwRQzcp=;-#5#Ir({=HbL2mOjVgS*aV7|(xOBq_MUxWGTGr+^uB@0p^;<_hpD z>)A#)@gDh^{E0Q<`_<(;duT@BdQln!^hW_Fwf#@T$KX1q-`e{`A>Z0>=^ycUgjlrMMs-6oRSy zcLEO7M-%M00T-YAcQ1)`G%A0obzTHMe9}}YVu$gaT1mv(Sen(A_CFSu}z{TXDZ1*eRtNfR% zMxns>#u^vKI5_W=^1wU?)T7CsnGsYJ*gTc@25l5Cjz{Q{> zWE*jxN=+9inh(4W7<;{U5%ve(PQE4wdP0%7ZD8sL^xOO%Vmjcbwfc4b3c#6dvcQ)M z@-HZLt8gL2fA`16cf@z4s^iqvW5Ay}niGM|d@E16GSNmfG<;hp9TgIQjeTZADwplIQ7q7j3U=90U>*oda!gcp|j zSF#-ib~hj&tL9>B3Ss|AmqR6ZzL#)ZPq4Ft zc+CgaJMu#OSaf~<^*&eNzMQ)ydZm|_<@8adSAhP^ZA8x;`_aqT4Cme(64+R{HDPR{fu1d<2vHKm=asM z7$d~dFKyh%67pqqR^N~x>|~1gk@gwtH;SdE+JN8lPHy(~po_cT80ka+?~kO~-XX@{ z(JKl^%0cJPP#rx<^t*16Rcr(facLHP#JU7#5o(+Vz`vQG85LsQV_MR={%PpPT3$xX z8UQ~1qgz9EL7eXw^!<1W`tsmRdSh-Kj2F+f-?;>NG4mlcd8`$_SFt-Abs7Bboooo- z4ddHJ8dQcL?!Dvm^dEo+^^6>$#P5&{I2(o zH3@%#pOFeCAz*iMS`M!t#992ijH?p(&HHfJi1;4kl-yq@S>S`>SH{!dAun%+CsAAm zydzRA4@pAa*;u#chXGz&kIPSbK;G>Wikx=!WA= z5|$OAU}ukcYG4-Fq3B`bFF&8{TTY`LG6*= zO;8%=zfzt6J(<@h{bB^X=HC92@vi}1HvdkifYRB@n)?I9b^1l|RT+r0Fx_uzL(mPR zw5^L5;@}t4$w;ghz!W6tK!af%#d+n`CGcw#+m1Ot$iMr+XV(%T|73o5jOT)1x4U@u z@qk^=8^jDJAP!Cm-M-`SMw~H@%zXrhhoAGC{KwQsMOtUCLyA0kcG^}9H z>Y=U>3fQ02=>Ak0uKz=yA<%Xm^u_7K*)ix(zLtbJSo1)Cx-Vta9JX)J?$e6`oLL<| zIZ%KOT;b`gas`}pp4EHagZi}5?$T~pQ0kRH|b%3k`WE~*u09gmfIzZL|vJQ}SfUE;#9U$ufSqI2EK-K}W z4v=+#tOH~nAnO2G2go`=)&a5(kad8p17sZ_>i}5?$T~pQ0kRH|b%3k`WE~*u09gmf zIzZL|vJQ}SfUE;#9U$ufSqI2EK-K}W4v=+#tOH~nAnO2G2go`=)&a5(kad8p17sZ_ z>i}5?$T~pQ0kRH|b%3k`WE~*u09gmfIzZL|vJQ}SfUE;#9U$ufSqI2EK-K}W4v=+# ztOH~nAnO2G2go`=)&a5(kad8p17sZ_>i}5?$T~pQ0kRH|b%3k`WE~*u09gmfIzZNe z{~zf9B?TP?B@yv|{X-Nx{*?igl>h9Y`PU`|#XluA94C&j@gNHNe;tyikYu8e+@St< zJzCgH9H*kB`d0(=|J0|ZqPn#E|7`N#r$RmUe_A09z0W&lNi}5?$T~pQ0kRH|b%3k` zWE~*u09gmfIzZL|vJQ}S;D0q8_0|*B>M#7a`IF!r@%LX2{(F=BFBt(c0%Qcp2#^sVBk;cw0^Hbo$H=J&^u$Z; z_4+^xGJW}E8&jeS;uhOk(!A=AtK;tN66%abM}-8sKV*6$ogQl2?gJ5sf#tzu;fgZq zsO)nH-4csrcbJE2zDP&y0c-ek2CJcUTGbU$VOYzF2wdt=b@`RJ{pg@ zL|`v}+a1Gn$vE|;oA&Id9JIJ&T+W+06;=53a2f|@qj!#q5>5;`=xa@Gi|BY9W{(uO zsQV!qbJ@C{jeEq8zg`*mRXyU47S7TO6^O(knE`uqH>Ciiqr#-8zTXIu&h4HaDNskN z#?d&TKMrlzd}}l-K|<~qHA44fp1}jxg?c3K`y+!l+f@4dO>j!N(ADuBH<3W9S(nG= z6x1ovyd(W632CUbIu@)3VybBeTT`|GoT0F}@bYL1N;|=!5o(@*nw~ar@aB4;&31wg zc^BgltIANdwpJ`Uq4lzVD_shbp}tn!yc&ybrc#4LMpAM2F`et*VgvB~^^MS?pV01S{8htBGLgGW;Ou6?3 zp}qs(Qopp@V@F-NYIg%CoGn0W9Ixw&5>ukKJ9?i)#T`>Ezt4x_6 z?S%xib8RyP`$Q^^t7P>3X_164^YV!^mjEGWz!BB8lgUs@56_V1L^21lM@%&%2uM&02uuBSe>MLi5oZYCZLGi!Ah3 zH1OURjYQOLI3&hW9fL>Cisf{Phhpb~NyV=tw&(zLy1^>ez|5rWzZW%0$c{VE(erjR zTC+d>oFwgtlGkhY#~u$yB_eXS8-l!1qt_*a7WXhDC!p+F<98eX78QHQ+U1HzCf#RH zyAD$Cmlb(1Z-whpullV9C!lN|GiBC+M8v}SL+b7RYv{|aiY_0MJNWwFaIRjHaMTcb z^rgFe9LoK^{MsQx9lt|C9!YXR$hoz&aBHwEI<-rFrA;9OyUP_Q^{e_LC$58!B25gj z_|x#!SqTkn_cMiVFN;5_56Qc4|2Y~vzy7^0*qng>j17hMZc0MNU&BQ1Ii(@dAV?C4 zI6Ncx`@p%1blmV+!7Tbr6em&a?3H@b0A(H`l=k&-qPT3$0vG!k`E%|jDecK_g zsg`Sm#J%DvE5*`qR95S$K(TntNyS*tnwEmayf0|7lq6$0p6@ni?pa|fyH0m49D*Hx zK8{mV3&GR2!8{dNN!auJWyTxsS?FZXUizf>xoCyDcC-GhFE+Xm)b*w@1~2}~mJkq( zL}Lxl55_hpqr<}9CJl=b_(_Js&Jy!j{Qg+@^vmLS9Ls*{M*DgMu9U8wGW1SG4pL$T zKLg{?{#I(aEYB?L?GS~3EoS0dUBwq3=f|TCytd|66o_mYJ;IzeJLAl=9}*P#qw!Y# zAjLxoBz)Yhs@+-Y8fJ9II+i{0*qw@ZfA8}woVL#Oy~xQ88wb~UXq5OM>O0R4K8uM& zhX;MX>=8J0(#Q>a_eDQ^i$dwcRkO-(2gW zl@ppYJ~%9OHWTS|C7Rzl9E%?OZZUfu7Kl0a3Q8(ln_|hqxF@j|nTSf%DcIF950OrI zey|@1N9eg^jC;Qsa=AL>{(w&&H~+l0XH@MDsw+A-d3rekZ3#F&wC!Oaa**ZEUS~AJ z#hX%W<@0W#FSkpXY1tFe*fHY-8sltqLUB>(n{Ehh*6b|$`oR&G6mq_6IO&EQsP9ir zKC{M$7+Bm`I^vLLgj%MHrZ1Y?_Bh9~(iYMFc;s_6Js1u0@ca}^O+XdOZ@&KM2t@bT zX$zisxnpf5s@rWn*U%~M9|PNeMB$;0dCr*>KfK+tJAkAeiLT93KYbN&2XPkPSNUQQ zg^V`L*WxU~kdvtD6fbQ&&I#Pvc=1^Po(Nf+Dy+JUe3%Z2MU)2OVw&bc-3oi#s%Rzzbm1pfzJ57waUJ-#~!H0#`oyDlL|i7G^MWh!4*~B zR5fT2kHZV@`*S-!CgE3>Z@;y42jGT<3QsroSbRe&{BmXd8QjCfUDbMMH~PG}<;&6f zX#9v|!Te)59W&I1ai(kTM!&x}f7m%1i7v9a$lLS>;!FvqYox3MT=MjAhx8?|+e75} zE^24Q_xIsSyJI4Hqqz7|;$9#+{c*Icy&)J&>#(O<1|;DgZXeZqdxCHTqv6WwTTUq8 z5KF{viCA28TZDnb*%jMN=NO$XkHNn$9+_;tnt?;_fL!JmnjpiZ`rQiL6?P0uHSxgfGPw(v&)xJ z$acqEiJ|llIB(;CP229xzjr|grnZ)rrUoMxi!m8?NggER?{3jtl7jo?XWMgrM&bLB zQK?~7x~RWQ;Jfe+KP<-3WZ+_*gKexR&F}Z+Ve3}@w0Xl=e|^@#UScb zbkuTM!Y5N>!B74*Mh^!9DZY zh6c8n@^%omZu7!#7Goq*3!;&d^n>QzCDF*o*Y#=X$vi|Yu24r`lZi5SO67LHjKx1? z>)g%1N8tUTv8s*!d8kpF&U#F*0F}SxXdiV=!0)eJ%+r09fSc+>LQfQ@pm;}_o8@Z} zsAr^ad0Z$L^_<8}cZ*F&bhlcg7evCa!`V1*;RGL?d*E%hxOx^kUus%nH=Tx9+-L56 zv?XDMpRx{rBuMx{+P6R&i(o8&tTa^I$S`6!OWLR%LjH2m)3=Uo>5D;)9IM>-4pSNUQQgL18A zUs-e_QlLdp-5;Z^Rbv!R{z(2*j)do34RqBv`v6-M3AY{=eD7$Gfj57hsofWb z%}ZWCNTxrEJ+u}*sOVj=@3}{uvtJYNN}Nn+&!!0MyuJOgky||evcHh~^r2Y1Q_UeJ zc}x(~o;SbN=IxIKWiI9@R$askx}6ybnf_Sz5bgCmx+E;ta;7PN_AcHpd|Ow3#s?J! z4G&~CJECDxRUGm(1@n01Wk)T?;#anz)Vp44<2u9E70WIXz7)B}LM;=A#V`6s#M7r^ zp4`ybqT?iF#iv!EC7FUqE0c}uuC9pwMv?#Au2_T_j$NpTA>myibw4T9<2pB4$GZ}ng`!bz z_fK7k+nMX^uw=aP*K-wJ~cxw%NGJ)rG=!53Go{zSl%|kXRosJ^*`ABwxgwko_a0=aSerk>sywFC|?XQuK zY}tP%&)vyJT+5Og*Z+oK&s3$GYO+aK+u`c7-x@i{?sp7*vqu)Pt_rzfcsv(vio3PF z!XOWwd@u7<;9NM4d(IjfR1%Bt)(G;Srb`v&u`PVDh zgN`pU+^v*bRO^Ph-&5|Jx*vvz0!&&D(C6S?J1ugRwbF47Hn)4@^AnUm z)>W@SeF_>&-AH^UoXuh%DF*K~mzZmM;exTri`JO#5KPy23QyXapbQhY=af~Js5(*C zA|yv1@BgXe)VNO%#qRG*iRVee?S8bc5*!k7+~v^YmY3CVB;AwKhfaB*V_mzF*2YMf z%Pn!kxilU#@l$>%OIAgTKN6y}2ZC@@)!zHpT^w=nXU{v+Yrbeu!uaj8=n%}1&|va-Ee*ftv6M}=)2q?7vrOe+h3gL<4)LVodl`Z=XM`-U_SOM2+2BC@^OSs3 z`(O_K_VcDXol-uIlJywy!Xe1>=q^p6i8w@8P@Cp+HUVpCF7j|RXW{~}LD3RYG=3oe&O9`8s;S7qH7H?Ak6 zke+a1o};mNL~Y^AY0_<+mEpI1KS&2hx~|@-bo56SnW1O6mSfSkzt1%9EqEZo13wwJ z^hTk`ijJYiO=)QF^zruBx&bJjTI=B#FFSlqc;#`)@kG35q2b{A&|Td7;6=X9hZywg zr^x5r0$-$c^7x=gi7}47DZJIV#UH&AQDooy%@?hxk8Jm1xr#UW@_k>Vg?zJqV0Wx8 z2_Ljk;#$`WMbg$ZJ``ro$Z6o?7YpVDd@1O=w9C9dR**O5FXgsEJ}k*HW}&z6VpRCi zu=+H-XR`r!=&fk{x$TwHW_cgPgzw7L8>nFl+k-kn0?BxUOP{ryClqU?v`tcI=wp?f z(w&)3cd?}SNpXL7UmWKTUKThwt6uAA1 z32LwnE`2?vhy6ragsKc8ut%y-$rW}_6#UI*-CZFDb>a|K#gz~=J~wjTUr8HFzObQ6 zedmvoH_~sHa$ZL}E_@w&JQsjCrcRsLkP_g2J3d#@(h8~Uc_(2(AB;cBHoaezOGHAB zwAM?C*~oEWvKD=~fpeXjHjNfWkA6t#ZbFPepvQv7vB0!&eK9iMx$z+J4Sr!iDf zV1mwd(|xIeeu-zVr~AEe>nGtk#Y-F$S8id358Jo-{62-h`PVY+)$zk6y#>~FSs8egC2Ifg;22CH zRQ9ZF*&h{eBWWc0DWfkJ6$K@)rs3XR6_2rfB&^?4nkRVQ1u4H>s(YSegUd<2Oc8Mr z=%(K;U+3SR$ZhLFb@5Cn#&r=o9j;+`H?3qdR-H6aZyZ&cEl zc_k*Dgk{&?AN0=+z{c{`i>^hHINdmTJfATVKPNdI9R3lDf)5t%(Hjmx_C->^2jkT+ zyUip!MOXm3RCvR%1IFP|mpFIJe~85_(!Xcg^fj=TQ9S+Fa3DGony;)VZh+L8ir+u^ z;k5R(iHJ(UGcI_@6HPUy8uJ#C(AQz( zRwg)qEn~l!m#h&U;9AXi9~y%+hTp~A>q$paJHsXgHs8VlXBX_9jbpKY@l0Dse<*r= zb=RJg0m-Q2<;;u2g3d_XXj3%@w?F>#I-rV`It3RUcx%%5JrdJqNhDo{exCFBmOq!j z#iAaqCsFy7QRvOH`_=ZOGbqIM2NT1ecvSU1Y-E#3Ac|7dOZna$h0Ndcq>=unqGZMX z1;yL`SnH_rfmbpxjt&gch}vp}UmB~bwFj$XYP{Q+lRFjJ30PNA1jQp3m#|wWg}m|4 z=aE7U+;_1?$HU8a6_Sz3+OOel%i(C6+IQ{dX@5jZyOOdxejV?BXDIb2#}Y?4-g%W= z?28IpbzhV?c%c4~DNP4M3v71c&!^yNPsAIMgl~U~L6@7Kh`-+HjDsF=@k%|`MT>vS zx!DCma4s#gG5u8%Ub>rWy*?3*>(q)L{@Lw?XG>kj*I&nABZvL(+9zP#nJUARS&)a> zqIx28PG#dqea=nIt0Z)*S=iUdF%<0xqdCJAn2q0gSLYkjWaFw8DyyCOK`4f+#QO$^ z4bm9zwBq`phnYPC-=xQcVIPT&rD<+I+{t(^GO)!4nUXFtT5XL*%dz>iehFTv*xsYC z@TeYovfa+=VCG%qz28}JWi}h#BMnF>vgV>$1=F3gqy(Hfa&9ALG8M})>~PR!&PKH5ZA(TvWyy>!(cTgr7R~H_G0M z!A#3@8JQom(U8R}qr~z&Bx2?0e@8P6Cm9*kT|JhF(MW&0^hyjWu46h&!;_8z7q{hC zos_}WQzMLhUxRVv2mGGz?q1a>BZ*?6uz;XV;Xmlw6f6P8+l#?2bJ_}QdO$aBV`ICnq-3$!GhC4*#GdLpg z&O2`cYbm4f`n#9@yQ^IA(2GD)jnW;Qd4m6|MPE4fqg$Mq5gSJB|r_%&ya5PYMBWokJ&4tuPq1Q`{&0AB@! zbk#Ufsm%4ynyv9TzNaeoiAyqGQdnFH^tp>cCTr`rMoD!n!|5t2yhMOe#)Nk~X1c zPebZ28P-Xh2`DehnAg8L6nhn4gE z+i`CPJb=!4J>bp2Sj?C|rz00Hb?%|**_(sY{)jL2Y)`?dfk}<7Low)P*I}m{<3VWk zw)wL&)CIV|=`Z#d&&5Q98zieb!|&|K2}qzcDx7A@ye1vYq*BS zp$D$b2I7>_Ncs*UF-pgy3ZbKX=30q}CWK|rwTdiUgqXgtmM7q)r*7$2!*XyVwPbeF zTna8oG1#|J9D+Dvrr5P2Zz0Z;vZ6B1*_dO7o?c5M9p7`S8yyXZMk5})^gAOWP$SYm z`z^@{y`}4YNXvg0(Lbl3_ILBcnp?MCNN0_}f%QkUN}nD@73|AzGyeGEh!a~yf1gW0 z)taJO&s7sqZ;0p#b?rE8lc38H;uwIp^EnyLe+ocR3Wn+Zx6YuveSVD>%foTTla-$u zdA|6;N6EjrY?0^!lkglI3FIcv`og;X!aXR-?CT{E`xbl7_ry4j)NEEj`S8{5MkZmlx9A8w0V}@O!wf zrJ55;d@Yc#yp@F7G$uSkAIG5c6otduVveY1oyITboIf_+7FKsG$_Ouab$afc@kOa0 zyE<8ZxME3xr}ii<2v;eu&cmO|{|)0QlszkvL=!10l~FfRR`!;?_s-sX zJIitGEuoU6l)94!Aq`Sm_DG_nZ#4X#-#_4-*YP}``~F#b9bz#q}Kc*8C*U}!XHER>j&oQ#)63F}V7@e0b|;K3|x-s2vb zO-lerukTKsu6dw8c)z_WGy!E#i!S9*Xy|4jJ>Viz0S6UbH%YTq!0&50HL7|A=(_A0 zAl_Ys`|kakdg#P}SJU$fIx!U3{f%m_@T!2=GY+jK712>NMd-!jFx7GT;ZscH?Sw_aT;34;m@5({x za&>iHPdW7b47+pqOBQxtNxnwUB;lWNDYt#Y)lmKM_{VFr#P761ojQ3w8|}*dPJTR| zf>g=av6$gh5Nc`md>N7oq3i0Oda{r(K~R)0VI~VZG=80#v&jd~&t@c^ATkIQ?ijUE zC*fUz_EYzJ;*n{9vLf$}RCFxXu7pVVE8BBUGKy;-%p00-RP zXwN(Ef>)NK-8rhEcv|3!rR8i9K4LbyPwJ=QwzenVb({;~{ygthJB31GA4_oJCg#wt z*ZL2n)QU0gx<%8E^jzRASe{OANd#s0W0uurMX0rsll7W6A63FnsfO>U0b@q2Z6i+& zcnaS3K2(s0hl*XePMH-T>(bGI+}|~zDkNMtDOm-F_6SzQY^?-2QL@GvVvbjhGsnoX z(2-Anr>*Nw8m3sq?QWtfz54&v)LedUf`aFBpE6jpJ(Wr5xT-wj(|W}Q@pfYZKjs)ge{5#e}`jS2(Eym z790wp-l9a>XL~-}EcFXA_6o-Oa0|B7Z7DeK{NHr!!94hKaI&>JD<8^3`+R$uVjwZK zu)yXa1vtN%%1;}_p;u?rs=NaY8~LwMciwOUosAvo-T#vC6;{ci$f}+10)KIPnznRESIAoOL`Jn2LryJ?8 z-Sii)q+~H%byPg|Tp<8PLQWQ}N(AD#!f`#$%LORDMMPL>S2j*`-*AXED#jSTCHKrE zIx=rQtTM1A5u2{=Guo(}RFYGaOz%HcwgeVG{% zt6YVFivGG$vI>2Z)gpX)szV$#7)FN_rT6>`!Z=5l17*X(c*mPO%HBhU9d{&{d>nIt>iWxl zfY2+960R9NI4ub))iNWW2!Ao?bd3G(@f;ZH7pa-woCd>hEMJiFJs^m=Gdp){I$Vts zu=}k@27cb>;gW=2dvb;4dE-{XUrWfLh1XhxzT43Fjx{^Hye`tGQd%FLO#AAZhiU-l z<_jOTDl$O&wLZ2+(V+QDbI7?}cCe7}{(1(FKe}x4CjIn(oCwiK>@| zj8QsxKd&p0m@0&T-dVBGSsE;F|8#keN-=PJ`}L*UE)(Rc_4z;j$N>4~4Vx|g7Gd*g zo0p~BMVR_=n((&iD6SO#(%PU1mxg|x!M;Ga;vu2=>#s3hJ-Aa?uqO|V{O7F`sl{mJ zCWUmP0$5*i|A5Wi3c~leCV8A@iJI5Xo;|3WjibF&>v^Zc!AfwhX*Zn$FT8%83Z6^C z6}L6vpmS98t^ak`nehGJo8B>K_9Qr(!jO|i1;r3D_};W*xc~y3V}6*wEkSN(+l?8I z%hCF9aOuIuLikNt=Xte*0!k0+_at^FAZE=e2&U!X;RiLk?*iiRz!oR18<(R1HOJE>N-@e9@=VW9! z3Oxy7$_9tMp$5hZD+Ghqld{pgpOTOZN~3 z5|zZRxu}PMeT`Jgd!-9FYN_2G^~e`2WWJ>He=P^?B%dpPWayyLG4*3*-34GjSD|3E zoCOU5Qd)V96UR zm$^y`&gOi-`H9f!+6o_%)>Rt;U;9s1kHiGny}gj`yTuUynF$@p`!63k-8kl*52xdw zyO%4EmlAW@BHjJb5)CsGsv-qn7O>(6vLZauF17WNMK-+wreum7=PZJ z+j=xI9`?EI*maIC29)*x%z8w);_M~+J4?MOxU!6frQwQjc7XYoi&g-+Ec z7S?R}xJCW$iuT7$yrI6UNuL+4Z0=OBE>a44F~9BqjH$8>#a!C`>Dm}Zyk1^GXC zQ&68Q`sCqe1`Z$8Ph?jL4__>_5l#GFATI%q93JVWB)uIksN zgTiz$X3yW2K1=xX8>g4HbClsD%`;yWqwH~HcEHUl!ZI%ip5-=* zPu6r0KE2tm=Fd0np=G&I<_aknUgRu?#M>m|0~e-Ti|2U=%~y{LM#W>OK{JbdOu?O9b5m@9)XYYGC^EGo$@69`Id< z9ClYa89wW1Moml=z_+Q`=kKC2Axl|dlJi|O?v}0dpFmq^^VMq~i7AJ=>eHM~<20Dn za{YVwLo6~rCKhrYVa1qYLCUT1T>LG#n{--|`1~1s8cOINX1hGhbtF|J3VH{| zz*NxaO^!nexZ6nak|iS&Pw^ewT1@ztL0j@%JPb2nBKZ-g7;6?NaW<~ZOyq#7q{}+` zSu)m(YUs1PZvKPH9D~D8wgnmpQXl*(a=HTrGRg8Dc)S-dbJdu z0?SSLVe1p};q8EklJNr)et*)<$2D1u&hlN0F&Y5Omslxki`j7Qq=DRD!D9ID`;n_< z-ZXHOzj2Jj5{$H;)=&Ox%|Uk~v3m-kG*pl%rX1K6fgCL4`ISQjnDju{`|{;nl+Nt- z%u>rjg;mo9M~y<%DEYwU2;REG$%j}mj&@`a<)(qxiv+UKA@7FR=vrs!G`)wA^PtR~4koQ2@ z0{gg_O*AN)EK0lKPUtjhGy9fO%P__GL1FfEKEBPr+sq*RU*Uf{v{kiiaSw0zE^EC~ zJQL#HxiPZ}6~491z70!<#=Lu<+Y%U%&F2)b<#0Ye*e50Q%clZq163XDEi^)Beq(I8 zs|Ym0PoCTI<~)4Z-Bs2W;ej(uvgAf;AwD^EI9$1$f|B#!y4asZ68c1#(058K)PfcR9TD$ch5o^Vc&@^JasG#}3-Fk!lc^DVFevEr)K;$Ac;X zxj5vdeeI%n3ij|C3Em$d_`d{U*2>=%Ak#f|u~sx52i>~=nNDb7Jzwt03tq|4wpsgq z|2jG>zj~g!e}g}$u33Ex<&Q`ITc2oKpC@B%gR_2!gfXnO?OaoI2t|$V*{34UvN1vK zSnaA&1?Jm1o%eOEKc~n=4T$Hn>f`@W-_q-e3$2sl?t@@>C+@%Vt~buRNPj?vmhkYFnQ=)1zNMU)J`6) z!S%wna|5NlX`H9qY@%)EA#~$7`@2jg*^u!f>D&iq!guW0!o{7Yi>tejt?+s{VExHg zMv1kN1mFC}B5q?Mx}V+kQ<}(XNd9N@V=uv1d^~!}Nl`QtIs;vQT_$oBIm|T?Pj1A+ zr1%f{#3uoerz&_*m(YW*Y0n*D5sd|L6NRr@!#i-Sf83&XAPcUyda<8VqQIR_YW;@| zN$_1W{)EXfbNrC`@|$040u*&x-YGPufZf81;q0bZAlbKE_0Zgf3o)CERE?Ep9 zb~6J46#1U4t)*g>V{z&YzkF0`I(Kr0Ee`6cemvLjra_{k^NzAp$*3?cc-}ddf|M)I zTKa{+gRT{v&^CFsxtj!**WJG;VJU`c ze$AWEtim0g3o>szQ=yc($EZ>}AJztTXe!(w6V!pQ_nc@6-n=}q!6DTTm1Z7!rfg1v zH~qUphvv&sUp}Bg?R7P(pHw$}YZwkUs7i_~%_Pvb3IDhMoihY1>Fw;QjK`$#7Slxs zKa3vlPT)|EfcOI1`UWWqGT9zv+#q~d;k0iJls{oW`R-*w3Cn=>W1;M7#S|DU=OAle ziUP3)k$=T!tw6AFY^f(b0Y2{j*jx960%e)s1eU&*LjT!+p(7rSP*=<2-(?npD=*Iv zzaV@cm6}aEitBS>V~Ry$YiI?$^XJ>;m6i{kI~gu+4TKI|`Q%rxO(Zhc18sp0QPXmix=xJoMp^u%@Nf1iy2M|w@+Dd%#a^YjD;^d#L(cy^crv0IZvau;A{oxOXBuv)+!7IIo zg!>-7$T)67!m&pilr9l_)$afNTXSz`L2jJO(`dy+SQ=2Q7do5`O`B9-Tp&1~>(i<< zNpc=O3Hou$tT+QUEq@(*xo87#4k-@XIYy!k)ActUTLZz_N$ih8P&`~4Wf(0!(FQe^ z%&kJp8Cbi$s6{?81W!L$VGsPChHW>GI9@I?#-WoA`S)Yv@bB|cu}VKOO7CCK`)e*djr%VHi2RJ(>eUlz!y(w}=(_Smnv4!n{hP2O6~``k=ojR= z!uQ$`0XK>>*5BH`$#*#kr^~o&+llA&K_~l?>M-K0GV}xa0e*KeoU-W5#CE3i zElN)b|DRVWxV)weE6-6#Y7gm9Hfk?YYRtOYqDVD9__P{dka9Zvi_6jJt~&U z#E-fV`4v4H$CJgHBp`jY6BubV#7{dt>3kzJEX{p>ZiMi;+a$NPo--(d4qkH$?QkmG zlO?^RXcgg}s}o_v(o~$u9p%aO3Io3m_F0Fv=0`Waw`%57*pm4AIeF}_d^@zM=A>-o)b?If- zc=Xm5y(PP;2JSg&F~{tvfXJ6v$>)`iGPEm6=65KVAJcHu(ykV|SGIO?tyF=;eLsou zLMry39&`>{p%6NJ{Pt;Of+MtWuK4SgkDsWdiP61PFl06Pl=o#Zl%L@&5f&zLav_gw z+IktF-QK{B>t!WGJn$zS$u{ za1^2Z>K~6}pe^;Lf}dU_)~F0u8#(8pQ1Qj(b4p3bV&AOMni>t);(Zq?{9Wh2h(Qqv1#nT;o;&E@l$B7PGMPS5zxO`Y31-s%s#ST;^;Aq*Y z0B#6?@ppZDI^318LgnvW?fN(rIhySlJeP}s4Kp@uxuICZXPG!xbrS7QUo3&yvxp`Q zcYUriFw=RYXH_E)ivshxx(`$FC)~6A+LDK7%;`<5`D&O^6#2YNo_Ox>C%n>ck4J~% z=lTXFG<@jF%PhTxge>mM204N{n0MPS*KhbNoQ_m2xpO-TU&nZvP1VpbWL-G@#GVx7 zI{Ea+XkQ`-tFk=p+#8PE^n+#{1UIw((~VH^XQgC-($22p>%sKr|@(x?3OQ;Hu+uw2Lf*2 zb1W`^wiW+SPK^{CEg6oxY!!fzebdBIupIX0mtO2pD}<4vVNR1sXAG^|E$~hx9GH2+ z)0dsXfZLA$VZ0&zpGr&APTi!VSEm>cve9O$`JK?}8BpAsCKCCJ z*q4nSatDhCqYV4u&d=ci$oaZIZr^+|T=!x>w)$+TRkMk_7r?nZtpu>Z1`c-uvRcVD@W`I;2r ze{rXbIzE=*8_&lF)l$;ohu^CJnJ^M;P>}b_G$4pYu}!CVD=P3t{?u#t)f@=k&}2EO zk^wR4-~LHCKs9qj4&fH6B zb|CJB6T|g^)s?t~`{gORR3;?YwH+M^roi0~H68kSH7J-nzf)wq3UAoGa@(Lo0=}9P zU&MrRVEJf^z0EcvciI$hST9CKzY{AeJR=#v;gnu+N+JisIG-i=GKx@l!RyY#X$Gzn zN;*+3PVhZDB!GoN@J1og0kYkuAY*iYkb9>Is!Dy1wYf@xEtij~4CIsHT41=;ked@c z=sL(z-H-%LT6Uc~w?!lCtBdm6fBNC0J5Muj9?b*mXW>Kaq%3d|Uz~>Y8t_%PdQeES z8Z@}A$)D`#=w6@lBI!v!4)xbY6e`uii>9C96{*z_y+u~adX9qgq#`eay@j~DRyLU$ zPyw?CYfr}zdXx+A(XH7b)!-?#KIE2s71Uj9;!3h(A{2jk z7%>`jzZKya@}Ug93jUVX+{(lM z9@bmi1f~E^3a7d$ctY}VW8vDHxlos@Gb4N_4J=-1CZuQ)IiT*ie|(3@NZ~)syROXz z4BEHu&Gy{_98F`FYk}~o*EalI(8(tDG`0^tF8;_q(DSXEQwML4_V2NLPH>S6f+@yE zWQ_dw@)`40XDm%0uxRO3!kJC|yH=?*WGI|y9t@;n;ppLQIfEH^WV@q^aW9cqx$b&o zlc6sN-?r0d6@GEAz`gHPwT=x~@?3JfGo*!i+Korf7ZjpVee@)4LmqlMU!aV94}#)1 zPYtTbaS*0bTuxr2<00*lKPu_@xHOzp@rD+Ik$Ysh#U4|!iRNtW80Q05jLBNB<;g(q zAL98c6M^b|COcJ4BBAx?xflH5k;v}aZT6ot3E@@8?!VO(xO{M^}WB3i@ET58e44gVvpA9Qq^u{fQU`K$@7pv5Y6?@hU^^BIN|=} z2;rwjSDGE*oUQ;v<=Z8zD-2X|e=m8rvH)$xoW@9kME+I%WnOYu4e&!hS6xsZ&MTOn zR{m9lef(7}c21QrnDVJF^8vvp{1=qzt(S%?F6Max7wNdCrM1w?%LV&ysdPX5PQt9h zMT?y()!-3nG&l9A8k8lZ%Qi~Y!uEvs{gboRV6gsBXO$lXoznZVqQ6t|(d)GX*O@E9 z__jGz<`|l!PTKH~Z zIhKv;b&rFpmI$3b)+}NoI2ZP9_Ee4IqyqQ3P7V%}R8&ec=G#|9!qvqV;d2coIL|OT zd&{#3voED^uwKc5gc?`L2G1!xS2tYO?9D=t zQO?)f8*?%3LBRR7xon)9d`4By4#L`&`J@fs60veB=e;%Yd!PRLXA3{uIS31=eQk9z z4lnCZ6$&PY;rC~DcmLL^Ip4yY}yj%)}Dg-)nyNt#gn0D#psUYi&(I5 z*Q(oCU5vDA7akZDGf=kkTj-y4WN1s++<4uB241g+OV+J1FsWoojcsEQdL&J?veL+S z!C{2=tV9ut*pvktTIHg_RQtZ^vlJ|puHJcuTnNn94ZdIbM&u5+?K_!TM}<4He?G>0 z7lQiR4fO^3F}N_T)AVz50D9*sFzK2xKtpkl$>;aQK-YLDaMzzi_;)(xC+36T*JP^L zgN0a-;qDlUz+CA3q}u%X5e1GojGvS7iGZnF&-g>Xl0l$fSzb%{yQJhK=;bg z6on!hIJqC+_W-id-8so@g~+#Zjfz(vxfTG0!cSFv#u7o*u77xEV>VX*duyCR@RF*D z5*zzT2_N z_rUi*1h&Uz!Ftp4yOp@}VDm|nqhfg!km(WsJ>o&^bs6rH(j0PI#*4rp zD(diIVJgIqX-Q2TNW*9j-eWqKgRoR&npUhr1*<}*5mq8+5H#(guOJhIGJ<;lmYXj? zZiHf&X?-$^v8QJQlq}#${!Ff+yJ6eBEiT2Cw)le@j@Y!L987Y$)fcU^{;8IL9Cz8C=WxG?%7cy+_Nlt@&&`{%FXVgX3F z%isAKL4w7s^uyOR9fAF|wiN3N!e8o*2>lyGfj`fz4@6MAo}7k#R?@ z>!Up15ZJm_IOyua>2@W1$5x=%!i3~c1<2V>7Cp~tPwt*>RwVfp-1aXKp*GB53od(9yFIkvqv z2!C0PqD~a0HwoqVWkn*Ry`KWq6WzPM_GG|=t)A8sETy=utnh`i3k7w%@bJ@HDG<)z z#C%%K4PJja$#Yq%0?m?|!xwKCVDmt-hoXKa^li>wNaKu#6KwMLD|5PyMzqJfQw%^_#8*^vD;=J?2Y-!H95+@MsD|{8dcj zCP#qi$;xMsJN=+R?3bvd##x}=dDCtP&OlaQ~p|`QR<#(_3=!6M42lvKAHq|?zI^7tkLmc(^qP$e!6i$b{!UZvPzAn62OP6wL!obSwT1B_1fL4% zO>hzZQQKqt>rTWyUGO8{j>sJ1{RsK^fUlKIoa@T}E(RxpN$1e<<&6}WdJr3wWLksk zl8$(8N+He#o&rt3Bk}l`e{WlKy#QY5e_GnPJ^@pNSz|<TTVNCJpu6j3m-5F1m?H4lp!$`pmk3X9n z8I1<Tpf3HbBhPtrF*J>=|#s{{Nw$bQP6 zzpRahPV*|ve{?gDwa=E&lr2zX!_!)IzCiRIG`pITn}Kq0xzCy@2{u*B3r|lNfmT{U z!NcxMOsjvxoM9QzDM*qao^< zGxOtuT;Qm;J$GG}jN#ma8a>$s_+{Oh)}kHJkT7O*cH3zo?|+JK;Eh=(CTHERPbE0CNCQ8U>+kYlr_(NO z8kU2~Yk}K`#y#=S`EsI`jEui~4IAmieS14a+p?F?4M-O`n8S)=A#~V%^zwKNhW^kI znJ3Pf?)m{cR&_VH_2}V@LNo(dg139h7F2*{gvw!en{v==zICm-v6}El$A;QRqS3We z?%tg$Dh6eJ`&)9C4AHS|^(zh~kl#>pIODAgteAU?xs;@#*@111JaH9pYp{2N=z&sT z`Z;IGCQ^u3gQ8z&Z>wV-dFP6Ei5+KXos0RVP&nMrQJnK6 z=5+>5en#1^3MVsZc`41#$i`>86fh79X&M@*Ij=_GC*u{q1MwO7Qc&P#1rH5v4{{X- z{EWl*?bOpRF8V-hdGu`h&rqDr{?I>x6{ux$@!IXf47^u-{?zHTWEd=ZBK?}v1#D*` z?{n&qVe4+gnO}tOPODV@`+T?nx$33LxwA9zK*yh(X4XW{iE+Dpmkq&N=}PGAYe+{! zi`Lr9gg+zF!lf_2Ou^zNe(j8iXe^k?3EXkf5yk#lY)lJH#}24c!PIwcoyJ@f6{Iq%o6j4v8`_@8iK^o|}zv|A|+R#9t@D{7A#927d$w#1!r| zB66-Q-KTX+FIa%A*Ma7Mw~0irTu;ME{X}SIQ8diwr-7W@z6S-%df4=_DJbk`5~yj? zMm%`!;Xq?s1s5e0UCQWhBYUavWIxHQ&#o8>(pqaBjB}B*?%plsl>*G(X|gOu@W9@B z7Zf-6l)#ZFmE3-LAtPtx2SGi8?9{Hu5E(x$HG@ieI%y`+RGTk#XY;iw@~Iu@^>-Gnecygtu!2>V|%mKs8EZ zvJA?{mu=T}HC`=3IewjtiX{VBigT@)yd48(@~H;0t3)oLghjy0I1_?bkLPrD_+i3X zA!V;K?r3>~;6tpvC~Iku50+!?oPZx?KNmqyF-?8xG$>+r|;P~0KHA#`+NCSS* z<-90p64RFaw}r?usP?)W;ej(_dfNba6%$X(2 zKsPMvPwwez7U_vq(p1C}h{WN$d3oF4}r=Z<6v zNd@4wDByzgVbbb(95wprV)x`C|C=PG=DIu7CULR3I7+2doYX<|23Hvtt57u{d+@?;(L16DYA=Qj?3! zK^cqno_+U-&rO7fYM+wfi{1|ZGoJ}=pC|4Doyid&91}N=J43>C9jA?IHl%`amU61x z01X@}jz)3jrlZk0hniv^1|BD8lzz%4dK%Tsu01_d4(8?0H}S^%fSTui%Gcmj?D5^; zc%!r!$}GEmg8o**$H2okjm9Hzq~HGU1!6B7;M19^Z%=|}&9l0HZRqgnreRUawm9(Q zQ$EKd9g0bDrexq zvcw}7uLRufk0SKq<+YT9e~6w`^S@)7Cx7_j4jly(tw>^ju6|D$nNCC9IqEpM)C0FS zR|RHBW#dSKw`RzE9`0rI99z`&fUN?otf3k?ueldn6RJ(z{RB`erJ5qTITUY_vJlN3;Gs1n{iRtleJeO+PijG%o|GfIag z3k%ok-jr^m!?6EBm6qR0u)Nwe!zr5&>26#8TvacJg9{!S;01D}2CQpQ8l%-k;roV_EQ!yxCRXAOg(7lk@`flTr5Kq=54_ zV*ti~9BWo2Xl%jjC*lG?O748PcTfq6uZo$?2ooI87q6D?ID*T%eEdk^=WO`5(}nT8`Q#WiO-qoty77YUNt7{jC((!T86F7 zUr+i=(cqu5Zs@(Q8NkcA^w~GB3ctTHE2lD*V3$Qzldv%fq{W5%UUnqFq+fAL^42`q zzWV#tHZldW2HzA|$`#@qTgBOThqG||2N@qR<5KMLpqKd=6Fs>PYo^&wgh7J*{Y{$u z?)a{M<&kY#8M@VYC@mhN)V@Ryf7#$+Y^hBNpz~pT7p(&nJ z=rNqT;96}5DV^;vweEz%-z%TF>62vm7n8otmn{koZDHZxaxEHptb01loV>xxAz#3t zG#&U{f~9PUJYw~7O!k3rcPL#KLd|8efW6zVb*n)fJU`J#`*oiVTJz=hgS~0M_blwv zrpSDFus@AoGJyh>4#lE*;eoJW?daQ#_#Ak+kDY($DxuF)RoCA2=Y!?%m9NJ8GC+;d zxybX=0q_6Vv_`2qjW+~V)oQKcVQs?QK5=s% zwmGNLbSdb0@{4AK8}S}o9#N*Bk4IA`^|IrXJlx~u6lblQ4+ll(YO-`0@K-!WCGd_d z_MZFKyK!d{^hAZJm(WN=F23XBs*xuMY*9<%c~A7Iwi@IlaHpd3X{#|$f=dL|oxf)5 zGJ)AH&d8EM=vEES8Z^uk@xZ-kjW)s;>vOB(I=j&ywbQfQ>WvC;VxRj1F{01NV{zd) zuQD08JMnt#nI?LsOD^x$i^#`2nRCbOmQygygDz2#kOi6-^Q=BF{E56=tmcDLRe14) z9bU5 z4A=fy`6vhEn$muX6MfmimPymk*=kVysRpOp@e1q`rEZMXC`X^m3=PX18j8B#l3`1X zgYC~$f2fHA7TR_w-nmkVn>eYBEkqB|2NzeaP+KBr!|~M6$sq|idabIAhlySt<1cUc zRrY{aeEunOe_y&K%>(2?}7;ok;B-22_V)}TKb9=Bj!OGq#X`HG17X2oHe#|}Ho zV=CBx-?H2^E*{R$qKnJU46xmKt?1>f1Eh2;=6d&;z*kkZ0^xI}ICz`+pcB!vW|n@> zhbxZ`dnabgug{R6>b}JI;vx;QzdagfTRQ{q>Kr5b8~kACzgfYq_Ea!sJpZzPFdxL2 zop(s;Btv&lHK)M!K(tqtAAg$b4bzSWA75Te2fGa?&IXoJKrG<_{J_ z?kHK+D2+P;Tf0{7TThVS80)%~FQ1A*@8sxpBXv$x2{PL-IGv3jd{=I9tR;ecT88}G z=0KdJa{n2sq{4Wg%{#T#Vz}kt=xWSE?C;NXj(_G%Ms8WvsRim8)FfF>r8*FtFX*i1 z=rJI|2A%fUq{Frl(_`=Mr{In+6^tJ{;?QM_dW10HBMfNjDmEN*g0q@lzn;ap<3yQn z787eWvRnRl^7Z`$5LRx!w2@AO9*kulb@L*)1upB*{kd4oqB!%$jS727pZmBM3&Gst zrJc5@9q7fJ<_$DI4^-|~+VwI zjUr*!n1kqR?;sQp7ms!~ONG-C$BacT=Mp;iFOw#M7d+uSvpEZeYZUJ{E@9SFlj(hcPOcWy|0G4bX@b$yyB4BmT%cOzpd=8b!#qM6yl2IA#x6{ z>)G2?r4rz*@|f?J@o*d~XPR5!aDcVF}gM=y$`MIqo{7bjb32F8bXFd_Fc1 zk5sPz-u$Y}#2KvJXG`lQI& zT6zjDZ9Hb}mr@K%bK#eqKKO;@ZUP>=7L=bO1;sI-3kYQq%C;YVyLus4dv@4F=ZI=1;B zhf1Rg8(T6SHXftCaG)WBpx5|(aRz=Z_vf)OC%yyy8EaPdT$~tb&NX)~$9q$s*e65E zaVhn%aAZU}xGse;{Ho4^&DVb%VH6~sNRH0_tDC8~HQ7^{y*UFajpdk@E+#g&y zHuFDMKzQ`jGe%GU6(GZr>&HE~D`92eZwtM3HQbICRr~cX2amfc2zf0Lzqk6oHA&NQ z=={1{J=3rZ_P>{FnQ@__{E0WSS8fzy)}Cw)OTv9R2|X1SRyEM$DmS--mW+v7Iq&9X z^Ds_xQ_ohBTCm8!qp%~m2Kd(^Iy$#xBl|sc6d~uPRzpA6*a-^EvJUs}8ZCqvjWk;$ zws=f?`FV=5Arp@e$6RA$s)UD1mwyMF63$fg=I+3-bYv(Mx-K70_II~^vo?8Yz^@bC z)2o;K;HsXz_v^Ye6kc9G#AlI(f-;ItFNe~g`K3wr?O(|-GTK=B_(~2gt{ZqxA5Y$Q zc8mYXNOC86p0recHpx<8ff5lx>k4?azm8SQkb>d$ZudK_i_w!O>y~45 z0i@_2-O5F{yu0e3596~!9A8-P_fDKf^nrWS(ixJUWL}lASFD9hj|%&PG8w2N`0%D;B$Pz8uz`Dukx zvSU1UiYV{oBzcP?YzH2`8LtMPTW-5U<*I@5H*VlB$(>wx`qU~%^p2C!dR!{EBjGE3 z$7b8PB#?Xb%$Uw038ikXZO-W+dyK(G_T~uke7`(+dva?&B|CMe@(^FsN3hjZ|b^;>#X&s?mot*k%U=?doj%S&=^&V$gWbMN;y z=VK84ZsvDW`h5wyslZGO|czq!{y= z?@pS&$b`sI_C5)^bg-WE+Q*kq!|GH2#tN!wc--UTuFP#!kV(7e&uLc!pNFT<@wR5; zaJ@Y3DA_Om;yH4^OtTtZyp1_BC0PdL-W7(`f|U?EGu1xsTLF~8svz>cU^7FrFLQeV zT6KeYf3`jf$+7G1rMbfNbf=ak8wDz6toKCwO9U1MXJM%g>#0=lzsY=dDJtiJjr(0eTv2DyH+K zN#_Dh^w@=a8L@cS?WBQ@M*{kI2KtnC7sAO?4!=!lghOe^l$;fj39ZZOna27BU{}$( zT1`#id3FwVfkK;cWHQ|J8`xmQl?7ssFKvpm^AX3k zDE%dR#@i`Np5CBhsNk6Sa5b9>Tg0-r3_mKvzKtA|=M=(M{`|s!zcLMM@|QmR*=~f69qZ%E;{me+>qkYMQJLRN`5NzL=@uY;b6N zT1wg^K$^ddjp184y6n)Eqn1_RjA_%y7pJmdt;aTvrzsE4bMH-$za5J9dIA^v&UnK= zj^MhL>vq6Bc4bMjfskS^%%@z@N`dT?-HsP7+XAcoyRnACG_Yt<)}oD5VDO2=nZ1P9 zITke1&-Ir0-ZKA{zG)!2;Q4AstDpp!cspDcM)Ht}qD=9=I%NJydR6fwBo)@~KV;Nk z$%M(t{4a*z${@leBvwU?3I~<1bv&s_!WwGR;!s2&23lzQn7d@4esNxcZDT4v_!@UT zd`CGy|F z0hEfrN^CC6g6g{=)cF+$+<5L}vG6l5%w3xbc#=njIIVP5&t<~plDx|&(;fm4^O#fq zQx2F)g$y>{I}ZoX+>ra?lZ1Zbzl5%uCcxx&_6>z;rf5u`sPQ5-97H>rMT*#S2wy{` zcTPK%4Pa_$A?Hr0w%ku&VJa$=Z?SYqDZp+K(Ze^~XmCBN zP`qNi2>ec3biaI@28Tr$SjIQR!i&PaVh3LpW8$4{JID7@aiMo&D~l@eX+F%np^;0& zxk63zw3}x@$JcA8c0m*#2=2~zx>$kI`_kVu6F&R({wuwU=ZZl5L95n^Kgo}Unyc~o zRpSNbSgYLvHE93+O8+0WLdZA$)A53-7!==sO_;JONB+xA&Nq6h@yr2T2JyY&pt(Uh zZLA^}nl|-oDP4BP-@&rw@i&v8?%JzU(#q8+^BavtxoS~c_d!;TaxFg7Vw&syR)rEz zZBn~M^WaSSpq|W74n!OXoy>k;g??A7PoK%C#`rMF8OJL*M6VuJ2{xj_cI8WVbpuN= zn0x4@J@Jtw_t9!zeaM2et{dHj_afoP1-rKmB!??^dqKpBeBVT6Pa1JZz)d-(Ch!qNVxQO1H15})v`BYuS zmj=NTsVTzbIn4O!a*&Pi#;)7_r|@h)@mV_Nwi*-?9VF|~g-UW?$vfmBphtXyDnDW= zw;40>>ENTx<-Zl6trL~=eQPQF_H!(FpP!7fCZ{$lrAOhvp<5fZo3g;w_vf!b`FJQA z4&dB0WslTOj*T}@d!g^ULD7Sx*P^)G$8}v-It-r3Qr?}D1EG*Nxy~;Pnxdmv(}qJ~ z;jO&+FY-Q8JJM^rl|+Glnt{g?F8z%Tk{>()#Ym}?mLqSV@_Qs(@O>*}98s2oy0q)vQ37Zsz$ ze0(7F-X(%QhS)0x zQiTb3BT)G5oP7{_?X9(F=`q9m^F#F(4`{%7;i8GJeIBfdh+W}hOGEAL46~JAGVs}^ zaqGd9BH-=c5n%JO0M2F>UjOqn5$;F1Fsi&KocZ+S6RyJK`*7k{<8RCdr=g-wgDfhx zeY8K>O!N!>M{{TCM$^$XjNY!fy8sst&V89|P6U=i+y62S#ew|!x)}9v8ovMDym3sV z7!|%`*L^2E@gs5fb#yz3KFA+Au3{Mqk?LE$G{++G%&*$1@E0cJ`QQ9kT_X-etEX$) z1KrTDSMjm@a3&f&Y99$!$-^V1QSGcl`IuAI$yj~12sXRL?}+mu_ti>;!(YEqact(; zh<9HZ%0@Gbud6Nwr+`~Mn$HSAnKHgdGQu9CTX#i%9Y`j8s`S!G29lTA_a?aGs0r#S zS6!EOI)jbo{lxkZh9!P%5zJeIab)u2u9|Bx(7riebM9RVMsmj0#uMHC?`VEPFg@AR zY<=yJca-E!C+I&}`jR}o^@WfB`THZ?sN0UamCTPW!LgQ|~bx{t^m|ezW zTF*=NdY3APo{Lw4saig0q}D*7ELHN!=XmT@myg$h z<1M(tO20Je8Eqx;D@ z({cMqvqW+%~yjlrIlJ{EsZ&pt#6fNBY>pe_7;)0DvRHXmKRI;0l$vysijQZ1Rq8pK5y&-*h(;}-Rs zoZtKGKwi^+WMm9rY;@|+`G8n(H7V}7dY=N-g=~NMG-Dy8B+sC#F%Dvf zex4puO$UiQha+0VUma#)HQ7b{)#f!9gI*C|`+WGOEDgfB`|^sm>mzS9$W05**gUNv zywT{r8!U6NOF?(Z*B}$`{gI9RW)c7oI49huGqa(&nlYxMC=(hzpGKP^+h+B_r7Xnh7m^SW8mRnYL(<2^r) zA0vDF>-IJ|z5ei$YfAlcus@zR=^dUpnvVv|3d+Z<5;2GSBX>JzH5~k|qx!tK99~{- zmk@nII5&zxG4b1>8AQ2q&iZ@!x?QcsKp$T|Dsz(Qh6-W+oPiUd(4pQlFMU&d81HMfAB)VR>zf zoJSOPteoQ+R?z12#Bqi}}rQEl3WD6sia za?R4s5tH9bIy~@C!mXO>@2_c?LeO280=ot+oU~$Z)MbmrIoT)vZG{E+djn7CI@2`# zBKFoNEj1J!3*vS9+OzSlwDx{=)?C!|N~ieV^u`Npzm5qnra;cH&iSU_!9accR+NFE z7`wM-_KpeXq15>bg-8W6^g6d|H2Faiw%>MVYRk{ZMx9&Tqs;Nhw=LhWfwcl+R?{ax zAE1G7rmXGsG19{*tJ^-bl!rGgY*jbQkUUM!+{J|Qa8OxiyM1mn7ws>7F_kKuh9#HUcm|2@OYP^~*QtFK;JvKL@T<1~ zf+uh6aeWho!>^|qZsZC9`+WoG3{HekEP|go*7-sr-^kdJ)0s%!ms5M?a~NK&ou~DS zW<%MNgwqt8U=UE(sHyoHi<;Y>^q3P}{OM0hWP+6u=AOJfFGRS$!8fZVGwAYBcXqcP z$D0IP(TSI2k#Zt=WNBLa)l~d;)6PuIF$?4KcPB13nIrdK=T%pQAYic3rDqR4g~wGi z`OX`Y`{_%b0zc6lT(~;gv%4S&#nUhTGhNQY(QZ0F=l2wR@pVLjmH37H9k=F1suqD5 zll&~NZZ2dE>^d_`IBZqx54;y3IUbE4HgQ|G#gc6XYjLc&5oR@5S$@4pI1Py+Z=6XE z+=9`;LYwUOu3cVpW~{6LdbS6VDc%&ky0js%DLe<++rKfdkeofO?v?K$v0AvI^{o?` zN};!O_plT3MbF3@Bm}Ky;$~@gkBw{5I3b#5`;;>jI!~4JQ`Sr3zkjvgkCXGV1zw(??9FdqhxYrjSF70=|Bp(M} zw};;}#Aid>$hcg)zaL&^5Ib?sJrrsywv8+GdP7G|n9qY2^7}fm&g1F*B-q}gcv>OE z25xc3&cA8*$A4vClqB1OK)G^Qo~eiEz4v|Xd)`p+!$sO^8u7h*sY>{LE2V*T?@7<7 z-C4+#=5@=3C;YkbHUCyg%4umtuYBq?Hyl^h{(^<~~U(`C?mGL1o0Uo`LW!b-0jz=W@T3P0* zu)<H^Ck$5OdHHd3KGLt|XOU}p@=q>kN<743`5YW@UTdyhE5<8AdVZG| z9;#5j>xAeM-S5rmq1i%1(nqO3_BU243p1!^D!(w2eT&$*?K|-lkU6v^*M!`2BI=KS zR3bCLjUKJ+LYH9t&zd*Ki9QfUAMc1yA-SxFQ`&`zSb(2Jgx(9?2?I-mcAEy{bV$9u zlQ}ll9@K9#FVvCqGIdk;Q*SGK{DSXqI{b`=H_ZlViA?2iq0pb7PND+TiZzQQrwCu) zEbP3yE#bcpi=MF52twzGHr?KoO#E6GS+mAn17@bQl`f?!(DLF{a3Sy8_HCE=FU|QP zPwU6CLq75FQj}g;p*$TPPSwSf8f8QA``z)j|AOH%clm1ii~+1W+#LPggy`-RqgeZk z1(54`i|H>%GQ8LEV7T^kFYKp2NUA(RboDjaU8lQa(JQcq*`(M8)uN95_UZN<&dw|tj}K>=A-;0f7dC@*BuDM* zk&Q12hiMwt_SgJM!h*fKc{X}m;=uaE=mY72xToW6_Y&bqUN}&p;GszY=|ier-_44k zuHnw+H#EYD_?5rw9@*Q;_=pQH<|Y!3y|q{LcfubR97yuK6$0O-WFHzgrhrSCh@0Z8 zJiO&6!xv6^)?P}>EwqRx;DB7*wwr{@5d1;OLifHqc$V89`$+SFzWIsWue)g|ap1tj znrQ`k>nl7+CArqNK*e7q7FANZZKwkQN9x`yE0hJ}!8E3i*} zi2{euuDs-#rGah6-jaA5l2@JK(}=s1i_40kyEMKPV!#w9O~A7p2WK*_K3EKg2e-H1 zFFc(EbSEy2tVjf*z6*y(?NAm9?E7AQL)sQ2mK%J9HRIuryTif;rFb|se|=-aZSuVQ z=quuQybn1|F8|J-h{t8a9G%+HIJkd=S3fD(7C&>(9W3N41y9E3*KZImwc!{4dnHDG zIDK)w`!T*G?9N=AJYf_KLboIT{aa3jJ?hhJAG}jwZ~Xce`y+N>%4qJb_UAORMzZ|; zTowhIzb4c+5Z_NBeU*LNO&a(g8n{nUO+sfU`;w5+0Ni$SV_L_1l5>34@+F~-@Ql=6 zNGVPeAK1IN56pMTJjKN`c#Aa(bn133|1#8vnx@}>TpZKE;l`GDjmB(r$KTQi?`;L1X&C)N#NJk z`DgKP!0Q)%%VBS%9tJ(Vlq{6W&-*2CDiGQyy*U?{gYeS8i;Dk#>f?l`nXS7&`Tk@? zSX0R!ytyg-QM5Dh%gC^3r6uQqLgzD)=yx;_*b{BfO*rhAMxx+eWht)5|N6Hs$3wc~ zDNm0e!jsZ5eh_f37}opnAL7L#usJ2SZ7Sw_`tN9A4@Wvir$Xf&Iq3HM{im|GVVFfzRbLb)y^Eh_#48yp zVFR;+*a`h=kem+s=KEhFik9k4a8~EywetIOI}hgq{T{x1)2YRvl+mj1u`3JL32qD! zIh2DQiwW9Uq-U|DJx%P;(*#s9YIfx7D1?Lj+l$^1f8o*b{mF3+r1y8djJw+R66D_? zB&|s1h8+ghPv3}!z|pl!qW|s8$D+F*lWcp-agS59bN3kO$+Hyqp4yZPPr5#PuP6Gm zlAsLr%~TQY+#ou#`E)wzS)J?K9$E+#F}e$uTn1oUwaeT`MhBKW_H^QZg}9|~rTa@& zK2|+ERB9b;3+K5~Bp1yhLHYpOlA2Bts_a3F4@a|6VB2Qx=Ol-%&a3Zq<6|KD_xZ+_ z28IJ~9M==aP*)U4LY8?U%kVspdk&)Uko`4?nnm8*Dc@7luqGoPp)lUgnK^`*)R8= zobw}ncI+I1R&WN~EuPE8A#lcWCw7;jF-uhBNPh($Jf5DI`;Y?W5wE}gK3@#MB|WN_ zEKAVZy_}dYs&M+Y+2ap{ParYPbX4a~5oEVIUmqoXIc7V4ZqnSHggwJQPK;I)erxx> zkwm>1rDH=UjMBDEXuH{)~&Z;1ym%7bxTzhoIj_@zj%nRS;U}hnAZo&3KysjM+dw`s4 zr(o7nG&%`2d3hE7R%F5i8_OH*2Wf;iT2tsuu?Eeq$J85JJaJ)PZt##z7JPJ2`JHPY z11Rv5?(fb@sMCs*<$PQK(jED>g|CuP{8wBJe_jAaF8?i_WUYaJX8U*Y3zOWCgLSdb zlncV#A*y_{FBllUl!yo=^Ls3}W_(%&>}7Vo)w;(QKM%!7-FTXbW$dq8uAB&gZqJXe zMx=9Kq1b=nSz8jC1ieY=$;!d7$LnOpchFE{rD#eZJs&6kY&{?Mj|O8q^Ta$KQbCh# zof&UW1qPYyI<(QF5=-A`^vG>20K*{RE4K_XzTzGd+D7;Aat`#ix z$EF*RQWw^zVAzwyvWFIWxT(TZ=r-Z`DJ9=}1s4lFG0 zsKk7Y1)AdxCd(rHk(xh1%yN9 zoXkE~@fgth5OsA{Gz((3sqSY#7f0q~gWdFzHrSeX5bnM(0PRQWe|mQ&!||@xqpWvg zp}#_+TH@*{-1)`dPJosO$Mb78ue4Dhr^ic5HzXcK{KU>w@1&yrli$T#FOz%vKhajcUcn$zIC! zz>fI*PwZRpNezb7QllHy#K%4|Sp6lqzYI>y3a zMc#aXa4uBaZ0xg@4JE$HsBN?7;*sv)7v^l#!fmg%dmF8q&m0S|Q+d;p|bpM=6k}T@ibxJ_mMsiyt#M5DeR`c^{4%5E5R6#6i!}c&MJ= z&BCjd3AXzVe1n6&;NE|Lu4yF!D*L|}#r@6&zUa&;^&uBrGaXl>cxIxdbcy?7{CKLt2GZM(8j#t{oB^63YqQ_=En;*o8{7j|bL;p&

6w}Q@?oP-7Cx-1uTFSe zjJgtIaicQ19PCO^I1a7 z-WY@uLzgHw(o>&W_5FrWG6oj=)wmLUR=%o5z?8EHyE^s^#=2Hw%}&nMmW6XzIy@@# zFD4O0Dn5OiCOpySby4xP^#*u$mDbWQlZ5=r*i4Dt9jrk~txacCWB>+P1`L!} za|LWwP_bSKolnPsUxo#<#!e0VJ$Gm|cpK@3d^IOq;Gax*BzYGCbphigd4&hBMuGl? z%g^Wj6oP(Xa!b*;D-+|wc9ybEh21G&dPylI%d!voyv88x2t#b6naTvuB~8njkyi2L0k1Fd~E-Z4kZ z;oiFLSbElKqN~-p>ILS5!Jcc|Xb=aQJpSrc>tazs+%rq!bK0$9xb(31M<=&ZgsnnhoFbHObJ zkL{@Wn`w|ta@TyLsi7I*c=^Y5qtz)SDRUqVc_`NtRc5^@Z$tZSXmHV_x)aj zmSqDd`VpWQLGH7D8rp-hQl@CoP$Dn(R2uo4|GjaaNx(Y(mlyTP9z^Zbl)iD229}+g zn++98LRYSb-6@}p@u?tZ{iJRgP5ay46 zzPJCude%(+3E*C*g^QZ%tk#0JDN#1sv5_d}(*)w1`54>~%U80v%;BFB|yUp2^F z9KBO4yJeY*G4JKlItNMb++exol7<($2Y*v}N_^#&Z^brBm1d#O-e6U$RG$J6I}MD)8w(t zM<%J;;isfyQ9U;$i1~3W7R09-JtBThIY0iFjOrvO(aUqeAvp>>U2joOOca6L*TN&O zPv$_)pE|y=h&cGXc9m(^g6RKKyUn88G63=%&kOK6!KQ_qip+$2!JU3D`Z2cwHvDJi zY?xjSee@UY!VL3aSgkeT1@UnRHd$tb6(yomwc)k0oi(uW>FwerX0pfCJGp36kc9&o zGQBmDc_=SkYIGUNeno^~V^3Zl@{YyiXY+M6`4ppByZZFN;qWu z3u;dVMf<_yYbD~(w@JfE`I|+oDN!(C_Q16`Jq{P0gg92ljWJ95yU#j`5!UWK**s`O z?#q|*LltBQN8UrMmD(5!D~E%`e6GmAGTnWLI-OWBHNF(d(Mka|#t+}*`{PM}5vM}Y zmuNiEV9IxCqdxM++#2PMBl%U8BC*?FVo}$K?QF`P6fkP^`Yx>L4W}&BKSh{k!csr? zf$@KdaBEBzKD85lYrUnCzdYeds_Ept-KP)BB8^A&+7d7(_U?14=|SL1+#fWaodQ3? zbqZ|dxPZU!weiu2c-(dRr;RSe;DM{#wF2|);U?YJ|GMrT$1`tTjy>Mv38&^HZH^Y@ z<6D2>5aqf8)W75Ew>FrHJFZ&{(H`U>8>>r9RNdYQ zvgb|O&hr72HUCdw=Bp z*By9w0V@7{k=NUs0XDWi_Xcf9k4bV07t1{1|BUv4{-HqZ>WMBnIOB$4JyFMH%X47= z*PU0K%?WoRTgq3w%@ut8e-B$7l*7P*vJ3j_vf(T5B6sEYD^pig?rx(+Ee=o=P++1T@-$Tc^(!%TPe?l8%_Zak}a zllbnJn$FuVo+^Q(-62Dll_GGG`EBnC(Vaw7Tc}AEmLTfkTeo{P32Ftl{$vwO#U)nX zAuH)DGN1ICK9?mvu{hV%P1k5>KO@RWf1wbBB<;s$AkJs1LMnI@^GK`M2sroefSM+-x=%Uh1auZUd52U7WefD%TYlQ z*m65z?>~cFw7q=beO*}yD0did51!2ho25NVM%Gb4PdCW#aFgV%O)OJe1xX+Lez(k@ z--}S|xH`i-x&$Cd1Ugk}4$My;X1YRnWzm*{EqYWcYQ|8wX1F5pqkH4`V-#E1pRoQ; zRdyn%49BJEeme+@v^2)%%KfH^rXJD8Q1DW9i&YPK1;b9#U{fBRbusYcAK`|*6TCX>R3lhCVJtC%_Rh#h6 zBE;k_ZL9*>Q3o%kdlm3xL(?IuPA0l1Ff^L3Wg?UHXIJ&71>kdTrJR@af7$6J?=~jB zwli0X9|sekdXwIT-v9D5FsQ`9@u-UnV$7=Nv+hiIAMVTg=7Fp75?C|eGPuQ zRC$%cUW==TE_LU&WI*lvL**asb3r7YDyi@@1^J>>PaOBm!_1 zbWqd`T?I>#$s$4VZhRRA3+)TJF`Eia*JdA=5)RIu1g25yoDXbi+bM99*A{i_wie0% zOa!;iLi1tKOnCQv(vL>^1aHf=d`74OJKyczT>v7MGM)~4`4GT$ zG)s&chs;HTw!I`5-m>ZQ`}3zL@J`WyDUFu~{F^53ri}$+|3Q^#3*v`tcBU_W8b||` zA4>f+qR(+oJ!U>DT!x!Rij+TYu0n2|{OIqc*|0g#mhaj!(a)~e(Fc@~p2IIC$9hOU zPf%;>)0KHD2r0D+?GK~DvRL5AJ)ue*Vu+~XDXc+iu$TGASMjikCwjN5DFwE*%yN{= z){s1dMWo#O8jLvCuF~a^4fKL%mh_0uMU!6Sxj^(_q1#$V4Ht^BQF-Zx(aSJUXHna7 zj^v)ywN~0!GAdEW<-H6C$x+?U7~ji9{59j{Ej(}Kvmi)!DJ*Ql2>Zty1BywHLA7)9 z-r~ctAm2CM){veC!;@#+*@$)%zOAT)?@bOyVDu9u<}CQn>;6F(?F^uNwoav=A{^Sv z)$U9#6Y{y<)wpwHLGTO3A12?3|L#9d#={l+k)zn?$df^H{P1%0VcBUi=N63rOS_y4 zp_>woXUr_o%}?M_c{ACo*@jn2{U-Uw*~E)-wuK6-^TN^RxGq;#^18+Q{8pXsAtpkT^_(VY@8r#Z;m*C-(jq z`70Oy)s(IDlY4~CUx%QRM#ZEL8!l~kuEe@GKE^ZeYH)JBU+;=wDsXW6M$g<4(D zv3Bh$RQ9;{HFmKI1N<7U{-)0ZiVTPH8A<{C`gG^3iXYhsAh=6MoVrvprsWMsAz<7j?B%d|d0gL+MQ%YGP-q zF!7oBsIR*=OnN`D)@0_9d=*To#C+Rvy#{XojV{qsuY}9c-@@mfjajC#mmbNb<9TLj zagnMj_F1?p+HGrt{3FUh`PPm<&6 zDr#F>m2rU^*LjkNj{`Rhd~%*E4#TE()J^hA$_Y=-hRnGjG?%-^QzPV*x;rVo&%*D@}Kc(Y~Z(SCb zFBhXw_ysSCB?_Jk*B(maNC*39@g<#iWZ$?kBkl%I23*^i6Cy2?0bUkvjO$j)aQ$ba z@<+N9{OnC%o+lOwDlP7xMqIVwSPG{_KuZ=%oPBY-&>;eg8(tXy6ElRno3~le3mnEP z?!_YdQG}y=OxL+dB@IJwe%1A;3dc>Jy_2Sv+So#`cSWis8yQSuR?_xTaKEh%ES!tP zyw}xHqoahwaO?=HskJFmQ6pu;Pc;-+v$oOYF(Kz&Z!6zI6yeKgevdRGbE)Zu6rr?a zOkjU6GcI5YX4Kh}v8medlKp9M1YH5@&il~IE9H_sk)G6lEg{&&!JdEZX%0Tw)5Ll| zE)?gFX&m!rivkl@N~1VW25#XN-Zx?5fMbg*w+!=1{{HYKoqUsA9BDM%cz1U!)Nq$D zhDX`MEj|;DLyu{=xybyumkQzWKeJefIG6`4{kPgjiH>gPs$0}#L~<@V2Sxoe3Ndm& zy+Au>4l-6(4s214qVOuhdwoLzD^{L?{OyTZ1#UC z%|qr-_Li314`JB7<%!bW|K}I2k8~z4)xh89P3+5^q>p%RgUfdE{Kx($_4Y9#Tcwxx z_~6NGoboa^+IF)Bq&`OVI+EYbhNN}Pa3c@dHO6W#B~kI<9joa*bd~UlhvmqlPO=wO zGP)AblY`-XxCV&-%`q?D=oHIf?G0VCC+$-Z2%zk@Pk_x!*SoUVR)2AIz(S z#t)~)#Fa~-f_GSWop3rv_TKonNcy#Y{YpANz)%eQasvUeB}I@KVM5FLF9Ta&+@3l1 zD+{-ji%8rVB>fl{-jnF7_pi<3 zpD&i-{?07_62oFN`C}2cjynqku33?dYz+K95&lD@Djfw*Z)F`$3PToK{aAjn5Xe#B zw2#yCgBHad-Cy;Ik2CoRW3px{hVW{-50d_jwD}ivW!oaLdU!#nEYTScD>7EHDG^@$ zk&mYRf?*Kyu<&fGzBTTr-D14FgY4DC|G764UE5ck)pq+gaonyl|NQ;0F!a&*raez` zat{vYpJ1{hdw>vn6Cv^(SpH(8O<&DItF}8O;XP!3kRN&e&}j-}h1=C{Yl=s!qU*0F z1~QQUP^?hhHU}K!5;#8_9|B4h2JSKzL`OVe;^xa+ig#1Bl@}~2z(3P;J7X>b=&B@k zRjJjYi&s!KTUIrW_h@g5yPXQ}g={Y7PWytrp8SEnwS0Uf*yZ>|rUG|-kE#DD~0_Afn`i|4uL3^D(%MA|sJz z3wz8$n^gWOMUL%T(4#dn>i2NpV=A%dGtML7AAUNyzR;P z33rl@eL2qiOD++jY}jOOa5#dUV#Z`YLccj6+M|y+Ez%p1MMkFKLkKFHXnLhaOth?n%UXG|F^o%LbW!-gv}8%f%-nVn9X!xPuuG-G^1Z>Y_g1lpY%KN*uC8t z^L5G1=XC|j*VePJ&=;Z;8QgM|^f}O&)LXK0 zG9c~NS=c}&+yemtO-;FC{B$wbu%Aa{qW%Bd)`ci~Kgknbvm z^kX}@?5w|@1N2v4Q?wRrVPb|&>83*>WVTruG+hXQL#q378caz4kf-JktA1b9p-*kU znvskg_3~GYFQns_zxq6z=bdot&d}SovdOqV%J;QgO9pCwq~=6b)nKj8MWF)YDuO8F zO4GQL51uaDTwWd^y$_`xo7lhPKqqfnwL4=n9L@hwb9bojUy#|63Hz=0yRBr$L-$?1>XAXh)tdpoG{RwOVlr)i`aThuwSKK{|6>c+rtNBY z+@0{~G5t3W)|&tiPt1 z*`j6X0r>lK-j^E&cx#nXE3xDP6#B7+(2Yf4d+0e=|JQ7&lT`0`Q{fHk&*rA|nx?|; zsP(h=59Gtzp)Dd|t1;k{&{y2BJ_nC%v?=Z? zJ#gxv7_aOvf6zR2sl2|c0D=mRafhrU9G%0G-urn8yV35NBk)=4Fef`_!SC{p`*ePe0I*F&BU6qg{GP(Y`?Fvv@;}+r)u|A4DeA1xG52bKDAUkus9)8s*92TzC@&x;QYIv=(K$o9a`Jt4A3Wb4`&`@EXt(a4_t~o z3~|3FGq%pU0WUp+;O==x(Mg(uRsOcx~iG?~8$z13CZk z_yBkQy2K+pQjlM&S6(*r(knb5L&-A#03A&9#+pS5t#z|>z)7!uDF5RO{iA=%bw{3_&}CoUF) z_v*t$L&&H8W!y~1?p+9tD`PZ?G<}t4>jjyv+rIU8 zmBWy1l(i!1B~GS)a4_u%$?-68_;K7Pyg;XW?B{h#;dU$Q(5t~Jc>n6@(T}c)xSvnn z$7qC#`WbS;j->ab(Q#^J+m1>wklZ?U{N8bN`1|(SyS5B;yfn%DYd98y-Fv?*`UYV6 zA3iPB!{zYx)zQFjOw~|IS$-tlRt5XD+FKPFYk~jtaoc@;spS8dQ+n@MK5Ex%h&UW5 zf(1s-4mpyG$n*Vg&~r}!25f(9ZxNJ(wXQ`E{bkFciifdkq=E{kBWoqj6jQNk;#;4+ zPd1+JVn3g#T!fit##k~{34eC|PufatIvl1-*Nlib!_BAS^@BC!eA=FPe`NxX|sY*WRj!v)!O-Ew3vJJCH`Z0 z0+^^aLP2AN*HJoqjYrCi=*Zr+xOqt0?UuR5vlE&X6&V{ieuhWqStWs3yw zK;~LiM;hGx*rOcZ~9m!7-B;ADhVM z@7rokf5k#DVjG>h&{+<9r9~}cZPSn~Vc(Xq$N6|>6>rr%$%8SB-#lsFO9R0W*)G4w z#CKFxK@o`{zPI~|k9O+jL&51|2LI3J{8Rw@7jDE~@88}WmRtgCJV#^#iO1u7QpvZQ z6+Nc{M}FQ66sR7_&m$Z=0xnBbAE#OYdc5=94^a7S~FJ4(Y3# zEUxFcQwT-Pmet>L$mgfVR&M591r)p-;tVrQ#Mi?=`?!5*2GlP84 zU|J>I*6WChm8Mg{(HVHr=1k`n9&=PnU*2f0oQelF{W;p%a4@*jO@pDFFYiF-`Pwdj~sHJsi!iO`$Ch+|EFZ%k&Z0*<+M+vgydrDXN>q2 zq9{23(Z%GVQV8a~y`ZU??*uWQEvlFPIid5#pead>B<#A^nHybMfj7=hoQ$X|!eG@L z_41ZN9Qo3Dod9|MvpGakNU!U*m7J=bm>2fd z@V~!%<0yPzS0bx5o`uDO&4;d$e&UK?Zm&mKH1^)J-gb{a1MLDka{Tt@;-g3A*VjA6 z;`OBWT%*hGu=mhR_i6IIBoXHCmZ@Bh(k&5Rcanar)!L+~U#qGxQ>SBcw?;LJbZ3Y^ zFQq{IwI}&&f9Al^?xRAsYYMSYP<(OQTp5O!GhaUPi}Z~hPdc`TQQ^pQHf#OqRG1vq z*mqsX1$?H0*Qw8td;#a+x{^GSdvn~l^0cG~grBe5QgM%X;!?`{rpf&}!8}deOWPh> z`&V5I?54wiU5B>S2bMwGma%x6cN~fpH}o_o(UB3;8M}u9~l4T2H>g8XSf9JzX>GE=?IU3xSbvk`;%}H$c*u}0W z=8Fp&GD+jIdC@g?K=|De4+H*Ccf3nwOHW6_P&&6D z|9LXK+*;x0elQ3;zHiv{VRa;AcbQK-uAoBQoMm0yupY=CY*+cfhy?TVO}tgE@#u4D z-{a4O%W~tVDMwU#CN}MgR)1LHk8iXGykDDS;T84E^jA+PIQEj~^tHBXd?S z3X2xmSnl~o^-)3sIDF$@X*0AWJ!r*j?lcJy;E@^~ZHY$tfZy!**^5x9#W|#SLl(NR z7wDDO`@_{_E~@k-@m@Kt@)#zs1D>QuS>2AuM3xHfU8w=+$6qC@z?^~AlIJb^uEkiZ&eU!hUm%Q(uy@{0JPXL2of1cDP(1G((xt=KD4vft!+;?@-!=Uo>x~+|LFz$5e z`m(ni_#fzbXOerH_WJdxvpEsnUHr10`3OgBOz(NeU?9lt8OX8C&jMGbL$h+ZUIe#X zFT=Mj1H5YgiJz1$gnO*5=VB!_^D!ux@W~^chse3S&Ak5Pnd%J8DyiNYY?BIR zYuEUlsc|QMku0xGd107)UpkfRsV^8y&N76V9bm8ihCOMgLowbodP#q%6y1bY76nPq zf>TOfINT;2T5R9EvHD_;%9%b_WhlAW{GoV2>uVg|lvVKJ7abOIcJ(>ZPo-b3)=7T^fO=I1|!6=NC?lU{p;tZXyJrBNhCBFF` zBLsWngO)e@cD)d#;n#{+dzq-opb-*1_Zsp+$vC6PtTPF`**ddQ7u>LG#9^UICI>!t zTYNa%M>xv|C~_NKo`vsmo0@1MX{cx>lgoO|8qTD5ZU0Mr;5#yxu6!ZAqt^Vz<`Ew| zDDXW+RBAa8aQ*ix#W$tEICS!faY+!w-2CjX6%>hXCLfzebFxt^UnTUpNhIzqS3c@y zNP$x`&FvLaJ}{gzAZGs18s6Guc~IQa03Ns<+n?+UN|fpd)677uVB7p#gxs54Mb?Dd zI^^KhweOl5{=`B@ckHB#&jZgiqK3zI~2f5arDU{mCqXJ4MCBHQ)B<9DaYy(2F8!5_jWaFpuZr(lr| zCf?Q)2TuhcOK%l>YFKO4ODA~mbU4P^lDKtWA!rK^HuX>B=ZyhqCn`;GaR>x^+w&Q#Ofo_0#Fv12lH@ zyd`toU4Pedh!ZMe3Ajr9(jHQamBgcThS&6@T1FN!?2Gbk3JTHZ=$N17+kEIbe#G>b zWG=)DHLj=zguuzhwzn@VXmI?T!%pi0!Y%2`x}x*I50z>|>NW&r;M-$Oa_Z@KOdy}6!eO)o}<#_4vj({RF3JS?c?fabUf=Op=)05MKz~|uKE&2K0 z*m=uQVsU2~nrNSwF}_K1e0@KPckRf6b;G}^?6pH6GCXa6Gu{To?Z>V@~!8 z!{T|vwkUMvyq2oC54dhO;GZS;mRIBN#GbKdqVf70Z@S3&wn=!Wrs!lG_V0Xd;Z8j8 zFS_N9oVQAcxkh8Xk@{4OF^u`1wqT6QUn27Mo%X?dUFUjF`DdWUiiz5jkyQM7^VoqV zH%qKY<*!%r41$t*c7mt(Cijasrv?MlLH6D2t1QFT@Tsos_v*b)7++(#d($P7C*0L} z{vt;rS{A!a6dbh0vpYXZBpRiIUV4S5TC^P`JagxE8VmDBypGTAWH?Z9lj(S>l5 zg7>r*WZ})TJgYjD1EGG9v8S8#Gb6U1cA6OogORD)3Ks8RTnykj|LA}r?!^&GgM~hf z8T*7f3l+j;^G6~N4&}f)^P$u&y>!@e@!%72!k>9?TTSmIlO5s1y|}~Zjf1OyLm%0Q z5-)?mjjBDH%HX=~1H}!gcJT9rIky7o1)a;1`0rWib1%vx%%2R$N!^_Fx zy}M-W2_LA`PCmQp|ex_<4NdJQcK*hIL%zYP(Q0M8gGF-G#T?i zBof}yg>Q7!d!YA5EA8B0g;?$R(l&t1qg&ir!#eJfdE&A-m)wU?7)B z64*iZH~cI9LD$m3>T&eh{UrA|Bre2QatsG62iDw@1u76SjE0Wrc}%$0J1p;Ah5mhls(V&a?B_^8?UmujBMe&nHKnhu_iQ^obqCBAr)Xcv_k)Ju z9jo^(hJo~It!q8MqLICsk7Xg)7MIqFI0oq?;sC)MwLGMPMgPTD3pew@y3+Si@_H9+ z+d9e3buk0KwkAkL6Tc|u`+cEqdJG7ExVz(`M+_L56|FI&dVQ&DTuzD?0S z0f)@BS(!V(->1(_)tc|sfyB4nr92nYuv<}at;k7h*VJ7@b_MLIQU$wIE$i|x{ zSMQLXkZ_5=;%-0eviDebNV^<&#m=k7Oy=PuzMF%hcSEs)mnlG5D;;Z%iiSMricn4J z75nkAGThv=^g?ET0*LF|{is%_!kgYZMK=z|<12b}>a)^p+}-5Q({zaf%}P~2`=>I& z;Aey~RusXrSroqYIvdI=Ee@3w=i||`eU$Z2sW^82vz7{(E7aGv3_c*e)(($ABinQ; zGP$Q6KfV+}xG8}uRiXKC-}p?zsrRJM7fRnJYh#Ju?iO^tr(~h?%PhMZazB;K^qcC) zDj;4jN$N*OGOxN9aaWr^4}%{EG*WZ}&``h2_ST_59M9&x^J5PKWOT%TJD&W8rE3LORj1N$Hu=qzT5~f25?V}v?IO5 zhDAXkW0KFc&sBslGsmyhC_1cUx9@A9ze+kjxw^L$+L;^3TGM|i^@Su}5N z$$_l-ZBG+F5e~%8j6LwV06$n%+$*z-0X^m6o6{_W$8H_4g37l4T}sA(&GhKYkdMDAz3mo)ilL2f4||Y%862oAw=Gl5!33I> zrlKSRUv(uLaQKwMa}{RI*XxLX>a5$k>AVnJbMV{50TUYX)I=Ia#S`vPs3K3P;(5@g zPjjDNPQwe^c9w0~LiSd_BQraI^o29luDYj?2cg+~Y^wsi(OTe6<@M{SxIup5*=N@> zm@k`ZR1v6vt#AG=t9d42@jr%k#Aqh^UUy`}<{)(26_RnP!yZ&=Z`|}nqA=Ii-f!*Q zG+ffJ3NJ}aBzt<-suGb5&`Rg>4Cl|pbgFvI$)GH}XzJDKwv>RDRGAdb`y{W`oj0pT zyj6-vk8Y6g$i|}v0vsj}+>j&PB#Kh93wPJf>q-;f!GjszzWV)iWQ(Twt_(W^(@ag3bjxZ2v4mkEEDAdih=N(=DZ%u%}oWN+U%(W{?AI2S)f4*EXLqGR)c z(HG;?BK+H5F~^--2Elp~&(@7sLG1WdF+al?Y;bg(I6(Nowd#^KeX5nPZ!(-&cd-ik z`~@Nt$v#H!?#nR$dwH1tPFq&;58;LUugz_4#82Ywk+K4mHI=?vhl&qt}7BQ}ude zP+JP_xFO;C>lEQ4IB%}|QaWR`#;zgg`QV#ls768&i{VV(ZH?AklIrIE;Ydg>&jbCRGS2V z#5ytF1>P)3ONxkh`{D=HWfj+pnI-U&KF7iSPG>l!`_{w0g$9Xo0d0MUQsB7#{mxUR zXJAZY?WyNa3?b}@RQ>AjbV#qYRv9Gcp|}Qx6NdFv$oYC_`F?UXG>*^F*nG}lMl#>h z4kvGjw~|vkG~oa!qDYU7N(7-R>;HzR(P1(~G&@Nz9xlctG#m}~!Cgj?-#!Z&;Bac0 zrv%|e_OVH>wiu+s`{}+-t?3j{dQw)pGrbV@^YZ1UKq#g>J}a^DOEzlv*hxMJ||@DF+%6zv;N5fSgj=>TSL!; ztQh`JDcAFH$a7=dv`H@dtMPWXJTHPJU;eoJrFp<1xVkWe>@V+L?|PVH?gU8}W3)po z3bEk8(n>}S;V1RR_g0dg(Wx%w&bkodQ=3pXY8>@L_2Eptrn)S&SSPvb*x534c@Uz% zcWXZ3rZ8mAkUa|Xg~Yk+&Mdfb;NM0gxoRA;NC^2FMfSQ=ys1J7gmiAued50dR4 zO`XU-@7=nr0#nxWn5+5e_i5I2;K?r2x7rws$Xh+CU;@%FFQ{89$3kJLgX$HB$rjrT=3919{ z{shBLo(t8Lr}J=Y$vHzlmK@9}<}{s2&4k7D(T%Rm>F{~vj_~>a3gM!58BefUE^LwY zeRAEbV^>{;Lc0DF-9$*L`UmIXsQlC$4b zCIKzrqD(6z9>k4e{bjUG;Jf(jg}|*8_;)e**r40 zlnbjA{?w)wGcZ<%H?T;c2>rYm8`F5`SmIUBf8%r!CW|O>rduV!Ip)!5leIx`lKLj5 zDKsBt_3rC`BG36!pX4hqR14tn<+5fra$k}WB-DLF9dJ8T$U@9@i#42~M3$ zY?16_Z)iOvJ9ae(zvMW*=h#q&f;)or8eLp)a7o)dswffWem!}vPx?^1MNTXyDVznP zfYAq>5#^|%5Xvd$P>o(`ZH#-m#8>v>zT}3LQfxfeHa*Ih1177D+#ObBz&5qb<)1&% zA$Usyh3Y|vquSeDPK%PBuYhw;bbSG)w9wy51*O2p{~Wtxk~2XiIxPD&ZvcK){QYyi zRw|aJ80G|3kom&kic|`Dt{Bme&&6=$AXA)k8gmcv;0XrGid`o8TXDf3s`o3PmpQ5B zFxj884$ow7|6zloYWlWyw~Ih#rF!s&S0(%&Us|(==8me#oO}N4qhdl+CVlh`;X~4$ zRp-+SV2ZcO3{Qn_FF`G z&IED(<{k_0ER>V5lsxg0wg%CLBEXVL38m$Xmi{Ur&Fbmb{yj`%3ej_i5%)-?>Z;A+4Z8C(p(Jha-IeGZKA;~ zH>012Rb3!z$ARKd4dP1(UD?7pQw3?Z!4bPQormXM0u}t%NPcOoeqHdCK?<019&*RHQHLr)cg?YIcd}>(z++zduI?1ZK`fCh)Qdw8IEt}-PPL~(2 zdL9Qqe^wRO>rqg5dZhHx-Z-qL&BtYwra@am?99!D*9}zT+=V zVD`CgmR>Xi6*zu9`p=k-`L4QO{-oLxF99F$_{AgTPxqI$uLYQRE#9&gL1M{H{%?FWtSCj=UpQ6;9EZi5wGF6PSHC}2(vHFSvM)p8bXfp2B#FM!T0K& zNlOoGaY{t%=%)ArWSv-f{+sw~?=XkbV^^QU!C!~qbA23!oUY#sp9p8{&$ExK$UP;J zo9(a7Vmic{JulL^p9F9EPhPq@P=w}XyyivW1vn~Unw-<`3fy{0(NSbycx>~G3fox* zQsgf>ScKECDR$P)AfC*D+AnSyf-)$(ox4?QMvQ}jtgn*SmB-O`E|m?vJ$ zCsPFpVkfy>kCy=mCaE`)x##cYPlk`Lrem_ori=?3<@yjWeGI(UZdV8KB1;@V~ zQZVhxz&>Ni)-w4j=#o?(J^BB+?TxsJ2i$I0+ShJ>>RS+gVRu{fkjR3~VV0}3o@Bx! zdZGV~saRmPRB)bteHI3P1Y9!^h=qVvhLs$$>7aYn-Nc#nmw5C&(-e9GVY1ML@l`bj zI{Fgo=JGS(uNnKmsmODvn0NC1&{QxuupC;yvgiY=)-k=;$}z#|oQQtKr}l71OKVXg z@85j}hrZyI3#3Pu{?R`7a4sAZEjX(4-5-9p9i}=PI^giVTOmh$GGWMLkrqqKhUyna z(G_Ce=&9lQ`QQFvC^cy_X3|N4hLqj+_eUK8mAC7*&mJNBzZ2dgBHv@-Yxj@EJmUZq z<&b;8yb=mR@;~{HW+vjYy)dgxM;cyTFJM-ZSO_=$EsEdC5&!VYsa<@b6g<#I9n>ZB zfC@38{ypW{z_DRRY(`xHq+JN))sjvIB|A!9X<9mbG|iFNWK@WA7oMCPf5O21nHd*a zwve9GSmqnx3txtaOC=AhkG`2ENA>4Rd` zn06q|e4WF2q;}t`c-cVatm;yQQ5AV$8PM-{ZN3y^MaK@FAl&y4wgyjRSQEfoNMp?F zUoH$TY`<&tp#nFU8vVDQqZAW^Rfew_67Pg>`mbYTUh^+<`JRqI0Vs-}DLbn`=CZci zraO}fZ`J;dn&y=x{C!(l3_vyGfHn+<6aI@P52^Zz8uaSKhnlFMf!^6>eX@FshDSNVmhrw z{C6#x<`yy(s1qERlk*LLbH=MBveoR#d87JM*g_bn2EE`cO-=>NY{#PZqopLjyD<`;%e&s_7q8EOvC@q`Pa9|u%pxN7_Dd`ji!rraSYry*OgYaZe%Pt2G-<@8J)OL zh+^XFPScXBu>52%b6gkk-kmlnR&&aMj(eO+sqBUDEN;2iCNmYz*6p>N2_l~KlBe=N zYtw+k`jdURv>oI>9}&8ERv-5YG<;Yy83AW(P7hVCF98d|IH&M&lA~H?VavQ<4RHf| zy*6H}h8P~{?mG=sY;!wob|f_y*}pQC&(;tQ_QCZg%0ZPtmtOuRvYO2COu2LbKt9*>Kb>KNVu4E`_;r89fh(|kAJM9;2?9~+k)YE zK=!RSDr%i@OO{o$f(-GkKJ`}?iY&wZG`3AmUKH5U*kzb9Q2^Y^A5&!5bMbb-Q>6pu zWq88m>AOIIJc#JH9551E1Z_QkGAQGK_ZMf}F5C`4=M^Ov_rv)(RRU5^8%S!U($=HD5({EeoY49jm4d4b2M3vC8>tg z9IXV@DGB^2^*9TCOxFa-#04X#=;uv_gcCC(we)Z2rW}+#*I093E*e+0IhlX@Lh>5F zC+2qsCV|Llk<=r^7usj3)WN0}3)er*YP8BLWBZx-h&2vbaM?FPH#aB^ey--0C>tZ5 ztb={0)S}8@wGVfWsWs_K@r13tmy(9}uCu05=c{2&hJ*DoZxy_3FFUD{kcdAyOuXVG z+%bR34+;;^U{ak!^<_*kj2&aYw;=Bdze?h@3|_{g&`CAC*~5do@)^*^x%G6fE*<`ut4T?xMq*jf<_&sVvoV?T`kF)8sW_S(m%eS5^hN`3 zY&+KGhQEh1CST8y+<^3z61T%Xcwf@g@M=;fIyB8pp9>^B9WF~zyIyNR~!y~X^Oo?-d`r2rbnOZ5D)8}>6IN%UEq*l_NY2}J<|@k-mGGc zLse_$%Y*en_-3eO>kAGV1|Hdl${!N&O~lF(WIYW5H0j6teaX4RQ{~sc@LW(XA7Yvy z=j4Nt5_5eKQ5d%zc#~I&0*#M~o=>(rfyr7v(JL2%vCDj6LjR=$2tFB2QS>i|#JQi$ z2cq+!{U>{g;_)-j1gLPZgbi86KX9V%^UKbGYr`xn%8!=z&5o`}<2ueV78n z*VEwlej3iUN7~Ki#p20o-9OIWVPGTU9xBQgfJ4VwCLbM*LHC|(54PtRgTp7`X%~YE zFfRR})UHxM`a2qif3_syfUEwnwXP2;Ts*QT~z zUq-_TH!a=&ID_Ca#ge~EJ{1~$Wo$(T3*aC3K5Y(*G`O`|cdr;b$y0o~)LS&@1{Wzq z9C!1g(b4Ed2S0NRa+cVr|0RA84rUhblV>wP_ju%9*R(|V>mJMf_`YbsXDeG%>M2m&wh>ik}|Ky3Xc!Jelb2T~u_9^a&w4c*s6>cf-+;Aykb=Va$F zOpMZ0l?gtDj@Z_~tH6{jo>L zy?$2gz>HA^wl!`&qM9fJJgpj*mD3bxn~!-RLAX-42Q+he8F|>9cCvBfzj)NCJK#be z2!zV+7-xgc+XIP0}_$_&bAq$u;_IC z$(g?asI^*GUi}OO4eh_JPEgT7rJaiw{;O@#Gs;AF{u_BcJl~nm@~1#SfY@`XuKjrT z3Rg$lSS}oR%v5XE84p?eKeagy(}2}#`$yr-VzA18d}ZBuDAWnLjuT)T+Rhbt&{bXV z(tnZVEUqN4qp-$!%}fCdC1HGPCdn%aDXlsgLp*lwIr}ahazPd4mWIu7)o^OMMm2F; zC9H>!TaO*R#DAns30r1Wcl_(1bztJ z8pAww672O)%ygGU1K*e!OGc;~;h;z9*MvdpYpJT+^2<| z?*u7MIzznc{=B1gg{Y$#aIa)B4?m>UIhJqF#n@E#k#_BTFG|^VvuI$xfhvB<`UKq95bTAkvlSd)Tw}Y?}x&DlrEHkdX++}SW;_;#{ipn$t8%bCZL`e4hXOR`4Yk$k zc7~t?!y5fm@~jBk!Zgppz^CIaAxVU9===DGW@!)kJ$q7=*Dm|x`D}9;&1;zuq$8t$ z`D!E_<+hT)d&LuHru!3O4V`eGztcrpNj!Xb!Sy`R+7fG35x zv3sNMX=LOo-^wrxgOVa1zv8cDuy@$XMf4yYp52P+ZKlLRB9Bn-WLpOC(-gOxk>0|X z!0&&b<|^RRJYDJ|$s_GLvcKS201a1rZt1w%NrUadjCd70;&D2t^H|uP?8WWJf<`Ba z$C%%2Wa~Q8<5hU{&z}4qs@xm(xbM(mhw5d%San0sq?XOY34Ta^n(p|%|@4f+~MQAnny-9}?@8b_W7pO$vJv(ia*Pq89 zi}Ob{|HebSh3~DnquC&AY?vV35`FnjISy-!|`|ICjX}qLGsrdCW4#W=AO8+Gk zAj{}~CsbWaQNpmPLmy&6^nBriGc_3yxAerOvV`(7DE`wdOnBR+8Y(8CtXwzIZrWL<7+|I+*?1+JH=X7=Lhw>3ZTj?|O@-S^KjNxq$R19!wnX!F zHCCm{MbxRJz-YbnK+Hle%=2X|E-aI{|wrUe=h?zDP zE#DLmw%2dmZKbB*YDV&a_)QxA+8yrs+PwfuA53d$9-+V=jTNODg*aqhN~t*ePZx`u z=B=$C&|rK?cW$sI1=y~B@pbx`j=wC`&la7D$B$EaTx+88(Nl|JrqrBHxR(8oP3Pjl zRQ{B;)z32^rhVXRvoPsvAGkhOMsog(z2{^;cayxvAGN5fE-CnY!n$mWpF8O#Nb}RS zh2fB{?cmPIWE5f;y?vs?22b~=j%wHH!qWF?9?qTqcw@W&>XuVfboQ3|=9(6Zbb+MB zOgn3wH;Qw1YUhJ=$X?66IS*IH+k7co>6kFnwEx4kObiKs**ADD2cyLox%wK`W0#Tk z6Ze;%_{ws4{k_3JEcsTsD}PfK%1@LlDad9c%jsp)o=AU;Gh?4@OtJ>WS3^?GqF(qX zoL}U}XX4G4vDG~BAszWqd#RB0sWzHO988vAAeAcSScMH}OTY+)|O-u@H}MVk>-@ z?=YZOIs&!G&;8Sl59-b>g6x+cv`j6DfBezt(GSb|xPf2#v$0|v$d*eb7!fawq2NUm z|JF?S6Bo|B@em!IhTL+!Zls}fwUU!=X(fDBxh1!FwG2c@vV55P3Q_h^{=K_D=wuJX zbDSxu9B%WE?++mR6V}3qgBx1nP$;bY)5`iBPhH zi)U4*bDW2sI~q>ue;#^N z0rwwt#wApe*KxhF_MLeeTJ?9W9zI+QVYRv~qGaEw&ct}Kh&0d-e z5rgZmL~6JiB;toP2cu^63@}4~Q!_Ltp-RnMvRhC#?z7w6e%DJMH}9{$xFnx|6Cayi zUfvLju~Ld!HCFoANfU2moAp6X<61%I+6)}Ip~I8poP{R@tUdo7w#Sd-EAn**a`Ac7 zJLy;VlCbw??eBfpiovXYBt?}{4D#+h-=reaF(r92=il*6>|4sR6yK5x-#;o&$kAz_ z;QOERiXa1};yX@WKA(xIg5xOlg$h>FYxta4VxV(*r>?aYnJcd?ZRQOw!yEFJB}cfU z!Cph2B74jRRi}18ZEY+=xyJ|i?+8?*6W1DvDv1~v+t_G)l!^3gI9ntfKV&0^?nmvF zj8e=g^6j3i&Vox*|CNm1%>!;pNs;f%@lf}Y{+82%^g*aC7MH4uP|KT7%tf7!>v;8* zv{z_wsWU)`c{v-5RsS2}5-7kUpWePbkeEj1#0PH3mNr8{hgYXc)0te)Wk4*y&cd?L#%HUFAH}yM_3t zx^^1WfZaM_U6K-#kug@PX@qf%u9%i4h6rvcWRGcD1hh2VPWCnS#VUoWRUeF>5Hl+ z$G#0G9!-8L^H2Om@MerH>NMfpExk}|u1g5Vle<}5H#wQ2%zh1*4y_^(@l)aAEzX6Y zns3&hS5>0a_Dx&v@KoYg%P{SL{CwE@EZ6cD;SXwsbo5kRDMB-a9wup?Jd{(GR~UMq z4{jW9ZyLGg!UX$X?xL(}{5HQ*)?QMLdn?uLe{Q2e^i}o;ED5>r=;=)D^mWps^kU6l z^}YgGel;zBChxPswf*;J(#gD>)-V;B3Go@z z4PQ4`qKSV$+g^<_^h?&^Vq%VlL@~9?Y`dZWhP-9At*7GIPf0FWbs4B$dHjT7vM=5m zGPT@mrVW%D&H>}_Bs^C3VkoemwhhZ?qegkp-RGpABp$HT9DOSgorUt-T#D21 zmFiZFlluN3aZvlw!P8zi8d1K;+7pXe$#>2Pm|3Af;Pw`c8bdVwB^Ot6`wTYME~NWJ z=Ucc$lboa%HnPrB_nruAPX?aN6HU(douqMjH`a*r{tyIlZ;mE9g zdG=~H;SS6FO}7(?1Ho$o*=3*l5OIH|i{B0kDW z=K^P^&_cX($7p(@A{EU{63CY6f-B4y!8JU_T*5QDUYGhj)-b?v6ELU^}FW0kCKJn)R4(>~Gf zkIPzhCSB*U;jC<&S8aU&gonIi`LvdBd(}2FjgAt|Gqba!AXY=|B-i29WWOLE*?j%I zS`Km?nD;E&nU721AAEF4E@nDREi*Q=8hkIDJ!nF>ZJ&CsB?<)QW0ARd?W5Q#*kN%x z;v0VjSY1x7|L;;Bs%+*>ijShB;gCv)n>^uRdS2MqHd2V}S3P@D6^g;&fbQp`ce5d> zBH;;pG2v&8nKTGcDnOL;Z%X0244l3*cPXKQio!3Z)OPUGfk*G}N5!iQsJ?ex$}>8S zd>2eE&-2qz%zW47Jul;syHwRK^%e~!$Nh9fVm!gW!}zh?Kr$qLR^G@J9Ev8=*PjS! z`@=K!H-mP2iZD)p$f0z+6d&AauTeUj1!g^Jo1Mw$=gl^4%?pGtk#pOXR(>rV2dbSn z4m~~r{zsrR=x;D|NuBC8`x1x5Wz)qyl7SDre?8EvFGZ(P*0N=L!h6=O{wZD)11W!K zVV72m0-k>V@iKz~Nj#Ti`@G}fMZ{*RIpHx(?0J)M*MQs;M3`q-d00WeYhSLcwKU4g z+!=ZBlJI)Xi}J?oskmvEjm5VS!dq$fTAw{b#V+nXZAsEoyvJ((YO}OE1bs-~uXf%6 z)qSESov$VkE~L~tJ&Sw{XB8vbz+~(h+U9;v>paT&`TzY+_OyJ%ld-o64{lysQ(J=c z!K$Q-lruh6AWcD0Q@*7N-yW6FwG5}hpLv_;FRMtu?Q!o3+W`}NY+#|Gpqz#|i-IlR zcnHs6lyZcL?u<58bL#PxS+MC0b7Ng{IJ9o(I?I+Cig)Ipl|C2>hqbZ(E+O`LP;Qf0 z#PMGVyxhjz?@73Kjx4M8&{`Q-xaUQ(`+szZV`@`taUgkvR}894dI9o&$n_kkDZ*d- z!<1%>44{VnL2c*DNH``fe49!*FqSVy*WYwaB>Sflc%`)ujTBTvMsoa(KJ`^<5j_LoV zXOH_a-u&XCA2=H|)JM%{qTP!yT{Q{9<91%>_ia-pcDHb}Zu-T*N2@v{FBDQi*zy6l z+|PJ0i06)}Ao+eN_i8z%4F}Nt{qI^%r*fQprm!ZuhTO|4guD|*lgQqVU45{Icq}=c zZRc%Du*|5{UQdgL#oW(bg-*M`*70MLk=Y^WW%5?QXEW&ymsyX8?V*L zuI7vg&ph6l*gOrRS^o*>^Od=s!S=oCZd!OSNhhxu5B0HjIQzS)J@fIo>m8din zKG*y6`Th&v*AJY-`H9=@^}4R-^?2Ox8z;-j zC0h`f%AP-Ek`xbLD0yE!ESMmA^5^9e!)SQc+%L`Su*K2IPp1d|#9+jh_dSsLC$IN^Q2f}pym9U(W~}8aQ|ra)ZfE&*l8E~H;$|;{#58#Wbm1Sxvx^L+*86& zwxpXJaAM%4=iceB*Q9{FtE!uUb2QM7Z;lzPNW{6sRc{aeu*0LVbWfYJkuY^A^_J!9 zb0E-D7-y(ci1!}MnDlvOq0IK`nN3gQK=i`Q-=#c%SSK-{Td|al(URXapNuDA-O+na zf+QcoNH@3RH}ygz>QN1WTncWR+RqRubMrBXxvTFurKh4DNA8sRoF~@a z$eUHL@q;tzypvAmmvG6uyfn=5U(WjNmp@2fu}Gs`Dv8AxL$-J;$l7@VKi|kFg)ka) zy|78}Qg_GxYKPI(Iv>zb^5yF1Cp~+q#DN>%-*zkN*;GA%(g+A)R<-&Tis-+NFbYPeS;z4<_JsQ7-y|^#P=*UR~W4 z?+Xno#~XULWWb6<(bwn&}ByZnki zpxKlKcfYZ5%efXp_+W14WlkUD2$KBCb~73ORBK0%@F*J-D`VUHvYz|m57qirS+Xg( z|NO1vBEGp$I}9nhl?ybCo!B?6 z-RyyDx(!*6&GsyMlkK$@Pq1 zG#8zBpDrI927rV?$oCcsnU60$8e!S{Yo{`3*APrvQG@XbG?;N(s26F=OONWbHd z(SkxQ1mv!bJ!`{+H!1&2`tlKu{^u`Z+}gyO%~U+%N_s}g7UJb0v_jC*sIqvQ!ocnb z$KYx=??eT)67D+$S1Y+@zN5Wyh zosPse= zI=4BdUTKk)JR6RT61&NNXPnd~xT?b21NGnYN3)zy(vPW29DS2jOOq zA@><=8sTkK=$FvLu%G9q@fK+dq(09O&AXq3&Ceb+ZT^ymP8~XKt_STQe6(w6iEw;- z*;$V-ou=UAnCk?}6JC$gR7&d|YcTC`q7^2F!`i5C8uiU|)VROQ^w&owE>>mF$jhff z-mf(S$x5W(C1pJ_Wa)~y^TzazOyV8B9`TgwMCQO&zZlM3xxju~Z>wJn@hQy-+GQWk zz_#-n{WZlH;C3(n!FBn3@NIB=p%Y^aNgk(mtjdqUmkAPghW+z_@`^LeU8fkjZ$5f) zMK1*N8>c3gsU)L9G3$mHXJ3$Dd`_A#P6CJRXNCeNoB=)^5WbyA@>laxn+?f2`;d|J zn5;}bC>)M>WEbO&8P$TyUP1Bb<85@>*qHPd4PsTY4KL%%fnP1FZz+Ipllnk-LkOB0 zX7fH1@`TS{rrEEQ#lYRuZf|a?W#C5ZRo?M`<50}w>d(N0AQ;c~v^$=65$j_2=x3#* zqQr@~8kyB(J}Xe%Xe1I3_6@lkDTIq#&wEbjmp=ns8rJCkyp#arLziY^Ribg`Due6Q zeHYBqxYydgDG>!VrsrGjqES)hd4!$|6B_#umWcQ!z;%|#jtY0&ko(EcSs_*8Q*}Jj zcEyDP+HIn=M|sYev-?%{;u$Nv^zHZuHnP4kH~s5Tax@Y4H#->Ic{Ao!8(0?z|AXSxb!N(daGxww1m(x*@HUjLA{eemcrm&gde)o6 zOxT5Q+rmR2npawEJSZ6)eAbx12qB#NSWh`;A6b;6_B5Z0HiFj;yLt6@hXPOML(lib zb6w^p^ZGGs7*s1Cdj3Ww2HbBe4NJUEhjZ$l)9Pg|pdlF)IQA+6FS1g&S}J29VB|`T z8T%nz9q?GMkDUsQaram3l4imO^N4g$%0<%C+2$ZiIF>QDP9v8a>4^@RRc-3@f-=hr zTLFwiPQ&We7kByMpq!-jE_Nn_=-M^>)}w;icJZ%cqiIm+J0!r%P6zQXvW-z*S*RA} zrNc>bPJZIP%Dp2@9F^pikjcryXL#~%secG~lzlb|aX*Wf1J?R!eWqhb$+vCS(t>cY zGf;3gfq1sDep;f90>=h_Oq7swpo%rJ$A;*cxKpI9zCn@+yQm9&^q+Ka=bqeK^EM0T zGzPa8bQj9n$pLb;vL+~z z^;K)kf$Q7!a*^IP&%XdvDCL&V$vT(_`^K*FhPEf8T6Ud&e>CZZw%xpS-6k1gi;MU7 z9|(knt@M!P8dOxdv_VQ=!xO6{6auwo^5H=-v*o506Sf{PJ919h6Vs~QOf;G9P(SCp z-&wo>q*cy66NoQ{HuadEq`_GHhuhXP=RE_bx4Lmy$uNPprI^Fai}3d(8U@Y}E}=}a zG;GYkgevO2g-bO-C3S>4h)K{-*ERZfxE<6LC87@?mFxcNFurU6!)%q@9213 z-Ym;C$(afobDDkY)Pmt3z0Frgj@e_caC)|5cOq`lT_3N$%nfW823*9IlR!1NntCqU z9-U=5#D8-nqPuj=rIsU!a4KdmPvU_Lh)pOh_bNBS!#wFVf;Ub8SMH7MlY}QB*Rij3 z66kpH!>12#DpT>>hSU8@>IvvpJMr}O5fezONM2iK;S2*KL8IT7WnjmQk&snmI*wMk ztIS8rqg(--`sYtBcy&@nj4eMC7sQXc+n-LwYbB-75uO6?#^i?zo6;aLQF=?xj|^O_ zbNcPyn}xr2{gC9o6b<9ARu}0Ghrp|7vuC9Y!W%lIBk_cAEk=$EEJ$h;!ZF7RkphPz zcs1;J%#QGz2ftZQ^>4~V(S1c{pG413cD&t~m+&lz}GBdbVV zwFqi^PnWKCD26)`Lfs7FI|vw%x^Qeh1!;HOT3#5IfNO1y;;)b*_!Vq=p~pN4g{eJk z*W<|kR(Rckp_>fQ^Lf5wa5@77legswq{iUXg}#=R`Bdck87sglR{#fZQ`fbd<$+fE zp%)&Nq}N$66SL@lIyPBMWFV|4{kRAq|p)f~6&#X~33pt@BKW z100O#IB6Upjn-;^xP+G_gRg_V#u(}6sB$=stg^d=Z)n@wgnC_}rewt9`6oI=i6-h> zIfRnl69`+A`$59G>@us^L@+DmJa@+06)TLNZa*R93tF2t&b3DSU`V^a(xG~19M1h6 zT{fNx&)5E}wIToi-~|Wc9emLkW+i(4K_<2M5)JNzaS2+rma%IuWYs zJKuk7A>Kem;raU%lGj@;T5sH#iTvYEAL{8OPdOj8W9&^TPA!x1^*$Q`C;iPAbdN<5 zk6(jf&)!T_7HvPA{v{ooLUOP3ZqI`LUZKragwJ&U{rq9+CwUm!bU!I2w*dFev{>xZ&9YM!dwNyTD>WSNjx#X$c8*YG zdn6ijs2+=$NyD@A8aA?8`S?ZSQ)XjDF}`06`(zbYgm zcy%Ee++$N3MnssPSWyrgD`SLzWM$e)iol677gznm>I|hCxN4h zD8*dd8+POb@hTDD^3=O*<$B^BzD3d*j`9##J~$OJH+y9@^JTqTwdXkOb}bi`HDddi?rnyRry&LvE}X80}9t7aAe5! zE{qbd_l}DtL8Pa|vEgXfZS&5R0p2TOy#g{DXWdKzMJL6Zc2*Sejb2~fc(s7}7 zBP&}W$&FThd%B<>fiFE&rPg00=ZD7kx0*jrLK88)t9t!5Sb9e}y!?0y2#X~ywT(qW zfJoYa$U+A8%GWNwtxZDy_LnbyJDY>m2Zj1bwJ_Y)s#+2-PI&wJOyy~ZK-5;z>Z_G2 zgx{+ym+h;_2JuMaty8;*FHEAYZIOkBZF*Yo**7vE-K3@Nb6X@VEIOR6*b$FQgP*&< z&w1h992upP@heBrv{NgdH5t5K2qWJQ^8b`qEzUNW8l%AFsN*mGqL5-f)h&u z2ARPWNSPSo$~bxfu5|WeRk9uGZWymuVj&)-wp;!GMAG1-sYG+lWxp(7=Fx0ovO70?F7Xfo&X%;_6+-w;uB=n9(=1*xGt<%z=Qh2XLo3szO z1c&jwyWxcw{P}F=v;E-Up?q7nKZ%g*x%%Q~ZBMj%eBMTf>JNR{inAM8mgD91!b0n3 zqj7bnNobyFC@2~jPpCv)z>L#*sx_4!SlCC~to+%N`0@_z5xz^sj8r#%C6W(PR{ga8 z+7dZ`GxE*jS71Ul)jOc(02R#@bU*8>W@7cO4^Kv`{D5ssn%R%tq$jXIw#&jUz$b~e zZ!|TLt6|}UMqnmZIMT&bh~H;(Ph!3PbO!ovp4k!ZM{@e!5fkg$UD1QicGH=xcW*q^ zcd8V~!@_`Td&+~0@UcWj-arlM3r>B#PR;>fO}@;QvVKjtE+^+1p-qE*a>w4e6Cc|} z;iQA#JyIaB<;dS3b)m@5a_aHU`vw@`V3r`6Y&iPJH0Hsns4 z1QVNY{`B_`NycEWvAJ>fOt@A3T#M420<#`g#qI_LNHLH0_U_6@y#l@^Td`Egc*JGf z${7iF#`G;(%96?RrzO3nEfaD1}Z z+K9}3s?Mp=T->s78+;#B6b(YnPM-s&G8g4cz!ntt6LG9;uVB@_9 zR$rsSDtW@gK1QBjLFV+liBN2F959(lpx}K+o*9qdwlMNOMEXxO9WAU`Pp()&p8LJR z6Lnf<;6D_9gZ)}0^yoy}O@2&)6$e=FO=eJ_F6O+w2B!-qX|;%FG2LKDCVltrdzR4o zfxdWMCIJ0p;zCH<5lou@+1qm>9+-0DCrfRvAnm@-*wQ27YkU7@jWBsDAXx$fDb$|z+Cit+6awKQaS5$qQ@T_~aUOW$SPJoB|2X#!(2E*~YoB7Z9 zCBmog>mKFy6D}o7=-nUdGtg#VL89dIc-$amSgx(+0T$io^Z9{r;k%e_^6E~c|0#Ad zWLqduT}LMO@LUFS&TYB6Og~ifN^FvPPx5H&d4XJR=kaiIq9vuz5p72m%TxE+g4c9M zxmr{*W?B0+J-N-mT1unNN#8K+=;ttfHbcBTAHS((`l56H^8pvOnCb$Ft*@XCWv2BZ(@syM23&o*p?64(02Aums2j` z`RS|~?^-Vf^bKa|zYTM6%e7q3-lAk2V%4#5G0}sT$J~`hPe>j*c|Lx<1U&V-A)Ne1`iqYlm>X(>JgY&+FO#@PdSJPB3vOR!?>%LPD zmY+*S6E33uBwbaJt)a@c-{UJi`L-JL0Tujq5R#`h##Y^~+J zY#9nMQ#~ApWG<9+;#+U)WC)h)$fe$RMdlEyv_Qcg!r)O}ml5D628>%VO)K^nqe;HG zvt&m;PD;d24D_c$Qnld4KilZwU=W#ift))D2FA2!R+0UD$ZSR=r4UcBw98%+&O^a- z(|)J9$Q;&gNKqBVi)|pF(_d+DhgfHR70Pi%|xwBHg2VoiTE~w_K@(Fz|7vp;AX8Q zZeUr@^6wS>->)nz*OvY71t*Iv%X0FPm4%&!hh-IuhuszHe}A)KCy;mlUdLIa*RV)0 avHkZA`P~1mjQ?H@z`|nme?R&E`Su^55v-*E literal 0 HcmV?d00001 diff --git a/integration-tests/equations/test_sw_fplane.py b/integration-tests/equations/test_sw_fplane.py index 2b3ccb38a..f387138ae 100644 --- a/integration-tests/equations/test_sw_fplane.py +++ b/integration-tests/equations/test_sw_fplane.py @@ -1,11 +1,12 @@ """ -This runs a shallow water simulation on the fplane with 3 waves that interact. +This runs a shallow water simulation on the fplane with 3 waves +that interact and checks the results agains a known checkpointed answer. """ -from gusto import * -from firedrake import PeriodicSquareMesh, SpatialCoordinate, Function, cos, pi from os.path import join, abspath, dirname -import pytest +from gusto import * +from firedrake import (PeriodicSquareMesh, SpatialCoordinate, Function, + norm, cos, pi) def run_sw_fplane(tmpdir): @@ -95,24 +96,26 @@ def run_sw_fplane(tmpdir): # Run # ------------------------------------------------------------------------ # - stepper.run(t=0, tmax=20*dt) + stepper.run(t=0, tmax=10*dt) # State for checking checkpoints checkpoint_name = 'sw_fplane_chkpt' new_path = join(abspath(dirname(__file__)), '..', f'data/{checkpoint_name}') - check_eqns = ShallowWaterEquations(domain, parameters, fexpr=fexpr) + check_eqn = ShallowWaterEquations(domain, parameters, fexpr=fexpr) check_output = OutputParameters(dirname=tmpdir+"/sw_fplane", checkpoint_pickup_filename=new_path) - check_io = IO(domain, check_output) + check_io = IO(domain, output=check_output) check_stepper = SemiImplicitQuasiNewton(check_eqn, check_io, []) - check_stepper.run(t=0, tmax=0, pickup=True) + check_stepper.set_reference_profiles([]) + check_stepper.run(t=0, tmax=0, pick_up=True) return stepper, check_stepper def test_sw_fplane(tmpdir): - stepper, check_stepper = run_sw_fplane(tmpdir) + dirname = str(tmpdir) + stepper, check_stepper = run_sw_fplane(dirname) for variable in ['u', 'D']: new_variable = stepper.fields(variable) From 599800a10373db07fea7afb12b98762ef23bf0fe Mon Sep 17 00:00:00 2001 From: jshipton Date: Wed, 12 Apr 2023 19:43:47 +0100 Subject: [PATCH 19/22] fix lint --- integration-tests/equations/test_sw_fplane.py | 4 +--- unit-tests/fml_tests/test_replace_perp.py | 4 ++-- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/integration-tests/equations/test_sw_fplane.py b/integration-tests/equations/test_sw_fplane.py index f387138ae..7bdcb589f 100644 --- a/integration-tests/equations/test_sw_fplane.py +++ b/integration-tests/equations/test_sw_fplane.py @@ -48,13 +48,12 @@ def run_sw_fplane(tmpdir): u0 = stepper.fields("u") D0 = stepper.fields("D") x, y = SpatialCoordinate(mesh) - Ly = Lx N0 = 0.1 gamma = sqrt(g*H) ############################### # Fast wave: k1 = 5*(2*pi/Lx) - + K1sq = k1**2 psi1 = sqrt(f0**2 + g*H*K1sq) xi1 = sqrt(2*K1sq)*psi1 @@ -125,4 +124,3 @@ def test_sw_fplane(tmpdir): # Slack values chosen to be robust to different platforms assert error < 1e-10, f'Values for {variable} in ' + \ 'shallow water fplane test do not match KGO values' - diff --git a/unit-tests/fml_tests/test_replace_perp.py b/unit-tests/fml_tests/test_replace_perp.py index e3a4da2a2..967780696 100644 --- a/unit-tests/fml_tests/test_replace_perp.py +++ b/unit-tests/fml_tests/test_replace_perp.py @@ -1,7 +1,7 @@ from gusto import * from firedrake import (UnitSquareMesh, - FiniteElement, SpatialCoordinate, - as_vector, FunctionSpace, TestFunctions, + SpatialCoordinate, + as_vector, TestFunctions, TrialFunctions, solve, inner, dx, errornorm) From bcea13da2223141bd02f52f8f2160c5e94287218 Mon Sep 17 00:00:00 2001 From: jshipton Date: Wed, 12 Apr 2023 19:44:19 +0100 Subject: [PATCH 20/22] reverting change that is not part of this PR --- gusto/equations.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/gusto/equations.py b/gusto/equations.py index a7bcbaeb1..f404e3727 100644 --- a/gusto/equations.py +++ b/gusto/equations.py @@ -794,13 +794,13 @@ def __init__(self, domain, parameters, fexpr=None, bexpr=None, self.linearise_equation_set() # D transport term is a special case -- add facet term - # _, D = split(self.X) - # _, phi = self.tests - # D_adv = prognostic(linear_continuity_form(domain, phi, D, facet_term=True), "D") - # self.residual = self.residual.label_map( - # lambda t: t.has_label(transport) and t.get(prognostic) == "D", - # map_if_true=lambda t: Term(D_adv.form, t.labels) - # ) + _, D = split(self.X) + _, phi = self.tests + D_adv = prognostic(linear_continuity_form(domain, phi, D, facet_term=True), "D") + self.residual = self.residual.label_map( + lambda t: t.has_label(transport) and t.get(prognostic) == "D", + map_if_true=lambda t: Term(D_adv.form, t.labels) + ) class CompressibleEulerEquations(PrognosticEquationSet): From e0e309675409df2ed7ec1e62ebbaed56df89f402 Mon Sep 17 00:00:00 2001 From: jshipton Date: Wed, 12 Apr 2023 19:53:25 +0100 Subject: [PATCH 21/22] also label linearisation of coriolis term with perp when not on sphere --- gusto/equations.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/gusto/equations.py b/gusto/equations.py index f404e3727..1d9832f20 100644 --- a/gusto/equations.py +++ b/gusto/equations.py @@ -698,6 +698,8 @@ def __init__(self, domain, parameters, fexpr=None, bexpr=None, coriolis( subject(prognostic(f*inner(domain.perp(u_trial), w)*dx, "u"), self.X) ), domain.perp) + if not domain.on_sphere: + linear_coriolis = perp(linear_coriolis, domain.perp) coriolis_form = linearisation(coriolis_form, linear_coriolis) residual += coriolis_form From 1054c4363345fcdea45a8da254e8871d34b9f991 Mon Sep 17 00:00:00 2001 From: jshipton Date: Tue, 18 Apr 2023 18:22:24 +0100 Subject: [PATCH 22/22] updated comments relating to the perp label --- gusto/labels.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/gusto/labels.py b/gusto/labels.py index 583df1159..59ed09235 100644 --- a/gusto/labels.py +++ b/gusto/labels.py @@ -140,6 +140,10 @@ def repl(t): + f" replace_trial_function with {new_trial}" raise type(err)(error_message) from err + # When a term has the perp label, this indicates that replace + # cannot see that the perped object should also be + # replaced. In this case we also pass the perped object to + # replace. if t.has_label(perp): perp_op = t.get(perp) perp_old = perp_op(old_trial) @@ -192,9 +196,10 @@ def repl(t): + f" replace_subject with {new_subj}" raise type(err)(error_message) from err - # this is necessary to defer applying the perp until after the - # subject is replaced because otherwise replace cannot find - # the subject + # When a term has the perp label, this indicates that replace + # cannot see that the perped object should also be + # replaced. In this case we also pass the perped object to + # replace. if t.has_label(perp): perp_op = t.get(perp) perp_old = perp_op(t.get(subject))