From 453d024cc4e328e78fdaefc628b3a55be6e39c65 Mon Sep 17 00:00:00 2001 From: Thibaut Lunet Date: Tue, 3 Sep 2024 10:31:25 +0200 Subject: [PATCH] Starting playground for dedalus+pySDC (#472) * TL: started dedalus playground * TL: main coupling done, need to adapt to new dedalus version * TL: first working version * TL: still need to solve problem with non-linear term * TL: working version for one-level sweeps * TL: minor changes * TL: starting cleaning * TL: adapted implementation to qmat switch * TL: updated SDC version to be more efficient and allow variable sweeps * TL: merged all sdc files into one * TL: added mpi space-time communicator * TL: minor debug * TL: preparing for release * TL: added RBC script * TL: correcting a minor mistake --- pySDC/playgrounds/dedalus/.gitignore | 1 + pySDC/playgrounds/dedalus/README.md | 134 ++++ pySDC/playgrounds/dedalus/advection.py | 95 +++ pySDC/playgrounds/dedalus/burger.py | 124 ++++ pySDC/playgrounds/dedalus/burger_ref.py | 92 +++ pySDC/playgrounds/dedalus/demo.py | 87 +++ pySDC/playgrounds/dedalus/demo_advection.png | Bin 0 -> 46551 bytes pySDC/playgrounds/dedalus/mpi.py | 96 +++ pySDC/playgrounds/dedalus/problem.py | 213 ++++++ .../playgrounds/dedalus/rayleighBenardSDC.py | 118 ++++ pySDC/playgrounds/dedalus/sdc.py | 616 ++++++++++++++++++ pySDC/playgrounds/dedalus/sweeper.py | 161 +++++ 12 files changed, 1737 insertions(+) create mode 100644 pySDC/playgrounds/dedalus/.gitignore create mode 100644 pySDC/playgrounds/dedalus/README.md create mode 100644 pySDC/playgrounds/dedalus/advection.py create mode 100644 pySDC/playgrounds/dedalus/burger.py create mode 100644 pySDC/playgrounds/dedalus/burger_ref.py create mode 100644 pySDC/playgrounds/dedalus/demo.py create mode 100644 pySDC/playgrounds/dedalus/demo_advection.png create mode 100644 pySDC/playgrounds/dedalus/mpi.py create mode 100644 pySDC/playgrounds/dedalus/problem.py create mode 100644 pySDC/playgrounds/dedalus/rayleighBenardSDC.py create mode 100644 pySDC/playgrounds/dedalus/sdc.py create mode 100644 pySDC/playgrounds/dedalus/sweeper.py diff --git a/pySDC/playgrounds/dedalus/.gitignore b/pySDC/playgrounds/dedalus/.gitignore new file mode 100644 index 0000000000..f08278d85a --- /dev/null +++ b/pySDC/playgrounds/dedalus/.gitignore @@ -0,0 +1 @@ +*.pdf \ No newline at end of file diff --git a/pySDC/playgrounds/dedalus/README.md b/pySDC/playgrounds/dedalus/README.md new file mode 100644 index 0000000000..1163e1f5bd --- /dev/null +++ b/pySDC/playgrounds/dedalus/README.md @@ -0,0 +1,134 @@ +# Playground to use Dedalus with pySDC + +:scroll: Interface for [Dedalus code](https://dedalus-project.readthedocs.io/en/latest/) so it can be used within the pySDC framework. + +> :warning: This is only compatible with the latest version of Dedalus + +## Usage Example + +See [demo.py](./scratch.py) for a first demo script using pySDC to apply SDC on the Advection equation. + +1. Define the problem, as it would have been done with Dedalus + +```python +import numpy as np +import dedalus.public as d3 + +# Coordonates, distributor and basis +coords = d3.CartesianCoordinates('x') +dist = d3.Distributor(coords, dtype=np.float64) +xbasis = d3.RealFourier(coords['x'], size=64, bounds=(0, 2*np.pi)) +u = dist.Field(name='u', bases=xbasis) + +# Initial solution +x = xbasis.local_grid(dist, scale=1) +listK = [0, 1, 2] +u0 = np.sum([np.cos(k*x) for k in listK], axis=0) +np.copyto(u['g'], u0) + +# Problem +dx = lambda f: d3.Differentiate(f, coords['x']) +problem = d3.IVP([u], namespace=locals()) +problem.add_equation("dt(u) + dx(u) = 0") +``` + +2. Define the pySDC parameters + +```python +from problem import DedalusProblem +from sweeper import DedalusSweeperIMEX +from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI + +nSweeps = 3 +nNodes = 4 +tEnd = 1 +nSteps = 10 +dt = tEnd / nSteps + +# pySDC controller settings +description = { + # Sweeper and its parameters + "sweeper_class": DedalusSweeperIMEX, + "sweeper_params": { + "quad_type": "RADAU-RIGHT", + "num_nodes": nNodes, + "node_type": "LEGENDRE", + "initial_guess": "copy", + "do_coll_update": False, + "QI": "IE", + "QE": "EE", + 'skip_residual_computation': ('IT_CHECK', 'IT_DOWN', 'IT_UP', 'IT_FINE', 'IT_COARSE'), + }, + # Step parameters + "step_params": { + "maxiter": 1, + }, + # Level parameters + "level_params": { + "dt": dt, + "restol": -1, + "nsweeps": nSweeps, + }, + "problem_class": DedalusProblem, + "problem_params": { + 'problem': problem, + 'nNodes': nNodes, + } +} +``` + +Here the `DedalusProblem` (defined in [`problem.py`](problem.py)) and the `DedalusSweeperIMEX` (defined in [`sweeper.py`](./sweeper.py)) are the main interfaces between Dedalus and pySDC. + +3. Instantiate the controller and run the simulation + +```python +controller = controller_nonMPI( + num_procs=1, controller_params={'logger_level': 30}, + description=description) + +prob = controller.MS[0].levels[0].prob +uSol = prob.solver.state +tVals = np.linspace(0, tEnd, nSteps + 1) + +for i in range(nSteps): + uSol, _ = controller.run(u0=uSol, t0=tVals[i], Tend=tVals[i + 1]) +``` + +Then `uSol` contains a list of `Fields` that represent the final solution of the simulation. Running the [`demo.py`](./demo.py) script produce the following output : + +

+ +

+ +See an other example with the [Burger equation](./burger.py) + + +## Use a pySDC based time-integrator within Dedalus + +This playground also provide a standalone SDC solver that can be used directly, +see the [demo file for the Burger equation](./burger_ref.py). + +To use this standalone time-integrator, simply do : + +```python +# Base import +from pySDC.playgrounds.dedalus.sdc import SpectralDeferredCorrectionIMEX + +# Setup SDC parameters (non set parameters use default values ...) +SpectralDeferredCorrectionIMEX.setParameters( + nSweeps=4, + nNodes=4, + implSweep="MIN-SR-S", + explSweep="PIC") +``` + +then you can use this class when instantiating and using a Dedalus solver simply like this : + +```python +solver = problem.build_solver(SpectralDeferredCorrectionIMEX) +timestep = 2e-3 # dummy value for example ... +while solver.proceed: + solver.step(timestep) +``` + +A full example script for the 2D Rayleigh-Benard Convection problem can be found [here](./rayleighBenardSDC.py). \ No newline at end of file diff --git a/pySDC/playgrounds/dedalus/advection.py b/pySDC/playgrounds/dedalus/advection.py new file mode 100644 index 0000000000..94169cfbcc --- /dev/null +++ b/pySDC/playgrounds/dedalus/advection.py @@ -0,0 +1,95 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Basic script to run the Advection problem with Dedalus +and the SpectralDeferredCorrectionIMEX time integrator +""" +import numpy as np +import dedalus.public as d3 +from dedalus_dev import ERK4 +from dedalus_dev import SpectralDeferredCorrectionIMEX +from utils import plt # import matplotlib with improved graph settings + +# Bases and field +coords = d3.CartesianCoordinates('x') +dist = d3.Distributor(coords, dtype=np.float64) +xbasis = d3.RealFourier(coords['x'], size=16, bounds=(0, 2*np.pi)) +u = dist.Field(name='u', bases=xbasis) + +# Initial solution +x = xbasis.local_grid() +listK = [0, 1, 2] +u0 = np.sum([np.cos(k*x) for k in listK], axis=0) +np.copyto(u['g'], u0) + +plt.figure('Initial solution') +plt.plot(u['g'], label='Real space') +plt.plot(u['c'], 'o', label='Coefficient space') +plt.legend() +plt.grid() + +# Problem +dx = lambda f: d3.Differentiate(f, coords['x']) +problem = d3.IVP([u], namespace=locals()) +problem.add_equation("dt(u) + dx(u) = 0") + +# Prepare plots +orderPlot = {'RK111': 1, + 'RK222': 2, + 'RK443': 3, + 'ERK4': 4} +plt.figure('Error') + +SpectralDeferredCorrectionIMEX.setParameters( + M=3, quadType='RADAU-RIGHT', nodeDistr='LEGENDRE', + implSweep='OPT-QmQd-0', explSweep='PIC', initSweep='COPY', + forceProl=True) + +for timeStepper in [d3.RK111, ERK4, 1, 2]: + + # For integer number, use SDC with given number of sweeps + useSDC = False + nSweep = None + if isinstance(timeStepper, int): + # Using SDC with a given number of sweeps + nSweep = timeStepper + timeStepper = SpectralDeferredCorrectionIMEX + timeStepper.nSweep = nSweep + useSDC = True + + # Build solver + solver = problem.build_solver(timeStepper) + solver.stop_sim_time = 2*np.pi + name = timeStepper.__name__ + + # Function to run the simulation with one given time step + def getErr(nStep): + np.copyto(u['g'], u0) + solver.sim_time = 0 + dt = 2*np.pi/nStep + for _ in range(nStep): + solver.step(dt) + err = np.linalg.norm(u0-u['g'], ord=np.inf) + return dt, err + + # Run all simulations + listNumStep = [2**(i+2) for i in range(11)] + dt, err = np.array([getErr(n) for n in listNumStep]).T + + # Plot error VS time step + lbl = f'SDC, nSweep={nSweep}' if useSDC else name + sym = '^-' if useSDC else 'o-' + plt.loglog(dt, err, sym, label=lbl) + + # Eventually plot order curve + if name in orderPlot: + order = orderPlot[name] + c = err[-1]/dt[-1]**order * 2 + plt.plot(dt, c*dt**order, '--', color='gray') + +plt.xlabel(r'$\Delta{t}$') +plt.ylabel(r'error ($L_{inf}$)') +plt.ylim(1e-9, 1e2) +plt.legend() +plt.grid(True) +plt.tight_layout() diff --git a/pySDC/playgrounds/dedalus/burger.py b/pySDC/playgrounds/dedalus/burger.py new file mode 100644 index 0000000000..c1c32abdf8 --- /dev/null +++ b/pySDC/playgrounds/dedalus/burger.py @@ -0,0 +1,124 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Demo script for the KdV-Burgers equation +""" + +import numpy as np +import matplotlib.pyplot as plt +import dedalus.public as d3 +import logging +logger = logging.getLogger(__name__) + +from problem import DedalusProblem +from sweeper import DedalusSweeperIMEX +from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI + + +# Parameters +Lx = 10 +Nx = 1024 +a = 1e-4 +b = 2e-4 +dealias = 3/2 +dtype = np.float64 + +tEnd = 10 +nSteps = 10000 + + +# Bases +xcoord = d3.Coordinate('x') +dist = d3.Distributor(xcoord, dtype=dtype) +xbasis = d3.RealFourier(xcoord, size=Nx, bounds=(0, Lx), dealias=dealias) + +# Fields +u = dist.Field(name='u', bases=xbasis) + +# Substitutions +dx = lambda A: d3.Differentiate(A, xcoord) + +# Problem +problem = d3.IVP([u], namespace=locals()) +problem.add_equation("dt(u) - a*dx(dx(u)) - b*dx(dx(dx(u))) = - u*dx(u)") + +# Initial conditions +x = dist.local_grid(xbasis) +n = 20 +u['g'] = np.log(1 + np.cosh(n)**2/np.cosh(n*(x-0.2*Lx))**2) / (2*n) + +# pySDC parameters +dt = tEnd / nSteps +nSweeps = 1 +nNodes = 4 + +description = { + # Sweeper and its parameters + "sweeper_class": DedalusSweeperIMEX, + "sweeper_params": { + "quad_type": "RADAU-RIGHT", + "num_nodes": nNodes, + "node_type": "LEGENDRE", + "initial_guess": "copy", + "do_coll_update": False, + "QI": "IE", + "QE": "EE", + 'skip_residual_computation': + ('IT_CHECK', 'IT_DOWN', 'IT_UP', 'IT_FINE', 'IT_COARSE'), + }, + # Step parameters + "step_params": { + "maxiter": 1, + }, + # Level parameters + "level_params": { + "dt": dt, + "restol": -1, + "nsweeps": nSweeps, + }, + "problem_class": DedalusProblem, + "problem_params": { + 'problem': problem, + 'nNodes': nNodes, + } +} + +# Main loop +u.change_scales(1) +u_list = [np.copy(u['g'])] +t_list = [0] + +size = u_list[0].size + +controller = controller_nonMPI( + num_procs=1, controller_params={'logger_level': 30}, + description=description) + +prob = controller.MS[0].levels[0].prob +uSol = prob.solver.state + +tVals = np.linspace(0, tEnd, nSteps + 1) + +for i in range(nSteps): + uSol, _ = controller.run(u0=uSol, t0=tVals[i], Tend=tVals[i + 1]) + if (i+1) % 100 == 0: + print(f"step {i+1}/{nSteps}") + if (i+1) % 25 == 0: + + u.change_scales(1) + u_list.append(np.copy(u['g'])) + t_list.append(tVals[i]) + + +# Plot +plt.figure(figsize=(6, 4)) +plt.pcolormesh( + x.ravel(), np.array(t_list), np.array(u_list), cmap='RdBu_r', + shading='gouraud', rasterized=True, clim=(-0.8, 0.8)) +plt.xlim(0, Lx) +plt.ylim(0, tEnd) +plt.xlabel('x') +plt.ylabel('t') +plt.title(f'KdV-Burgers, (a,b)=({a},{b})') +plt.tight_layout() +plt.savefig("KdV_Burgers_pySDC.pdf") diff --git a/pySDC/playgrounds/dedalus/burger_ref.py b/pySDC/playgrounds/dedalus/burger_ref.py new file mode 100644 index 0000000000..d04db05046 --- /dev/null +++ b/pySDC/playgrounds/dedalus/burger_ref.py @@ -0,0 +1,92 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Dedalus script simulating the 1D Korteweg-de Vries / Burgers equation. +This script demonstrates solving a 1D initial value problem and produces +a space-time plot of the solution. It should take just a few seconds to +run (serial only). + +We use a Fourier basis to solve the IVP: + dt(u) + u*dx(u) = a*dx(dx(u)) + b*dx(dx(dx(u))) + +To run and plot: + $ python3 kdv_burgers.py +""" + +import numpy as np +import matplotlib.pyplot as plt +import dedalus.public as d3 +import logging +logger = logging.getLogger(__name__) + +from pySDC.playgrounds.dedalus.sdc import SpectralDeferredCorrectionIMEX + +SpectralDeferredCorrectionIMEX.setParameters( + nSweeps=4, + nNodes=4, + implSweep="MIN-SR-S", + explSweep="PIC") + +useSDC = True + +# Parameters +Lx = 10 +Nx = 1024 +a = 1e-4 +b = 2e-4 +dealias = 3/2 +stop_sim_time = 10 +timestepper = SpectralDeferredCorrectionIMEX if useSDC else d3.SBDF2 +timestep = 2e-3 +dtype = np.float64 + +# Bases +xcoord = d3.Coordinate('x') +dist = d3.Distributor(xcoord, dtype=dtype) +xbasis = d3.RealFourier(xcoord, size=Nx, bounds=(0, Lx), dealias=dealias) + +# Fields +u = dist.Field(name='u', bases=xbasis) + +# Substitutions +dx = lambda A: d3.Differentiate(A, xcoord) + +# Problem +problem = d3.IVP([u], namespace=locals()) +problem.add_equation("dt(u) - a*dx(dx(u)) - b*dx(dx(dx(u))) = - u*dx(u)") + +# Initial conditions +x = dist.local_grid(xbasis) +n = 20 +u['g'] = np.log(1 + np.cosh(n)**2/np.cosh(n*(x-0.2*Lx))**2) / (2*n) + +# Solver +solver = problem.build_solver(timestepper) +solver.stop_sim_time = stop_sim_time + +# Main loop +u.change_scales(1) +u_list = [np.copy(u['g'])] +t_list = [solver.sim_time] +i = 0 +while solver.proceed: + solver.step(timestep) + if solver.iteration % 100 == 0: + print(f"step {solver.iteration}/...") + logger.info('Iteration=%i, Time=%e, dt=%e' %(solver.iteration, solver.sim_time, timestep)) + if solver.iteration % 25 == 0: + u.change_scales(1) + u_list.append(np.copy(u['g'])) + t_list.append(solver.sim_time) + i += 1 + +# Plot +plt.figure(figsize=(6, 4)) +plt.pcolormesh(x.ravel(), np.array(t_list), np.array(u_list), cmap='RdBu_r', shading='gouraud', rasterized=True, clim=(-0.8, 0.8)) +plt.xlim(0, Lx) +plt.ylim(0, stop_sim_time) +plt.xlabel('x') +plt.ylabel('t') +plt.title(f'KdV-Burgers, (a,b)=({a},{b})') +plt.tight_layout() +plt.savefig(f"KdV_Burgers_ref_useSDC{useSDC}.pdf") diff --git a/pySDC/playgrounds/dedalus/demo.py b/pySDC/playgrounds/dedalus/demo.py new file mode 100644 index 0000000000..d3c555e219 --- /dev/null +++ b/pySDC/playgrounds/dedalus/demo.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Experiments with dedalus and pySDC +""" +# Base user imports +import numpy as np +import matplotlib.pyplot as plt +import dedalus.public as d3 + +from problem import DedalusProblem +from sweeper import DedalusSweeperIMEX +from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI + + +coords = d3.CartesianCoordinates('x') +dist = d3.Distributor(coords, dtype=np.float64) +xbasis = d3.RealFourier(coords['x'], size=32, bounds=(0, 2*np.pi)) +u = dist.Field(name='u', bases=xbasis) + +# Initial solution +x = xbasis.local_grid(dist, scale=1) +listK = [0, 1, 2] +u0 = np.sum([np.cos(k*x) for k in listK], axis=0) +np.copyto(u['g'], u0) + +plt.figure('Initial solution') +plt.plot(u['g'], label='Real space (u0)') +plt.plot(u['c'], 'o', label='Coefficient space (u0)') +plt.legend() +plt.grid(True) + +# Problem +dx = lambda f: d3.Differentiate(f, coords['x']) +problem = d3.IVP([u], namespace=locals()) +problem.add_equation("dt(u) + dx(u) - dx(dx(u)) = 0") + +nSweeps = 4 +nNodes = 4 +tEnd = 2*np.pi +nSteps = 500 +dt = tEnd / nSteps + +# pySDC controller settings +description = { + # Sweeper and its parameters + "sweeper_class": DedalusSweeperIMEX, + "sweeper_params": { + "quad_type": "RADAU-RIGHT", + "num_nodes": nNodes, + "node_type": "LEGENDRE", + "initial_guess": "copy", + "do_coll_update": False, + "QI": "IE", + "QE": "EE", + 'skip_residual_computation': ('IT_CHECK', 'IT_DOWN', 'IT_UP', 'IT_FINE', 'IT_COARSE'), + }, + # Step parameters + "step_params": { + "maxiter": 1, + }, + # Level parameters + "level_params": { + "dt": dt, + "restol": -1, + "nsweeps": nSweeps, + }, + "problem_class": DedalusProblem, + "problem_params": { + 'problem': problem, + 'nNodes': nNodes, + } +} + +controller = controller_nonMPI( + num_procs=1, controller_params={'logger_level': 30}, + description=description) + +prob = controller.MS[0].levels[0].prob +uSol = prob.solver.state +tVals = np.linspace(0, tEnd, nSteps + 1) + +for i in range(nSteps): + uSol, _ = controller.run(u0=uSol, t0=tVals[i], Tend=tVals[i + 1]) + +plt.plot(uSol[0]['g'], label='pySDC solution (u(t=1))') +plt.legend() diff --git a/pySDC/playgrounds/dedalus/demo_advection.png b/pySDC/playgrounds/dedalus/demo_advection.png new file mode 100644 index 0000000000000000000000000000000000000000..9d181b11b6a78bb6ac696548632fd5c531eafa31 GIT binary patch literal 46551 zcmeFZRan$v_cl6!f+D3T9U?8#A&sJxNOyyzbR&(DlF}k2B}z9)hk$g0Fob}BG)VWi z=6(Oyx4)C`Y#;1{?RCvH2s870e(PD!TK9dgwdTE|yc8}r88!lez?G4Hu8csS`ymi$ zzi(i|Z+M5M|G<-gvxKIzik+FWn~|d_Le9w9-rCOD`t_^(uBMJoukCEP+4$Lb9^JQa zcD8pCWM{Ye-!EXZb2MkCrl)m*54mYCt>uJ3;2WVHXkW#0Un3Be&t;xJReh7PG2`|| z)%37=ce~eQnEZ7%!F@FGZC>ek2ByO;w%W@dYCp zT>1|IAH?qrpT(V>8Gn5H#P6D!V{Xd8-i`%>+m^_lS$9%(ibfVw=sN=O8#DnsB-)i2 z|M`iaPxj6engIA2pZ4bG|NJUA`kErVhM=`g8l4ICmw6*Pw7dU#8%GQR8};HhBr=FV z)Qg3~>HJWCA)u1{|3Bvc8L=KmNp`_8NRS^dDdU+QWLMS-ntu_g@ctxoHnk zmGxw6PNIIdnucC`nusTvu&9Vg>P$#@R!aGu-EUPZ#P|ulJi-Awwfq1r?_gtq8_Yuu`(M-!+o;SpE_A)daJnMkZeNF zP0uOd5o!Z6$f~6g3Z`61$8%SgQl_T#U0q$cROLPw9vzXCI@?taiygbGx#_PsNcwA| zU2JSvD=RCj(tEnQ5j_mQF#TKPqYzalM0;~h!_uRcB)VF17Jf7Gn?waRHUoTajKiYl zj5bcT|MOHdE?RUam-h8t38hiIRw~!cDG5hME?!>VYnYf^cC*3^FFy|cE|rmynRe); z6uvxOtU`_#>U7d4>WYJ3adL9v&Tc1}r_E2Ar;98ey6-K|Z*SYBW1)L_O81<4_57`` zR?qPKD1LP#FkA6@YKLti$_0g$)zt1)&w8-r?ytlSE= z{C;tW)cWp7uUG5*dgo2V7as=GK7WO#3sr9Oev_$+*Z;ZUkffGYBW(^2=VpsH+s&A~})#@EW=YGwZB=BCacQ`5cazR67+RLexjcl8K zx&J=(-QLV%!q$n38!VQ)RhP1@I(}PNVVUpd=j!Umo5*Uc&jw9C*0c3Vu=fZ*Wt12R!j(c|Oe6o)$f;_7vB@wx;WVg9M~!^+n`4oG^7 z>mql1`gO5c);`u<-f)_ro?p{1t@klzch(+wiQaO#xn-N}-7i12JpDgU8rway^tScm`QS z<1XJv3%g4Yt5Y%=wj5d~Z1fu&y})bqpk&uW&KP^5LqM;|jZ&vvm{|%&+&FYe8!&ZS zATl+i@;!G@FV>fnl*Fua-CTIjm}b}>dT-hlY2#=5*L65srfPBf`wA5oH+K@B4cWJE z-@Fbd&Am=nb6=*5>JwBaoD5dKB>o4>^-ezqWviyzg)2cs5icXu5<$`;+R3}I7 zEG&2wzNTCltC5=}hq9dI+>KtG`6ZTCMuu^v zAJ@$ydq}vEzhdV<^_4OK6&0VVQ!)L{*sjWAvhVR5hqY%-(R}Bv{)@J+>V_-#=TesS zmVwztMRK}wnOdaujTJK%x+p9G_Z4=n+kdWPC^X*LrG&S{`p51k8DpcqXZgHiF$VDi zd*o>=Elc>yx+(?AKg8ce?zg^R@UD9QkBs&LS@rt@dAhJk|4qi}1J#qqpT{iU90aAs zO^}>88HgVBAv+94q<#*&&+m?}6s$}h48-UeB~dx3EvPxOY;?G2=msYy5=*vAr+$kh zGsboKBA8&twp3_iS1N2JQ0F__ck7w{eSxZZR{B5!;!YFwkAimOf%ht%-``x0MV>9L z6sib4ilTX4l{~m!QSYfkeTNo~_bhWc{eF3Va49Dkm7Yjl9%H^c{q8YZaxr#(ETx(b1;5B1 zn^6o-czVA1V_ab1qLw&T%QEHW(HAD{^6~x$bq|tjn)&rfPSG3n4@s3n|3xia3~i6F zlo8(x!Zs^Nd)xBJ2|chzk^Lg!MQT^?Al_}=DdrcO_QDBZ^|g#4nSkTf$wZBQx8(9s zJgFPOr(Odqj|Q|O_48K_ZWj}7P+CvydA@Eq9w@mY$i{zV-XcG%;67KhLfooi&F_^m z=em+H)*PGY7VVMzFy`8zT;+|O<%}Up#2fEbpKisx$K~+Ubxk>sQNbNq;n!VhLiJ zN4gXq_t+1Iqj5{)$Yg)6_0zI4?-Fel=Mi_0csr4h#&pId_`U>T2hlJFQTsibGK+-G z)F^=b)M7%aphqhh+{O;aH>)d0pDhl4F2?b^Zddp?KeFi6Y@>UhqaKa!q)^WY_zo@j zP7}UmmMcHH3__$ema)6_5z;ffQkh7y=Iw#Ny3H+4Gu%%SXGELH2d9*6f?N$(&ItXv ze}dmSmwMs}%wyf_V+5?Y5X@%_FaH$(MiZ>SDf3?RxSVMCm;+`qg7_$0b;^17Zqd3u zB_UEyl|C|AXtwv~aWqLUCY9a+lUmS2wFQFK$p*+H@F4+oX!`{C5nE3_6ceBQ>(+8F zd9rCAP-AGYbrF;K=>Z((H`dL@dF>fP7H@2KoSvUZL4<2r*U{HtC89}DlQg22J^Ygq zUB)gib>|8G?X~T`h>HcZ<5zk?wJkRs-iBI*p1c|$MhFj{q1D~%&6!2LyuUd zP8o?d6ba^@>zNPKaX@Ck>7qvHmyEm_=6N~>b?yvI>`{{ z9f`Nad?u(@$tf$UH9tOUL^HJ9{wB{G$CI;i7vUpv&XQ#MZ|V-Z;rt;*w^ew>%_H+S z!P0>6b73WI(`nO0e%sY1N^0h+k4b1g{f^a3S2WJcFgH;M>cCYKO7WSSn#trO#L#(sQO{D_^jKRM!_%e zO8@RWKy0mhr*{wWL{<{-RPWF6>w0J|Sp-<Lqj}fW!$m+c5g1O4rcnp zjmWL+^y`xOGGk*?Zycjel?9R1!~694cN^r}`Y^Nstn=f;?R?bFU8;JU@ZHbJ)YowS z`60zzTV7WF;A`iqSso4J)LS)l#y7!4-(W8ys|5#QlyzrCf_~no2HUNG@(^6-{zTX# zHE2`n#V_}Y=WrRs@3V_J434Tz&8eUYwZ|WCjWQPyW@$&oz=qt8)EF-<98~_ScX;7~ z3l3HD^}Ydva}-)OUAR5oJ(qju376!E#|$%H9eId5r#!#&I5E4&bM=YW6jN838LRR* zhv9y4^ZGh9A0qkTn&Pj7wX1sC;LL787ysK-c0=egLI)}90k$$lh$j_gYg7fJ3gwAe4Zy32XZIn1u>|<-}46Un1 zO7cKrMEC3Q)V$LyM=DXq9irHWkL8Pa6ecF_C9WFTA?nrO^<8+^@%Nua!d5+!mu5al z8ySSB6cb$U6L>Y+6Wy=W7gH5alnkSBFmBz8n}|me`DzBs%>*9oLT2sXU0Z~g23!=# zuw-X&W}WrEAFhM| zY-lPK$a^3u~mJ%=se+4f`h@{qbvE^=3DpB4Tt%MRuSh(mZX=Jt`s*WAb5Rc z9Gy!xl8N@)M-N>#(mNkw%cg9|@_@ti%q8ltd@vE^SKhHQ2q$pN7v5b&dxVquCsoZ~ zI)~i?UvWU2et_n+nv?i8D64V57*JU5JX@`hL8qN*2*Eev2vWX#h*PD}nA1?f9la^;S3Vu6Py=$2sq# zQlmD^1Xi7(Aqx|r)T-+0X-AD`N85~#A5R<>KMTg|BPewuAIYJg)z`)&RxNv6{Spo> zIZ~@87NZ?&XXwOUV&R@;@W`Vg>3Vbk5m=%4C3W`I-~&S@Ick{dR6nug6YvR(I%PKz z8%A9*^anm^+~#;nN=nuDu<(){6TVIqyznZ{qr@^Q{^JrI47O9(R&x=}% z!|5xkwF3gQQ-s=Lsbj*SgXAzD<3&mgZgYU`ka+V#k7TAJNK>8g#S|}Q@z?HQQ)$8 z(P7;uAJGdSM^0Ab1HphfP-vtMsxmt_Keu z2n(V)=hL8cz4Ub7%OtgYbzXrJ3=$p-0t4@ZF)iQi#YO+9r~%UI^#hM!k$Fn@t-bGg z8Z17ONu73eQM;o*Z2zVd4Ca2aS;S8H?v0&)Ycin;*YD;g1CM74taXA;udPmK*7Rt1 z#zf6Dnk43{rEsoe(4m8R+ALp~n~vp`-Xu#^i>rP7p+g`EWgn|L8cS&F=r@&<-A?M0 z_2MRJY3Z8Jd)9M}sUJU*2uZEf1~EtwY^+Q8xU(h?he6h@i5BC!g4D1|HIxd3sIpB- zqQ-#V{}rnUlPq^a9wg0`L>=wEWt<-KmL=u$$))cS|Jjq>F%1JwkmcUJd-t=ahnAf^ zYHiID$%)d$P{K^Q`9NAydKdy{(8#nQj2Zl^U)ja$X6#hM(3T|SFw>~!71jK8Zu~z{RI&^+0^%8Wta1zdU za5KoJ?Z&fy*O#`oY*LXFr5Cyy8lC&ABdbGMSP&yj`#$kL(yj=xDsE7Z2>xABz?JNd7vJvHTEz7I# zSk)gfetmR#>`_%dyOa4;OG0hq=fkJ@Z+B<1fKdN>tFw>5 z?%W-IP3rcF*u(2B+c;2SG-zIJENJ3)5TVgN2#EZ%Ig&OPkusj#; zPd%v0^!I2KFUGjuMx0kF1ift!EVsQ_Ud{JaU_>pkP(U-JCHn)x?I*8H07~Bwf_e^5 zcTO&(a|fD#5BH>*rJwqDENWb(^$X&7eLhBz=8%H5pNsgW9WSsq_-;X39%9N8)LQno z9(XwaMq<-Fl4#-O_-{5d**UOAoL)^pT* zhiDp|BXA{^3{Ni6A%U^|J$lJ}>BTlK!3kNGKQ=Cl9`=j&XJ^=@O&d7>uV~1^Lz1a~ zFR!f9aWKNZbbb421XP{$QY2V*kzFEpsI=)#6D&V>cjEn6Q=6b<=7dc5Z}~sP%w(Ix zggsXpluX(EOm*)1gnr zwQwJiF@s{E7g#rFSrKl_t-WcVGWD&kImt8$1p!hjkBSx6i2|;_0mf1CSOy^xtLyQK zcPIVZNf}Fs#a$~q^`?Z2kex6^0kRdMHF6I}ULZ=?@fR)&VV5~&rxUC)=uDqXEP`*; zyJ$z=Rde|ISs8Nbtq=2aUr3S5KR8NyG44Lh%i}IP{9t;7Z$-*6kMbkAJHHp(KfBKR z6<`vrB&wygr!9J-)@LaKLWyo`P9bYs^I7PjL-s}pcmsYV+a^C?REG>7y2`j2v1yd? zy9)AK3X=VAIq;1e^a8Fd30}uxjmhiECS$+6hLO4Bp0F2 ztV!*0>ZE`=`O86DLXeDTTVC-7y6~2`9P5TGXETM`1V~l+IvD_9VhF?EWBg`%pP{@p z6azkwOkG!{H}V}N0w;>P$dRYWRCQg63+It)Ib&|Pv1#S8=Q>FolgQE5i#tyoo=T=N zVP8y$fSst8>Zs5a(XNNzf2$X7bd0P4;+8;cN0~qO4aYdi{cM zG^k#k`1G8)N2$|1zkUW`-kuJ%f^t0kb7lxE5iZAW7SoCqWbVn|702hDtUR&jEyhoL zApu-;36<9vJ`;MrkTdon!Qkn#!{oi^DJH4a4Eet7kuD@}_7X)At*lGHf3RY&DeA^Wr#?W_^e%a_Bppy*f~Rq4H<^y`_Uj(1;cHO0lg$Jg zz(19Tk#MwoTc+$Rm=PRCH;!hL_2%!OHi=O(IbaYeZ+VrX=MBzEsOdK zlMcQC+tK$F)~O^Yfpz=#-)EIAKdO4~fTSR+G)YBonNzYGuDW-MS5Ws=LBkky%hfDVSy?UDkw=_WvOIchGXHJt&9D-DD^J^Aq;y27 zBuD^tN_X265P0$P22nJyItx^g48eY`WH^(t2$ohYD>r1!sm!_Ni8aDqmJ7}jpG_oF zv2Q0+S^BBZJsA>6P0~eYH#1XlyGC)29g+y%^yxjFUn^UXYe|BCG~A*o@*-KWx@vKSTaF6 zL41~jf{(rrP6!p^Cr3$1x=m91C`s>g%#*=V(+|MR&FF4^E+JpvX(5{=>R;3h_Lce( zrTCcOvD+F)3BjS^MXPAGcYyWPGO9bl)mZOF+18*#I74~&P{^|K>#r4sFXLy#D0-$E zg|J-}l~porIa*l@qECFOpy9#j@yvqwjolFyePff(A!YEcK9gttn4mxv`fHEy)f9sH zD?c`q2v+>{)8Z^=FpmJwC)QBa>fEB}@jM{-cD3Swx@y*lWzCOT)|5E^U6QI9F>Tyf znryTCb?Y^#rdfGB%Cu1kaYhJ4(Wl>Y}ZvwV% z&)ht?nFIU-zy``WzN%FFg5=B%Ib`FZXm2{-{3S^rHo1uoe|0sRIVHk7>ue!w4}jLd z=aIDPfE28Mr7=jM4GIIcXyZ1GXCGo?wpQCano!Rg_uB0 zV^ow2K~{KLwPcqyt|}Ynd-pxfdyvzW9-Ke;`zw>ykhzWaxe5p&d=ynLhoVn?(|TA> zy%Qg4DQ~X`I4u+0yVni{Zc@I2`S|zLUOQ-oI@M-}8;A@Apx9L=0goM{@W_IbvUV z4J&Gm_ao8MY-p-?17DAP;WYV4Q0ubpI==Yt-#-f|o;M)1tE}G-S!9=#eE9hBBlp%G zxWn#7`9PY(*>H~R?$J@x@o(}H>q)NQl!+l{Q-pqjbl=j1gw&AWXaApS$S>dd%e=e$ zma^k>GM+v=DRLy**4iks5_l&h?1mGHCgOpeckAQXx~}$`cG{JE#2>uPMW(GmzNW2l zG&1o#C(U|(B-~{D;kUdzWua?h-}TV>UDi}xx2C&3@vh7TQ3p;(u2s;+DW_&`~?>@)xPvk` zb43ITBG9vl^r`OGbiMv{%S~g-mOT^Oc*Lyl)y4Vw-avtNWhiu(r0)It`_}^Oxi*uK zo*o62(E0g!+}X~~PAuwlZmZPQ{~@imQp4*s0-Pyg%6nSr`r^fWxYeuedTqS_wILrB z@<#LTnijT`$M%UU{qzorC2Ve#p-}a7)jUD)#r(f&Ld~tdZ%NUfZ`Iy*%_y* ze#s#B$>i)L&2x)}n>!Z#EBe!?ZcsNmpdjB3JaX*K4cpf{?;n3;O#u*hyOm$G|10vfUM~~f!4~3gStDI_GIZ)8l?@nnaNQIO`R$emnxsLF^LyxxaY|Du? zSo>yCcey{+v%{t>7~dROWrG(Y?YGqSTJaTs(^jwPMaCk|hYN<8)^7(_0jF&5>jKj8 zE`HRexCr;6sBg;0%3~)e5NKaRaS6p_6g6Fm_j>iw=)x!o>n~m&j+%8(wb$#iEh&~! zUM4dFG9C&h(I-#it(%zoZ7<$m zQ_~66^ZOGuvi_x*C2g+amkJ3g5aYaJ*IK^137HKq0qA3oA@$lzvno>}Kua~!CajFj-|@HC8^87}$Tzj;+s)5f zoNx>FlH;%W|C4lM+jT%^9>M|5b$lLc{*&WF$Dd3+PjZJi@0{R-3Wx@BUDP^ zmOex5zBw7>YkHZ^53q3Fe*A=iS&LqEoD5)^E)#&ZsGn5V+@ata!5^=mYyPByo;lIT zaEQR_qUkBul|l+2LHn9x4&Mdw3yx3ss(uPTNi>yvRN%4_;8c342e+Er^DuA0!hu&` z8}|a)eJ*AM=*)r$5+E1OM$AjISU?zm3Vai7I;)gy3PcBnT4Oa=42~@BiGPmLiMzVz z3Ky^-^0zrafC?^A7q$HbQu6)n`?Y~?08EGcs#Vl%eU>`?F!OV3a;H#mrs|X9=As zoA&IFA^4$;Y#j2Sq6!2Y2mD&Lc)|V|F?pqtk_u~@)bvCgEhi}b2mWT37UcV4wc{uY z+!-Pn*VhIFnaNakRg5aQNtD0G{+HvedOsAM9Ir?CiH4CJ0;>Z@^*=9t?rl6P2~B4L zWEo#7&atssEO&Q@hwEIRM00}?McCEHjY4o?OYVR=C&MCxQL@QHlcF>5Gp)*I9xs7H zRfrO#VqC;&4@o} zb#ax$b$mwAmzAYt#u?H3W&4sdo$*Mi{IQr1AN;DTtEV#fuQOjg=!ZJ_&;}+$K-N(6 z{JR%GzXrK}X>5M@1pEDGc{y`r6?*a9H0=T>t-oYuzH4r74y%>4HDwcZ5d}JYduxS? z`oWiw&_2AN;^6_0>b=z-~syU*lY?H z6wsY8!V^2PA^xoxl8aIx0s>%}fGgaO%WG~GKwIhfO*pH+yg(LLy zrHaR(Z&0n{VP%^G|N46O*fN`7Ywfc@1Td*K+2yPEOzWGq7|DnyIHFq`T3ONN%s5zCYYm(5lH7!v5N!9zBX@DfM;I+^G8k!Pg36v$b zlmQ3r)X%JLnR`I^S+Rb@r%jK()e!|TQ&bD2h&eQZqQL6k53k#RN<)*!YQh;1-dM1>Gc#S@+<+4T{}zRucc7#2u3jXIc=9O-?_QhG_hg-^by?gS5X*(b z1Ndxmi;#+{A69g-oNU+c;l0k$NavSk)b^fzliA z=jzFV{(F5_LSDxlvI&pAX;>~o|MN0*BQK}A&1OK24OR5$ie)4uAZXsKp2K^hs7L@J zK0O~_f?k~~r~5w>pp9?#hos2wO%!OA|EhCihvk#szI{78Cr3$7&-d~#ZIp29#$-9P zncuRj-^GNTre9rMh3;Ai6tBnf?twuw=wWZ_ zJ>Fe}YL@BdRHYw1bKUvBo}0IB83W`C2ng7j52i-7m*P?IKhf0G^mjlVB-pV+rKq8Y zKr!p<>pk``y}i6hM7{XoEG?#Ao#$>IqG4c`cpR*ir?DtfODnnDz84hA(0*Sxs;i(x zU?v8E_Yx3MN-_*b@U}F1)TVvKs^Oe*T0_%j4i&F#xuj?5HRB$0VW?jVH!oaT&Yj?` zu^cwxQ?eOl+ahifYXj zKK+{oaW2h$nH0$M4LrOrup0$B)sax$nN2`kp3Pk?gJG6|v{$5}>dJ^m28}83nrI*> zQBkIv2bB~0#`CR#=R0lq)80oj=08-FmBnRbZkCqHA^&l(50ve_08+I+Gh1dV3$4NH zGepG1Qj(G#S_4blOHr7(n6YmKRVZf#SJm6)k8HOfPs}YA~X~`T#~k>CF1Hs z?a{2~_pRFXUsX2baORd!BF4_|mV-zcdHJe~CZJ?MfASxTZTaJi{ziHD-=Qp34Gk3y z4Ram$l=O53XnxO?QEY;PbrSu zc5iPlD{pPX< zX5r@A+8P-p5{^MvZ||RsV5%~D%aL3KvF(e3ni*hu8pU@3St*Li=9y`v!2_iz)@4fk4G)) zHU8Ze0%8svec|07gjGm+>Tc_Bh_hQYi+HF<5FTYWDt18e(d_IZ2$!6M&zi(=_3d=6 zOTNwR&b6_8QsFloZZj^#paWP`PI-v9K`%vx4KhLO>O3}@Rw~V5Q`GzPLt-K{Cu{h( z4dCt`9)5uY0B5gy)~gW(7@VAJ#FAe!xNgtd)Y(Ea`q0qOBb}=7*#=JxMW}>=0(LUs zd#SQb2K!S*>BI-^0N6~QTNQ^4*)#3E6E1vB(LyYQS7F1kcSlB+o)_lfuBNHsNoEcG z-rijL;OjG6>vY+eajn^cPuZ&AdVe}15<+^BUfsWiFn08W2mR~~`*-Xb4v4VuC=4fx z^g!h3B3(C5C8chntIyXIoj2wS>0oM~K>>1LPPIgn<2TM7dettG(sEgRQGx?tz`t^nC=w=N?8i;77VJCn)QA1sMhmX z%|X?LCBleG?dC+kyi9`7@qY7&EFuWfRS2E;L?xa_wy2GSl=>Z#;Cian!P?jg zE7#`AQm}efJH$(?Jdu^sR^FUt?ZYeF{oS<xsa_y*B;4CF0_Y{f^PzLVMI;V zcn3$w+KJv_i%vt_s_tgtoSvZ~UE2+~yXPYlS@w}gJ3y%Hp!Cta>5mgI+hA;9X zj#=gJL^1W`cl%qEq zHp!npkqwWG$Us}eZwa&rBfexC7QF$!YttV#Lakp|zPNF8i&60Yl+mt$I82RzVF3kP zq#(QWERebd}t3@`yTN(4CAR`_3huw-Y67igXw+3ONkyxc2#Yo7<(22>jbas@06M=+@+*F)pJ~J^V`7R4<9@ z4<|dj@(2Y_j0ftmNB!#>%tWEK$@gO2AVP5V)UA+9)kwY9HA3QRkG1+Qy1lFCgeV9| zytXZuYti$VW0FXo66O$~^hd^X?0Etb2omI}h(hbg2EbrDfJM{X4$9(-`-uUlYH_L< z$FXZ~(eW0C4_jiO%?j4aT!{Urr7zcThNxe=^92vP)pmONEPiU${Jqxx0vlu%Q>bJA zyn#aQkcD-FA29R=%_jaMD{=-{N>&%Q1b<-ct*6;oinv_o>!jIKpFPXi5Ku51x(oiD z%kjA|2aAKX>0K)}v&vA)T?y#&eWrQKOM^V0GikF4H^aimti1J0IN{mV3p9(NYqqh7 zA`ly)Y3B|(c^C*akh)ideg4l@F?V#LgNhqFUsNh`>H~W*to1=p3pb-Bi@+2s?&xDA zERq3YGqI+0H}K)&WOWo_kulKa7kM~Ps{Knv?06}3!{Hs0&M!Vb9%X0n@3MdCDJD*g1YnVH-*lY?UO8-rL-+yRjzxuici z0p9yf_;Ou*d$!^I`}anpUzGs1=z=>;`CaQF|Mk}R9z~+NU*7^G)$tRf)`shEv^)@2 z;k?kG2-Pv-HfwcjS_XfU%XEte)vIXMd4MHFPyf=M9&Uo*&T8OuLJvv)6^Q^9dCZ6uUT4uUGcxbg24kab}HK><*pq^AvSJ7mqP??i3g_Q8nl*Km(kA^{TVr z@<93kbozMDI*Fcbl%GQ_z#3UqC}Xo&DoAZ9SE)L#cIS_`qt$m6n6M)Z*) zS9Eph<)s{Z3V#C~y+!oMVl=qz8-sjGh-B70Ojz0O91RHX;S9h^NUH~uALIM_l8TlG zQ{!1_dB5x#wBVj#teL>70LTQ%@JDCo;^rn4U=VT1J--m>J_B`4`u}mCl(P`UV;fwtdR2)MI$G~kK(9)9EZ~iJZhIT;E?ww|x)Kyg> z!Afa|R0~ zhj9lE#L>siT)Zps8u~PJfBq~4G+kY6qa1g2B-dbnkzY+uk5WrZ>(A1X5s!!Dz;d20 z?jE-@UEs&>#bKTe27=oEomkA5vb(-gD&mJr`oDY4#2~H#?=U<24pH^^_zlcxu77Co zOJvpQ07k0&%0R+%2uAM?+=RL`-X(Du+530aSR$bn@k4yPgrVW@_DwZQzLw|%#b~Uz zLxX~(^31aZOxR0HOV_Snhw+7e_q>>YcQ~GX^<`=i#T5(Y=HjZ^AC|p|i@OB(k^w0S z+$@1juN$TsflX}(iY1Ith<0g>bI{9*J`r9Wah$AEy%cv>H{DZpx zy};3a=D)cI9eqx1xQ_{li6^r6W>LAw+#GqQ*>mfAFi82uvx#eUTeV@#DqnWM4Uv_( zMPI%EQWOJxwOrERm%8b_Qf=pk=&Y$XxI-vOmq!-HT{ZAsQ))-}(?Sr!wB_A(xyCt{ z3BBBK3W2hr&cuX-hknZ1+Eah5QWf;{Qp@b;#dC^^`nOhG&8{!;zDOo*dT3*Io<~vq z)nke01jiwPHVdwY9mmxI|!uZTLN30*2IBSXhV)4Grb?(a_bkKD+AeeYw=lQml%Oj*d#2OOTEfQI#0V zD?mLg-lzLYUS3{t0Rh*RV5(IK2h%^3DRlUY9QpNi3E;H#KDjbs%;) zQ{qdFpoUA0rIu#g<|-GBPr{0cB`!3*t5u2L%K$0Vh|2Ey?r-CovGMOz5X;zK{s$#U~@a2-4<;h!*+ z{ZYV)VQ*!?98jEP0K0t1js1*gai*ND49r+6bF=C5XEKTjDp!YPx;>ko}bN zXqH_Mj4GPX7|qOlknl*W6=?0m?-bT`Y!3;gn$xv>=7$eoEp^AE06w54UYL1SoqFRs9F<6T;K!tTKzB=7T3jI=X(g!@h-<}mHQU=ph3U&`u9G99bg?K0*=|6=Ie@^ zDA+vx-Knduukh%L)_>beyZZY(1_q*l2bDonGHND$qcrRklRn6aqv;)Pa9B-&PWaS3 zKUVq8gx==id1t}vgUvCWpGTGwuIvI-)~`}+_<2}aaV3=gj(!^v7&z#p_APN*kq1co zJ}L_9H%La=KrK9 zkep)^qot(t2_E9c2b>0thiT=G+`f0urN7g)7IZJ5`S_3kp(&dTrEO1aj(hj!x=qbJKfIbXiqp{BNzy+qVHuZmk(Fv0xz{`OotDA#3R zL7gaSzb=!^FJFS7DgBndbbL$<3fz%b{>=LP*>7XRnn5A0SM2IU0l6ueRXQtv-4BB2 zuS+q8j(-w#Vq#bIT|3`UA|Obr?ZDQ$K(_flA4G03QIxuQ1b24eF(v4-W@1$idY+<| zRu`zOJrKxVfSdz@5vrFuaT21gv*3&t@d9cvVW}q(>2XaV-t9*WLY?HH;@$Sc$3;a& z9i5#k@0O+coYPeXY?fA5q}1X!o;z|D99&^Cpd>9LFxB0yw$q?0RBfnRQ2QExXQ-B% z8b$fK14H?H%V)gyTh~e;s=(KBz|kN0p4XUxe==hCIJr+45MNbVwTfPEwzpVQhn;~r zCDA2B8BnH*j_O8;=~4t$J=csaLQ#VOI(7~5K&p4)J1ur0J@3KcaqQ!rNO9Ow(X#Ud zZ#sr}<-S=tVKLVz;^X6^>w9r9IyQ#3;=2LGD4+-bK;B#lpQ{YFl~K1Lf?hTd=|uGt zUL&9-0W>xx4Gj&u%36?q7G)2=JdAXM1O~M#>n0>9!_amPZ8xNTpvs`6sHti@-wv_D!`8V_kNv+- zAu>uxtoRmHSWm`4B?u1lKnO&gjG>AB!@~{}F5b0YpwC|(6#JF}>MUwF)P{C=3e?a9 z^bAzjfj9L}-#)FHF_2Hi`qkj+=Etq>20CX zbj){)k2-n9je~p(K?oCh$v}vKiDP(qenKJq<^u%kS)Vf-9rh`3aa8wx^^7YyEYBFO z`agmJGKfFjV~x4-c;OaG_eHJW(()k;c2Fvao|t645>3{z{f?>)7B!w)ju+|mfXI*} z<|_)3+-~O>YB@VVB=*L2>>NO3g)ataPZ4l>vl#K1%)1X1AZ5sEkf~%~xW>ON36#&5 zsej+WTX(iPS1q@+JwMt)IpYy#EL4co7OGrqd5~1}}#d~P7Kyit!`aO9t z%;hQJ&Ck0g=p?sqKZTH|S#BoR(CGvNQ~&yfwo7b~e5fyR0hLCjII#lsFJ0G3)0b(& zNx0;EC4kycG3B0y?o-IjHl3y>CLdB$ZIa8GK0#l_{2Xx9-;1VKnG>ztP{U7&aWRlx z7bcF%9cwWiH>w|rZ6*{mICO0A@LL&t7rw_3YKF2e2Wjo5G775sRSCnL>vZDjYub$J_-JF$l3tvfYbs zc6RoE5%$(mRkhK-@5TUwFbI(n0~MsZQ4~ZPA8>?Vf$%2& zt>}c4rE7Oky9?2vpdJ?bP_}q?psA@T-(gh^o(5b?GDIT+2b12&>DF&zK7v)eO^+d1 zABRN1?R0l=H+vhB!J9NRU2y-A-lJ&v#pK&uPzW(a6gpTE&20nm?Jb3vGGy@u)yBe5L)A0BE2KsPWt9z( zi(z%S!3s8N53jrh5|rIy_j}8#!%vb-??&~0HRZQqFqGft+{NXHkE%EDpu5)%+w~Tg zr-f3)=v_jS&Fn%+#K>~u;H1nc=dz;xtTnqdE$a2{%6i`{Kd_B6s1?_dY|FL0 z=X!?vjUO@fiY6{d{DJ4R6qKYr(ZbLQ3zY?Pb(N-gw_G=Xt<3Wspq2X4;P1%2f zhuim4JFysYn*Zb*Qe@ZRgeGS&U~o|iF>xUuqG!B*?_);nzIaMo}d{-o}2ML3YyLgjhAW1SRq( z4D~*?p@*LTHNeQ2PWp$}K>zlm^xyKlxA%2e$kW(WhkMJc_D@ASAE45H#@|^CgmCef0Gj^%3}Ob5BVHzUj^n)th3*w_FHHS?OPzSC(pHV2RZO3a ziACXCmgd@Kj0qjJpk8|3uuBuKTp~kk$eI&EHd|a*iRuijQMq;9n>+% z=0F#jscgn|ihGG7kvzu(n+Rwj3x76l{{=s;2dLdHCFt!=ehufj@DA!P?BIT`8wazC zjO;v)mcR=+v&)IMx3TmL4AHFmxFDDT{OSg54p7EFfY;f)8G5|V!4-)m*ETo%Le&Ly z6elq@og%kG8)V;ca^iuV!fvV8$Ct(7b>xm_5tr znZ_7k-QRRAfUOV_IxFE8tc;WdJ&=sBf-5nHfHaW9CPvRPVMLgdq9%t45jvA5?;@Hl!5 zbtL8Fz(LH;!QvD3TJ-5e$zikbneP5+2m=Gd(on&}b@&>pL1j>UhG1Y0if$!i^D{H` zStaAGNb#hiX@>sEHa7FB2*{7QH(zMvf?@IEYtZCZJ=VaYYXBGW3Ivj|uGP|8Gg8NHj+uQu66(y~ej}Llz8a^sgPrR+QKz*<;4u=zfMYR?EN2?xS$TVo+fB#aOfjZH-iY8$Mawc z51>{|pTUQOh%j?dgwAjr=Q1CD1o`%dtsu`6(aZSw3Lp`NBS@~yAbB}r1Bj4bJ%9cj zJtE_Y|2gRqD4MS`X)f~4@wP6zXOIMlC9jO5E5>(a-6_q73walcD9*noeWKDmxsKGH z=GHy2uILfyz5Rh)k)om^7;~W_!8zgByVVxX#qOcVCy=L-p3Tm6yh~BvKORYSc(_Ur zx_#*dM@n2JP|HQbMb8u|GBBrZl}9|_dM~`4Wzfd_>z%rX``O_F7{BtYXT*@;4a7ES zM0`OZQ4B~YZSXRiH?NPDg@P<&+f3z%0y9bn6pJ1+Gh5p5kHHNqH14V6B0d4}G6>cP zwfjyFZjkVZ0Gh_*;5|qp+mPqyL!?Z~r127>1oP_ZnLvIddxr`qJH&OcJ;96YRcob- z?nJ7Hph*;Kb-+stcd(1c0z{}$i13ht)XH$t7bL+@EihztTpM+1MxPyFAlDCquf}fB zO1uUh7|3-U2QWEMX8?hT^YJ_|l_1`&@0th#@pZ(a||;R4nN@nmx6lXzKO_P$XcU9$PdpSzCG z#8;mTwqx8*inP`u?(y+0bu2Pc!nUWG8d`TPb|L!nZB8b_=*g<2C=i=g_U@92cR#H};!Hg6h#zS-?nuS9{ zL!&;FXuAn+3aEH5oH&AYHz9!rGTGtbVRnn0-@i?^4mF~a9+`oOwhrgNt4-6<209pz z#_N`duX;uGc~^y*g}KJdd)5{>V>`kkBX1HDe}o5L&6hBv*nD~OGA$R&fP=Yf43G6) z5PF`dR5P5yH+;4PKP_bUBRl$;cE zk#`>W()xbt{%e9l(Jl1X1I%fUNl{cr`+>{hjWsfl+QAWs(8+}z7rSxdI@xG|clZV2 zXM!Pa5)LzX-@aC3%X0dGssPp+z(y`ep)s{M$pnz$IX2SsSvFR{yq1hTXs9&2?Ip1Xoxp8KG{=GLMx<8 zks|-?K^?Z(G%khnWXdQ$28J(C7TyGtM<^tEo3H@(zuGWBFc8V8s#f7(j-Eg~ct~Om zd*XF)*TL^bo(?&$&3&+)6wiRJ1vS5VV~&uxxQ-y08+_;Efew6;!4&+nn&hcD+2clV!>(1j2b``h12jFd=7#qg{TS9qc2~% zG7YQq9urewMt1q{TCOrE7&-6uD<6W1}QDm5zY*F`^Tr?3!2Vqr|Lp#rw=L$9Q}Pz`at0lQuUv?SA>oCH=Nfj623-^#kgo#pdc~C)=`2WwSS+x?@;Bt=x(@`uu;-M#n^pS5Yk7 z+xq&zbk{haWu?8*x&eQJL?|722s^xR&X)qFpcH&!U1=SIP#=M9S_;y8OvkVDl5AUu z?(VRng2<|fI)SuXfpm(Y+6R0C?E0c`l$6fLm?Gvxh@HmBTe>0BRe_!ua9jK*5-tDr zM~Y%lsRW6Bd{X}sRO3J$o&)heSvg1R1U>-bKIV_U>Ov~W88}hpUN-%i0q(t+95H06 zZcW{o^q_ZO8bwM0NW#ChWdp4-3hPm@E|dF~79l`|wGXDA*{|EMR(-*z7mtYa5U_3j zm6zQdwe9LxN7>Lqwk!k3z;cnbipnj3Uyy*4BL$8ZIKsIRr|$0B!rWYAfM-SU3*jgE zf}S<`SDZm(04Z48ig^@t_eC%_t)78ul#IACcK^y3>*p5~sFqkzA!^O?*gmXJT{u(8 z>@<~s{msBWgeFbqfZo3@*XmOKw~Yj|tM|_v;HqG;l!{`R%A^sAU+f)Jwkk~wso;x= z+AP^*>gR&p;b&?p*jTz}TMnhcLxNy}d3d7q6`P^IC&|b1>c~lT%56Ak(Js4sAXhC2 z$jA|cT-$&MzMzCSJRP(tlpkpte-Rkh17Q(}27FNBwD>zyR4gkaK^VD%{-~{VQul35 zG={b}6r5|JGyA`+;s zF2v!LAJ${lpcg7%1QMp7-zCI*Fyh8OMM(__2;Btq-E{nGYu4kSNM}>HzBuM~EH1qD z=bq81!;zbaEOdd?1tS}JBOeh<#SLg*e_mSi+FOl20C;HH%-WTOWM^-*^Yj{>gMy>) zuRRo@5dpYk_@`DwX_ODR_RXXsA)gJz$GIALgJeZg@1yXQ+gueAf+FJYvZu7C;xoHq zn}|P>?m3p`hUgo-!OHRkLi#m_|Eh53T(&T)u-JDAIb!v5PONj<9~bh4EKLwwMgB@v zSr-4>rBlpciVq~Vgm6Mz@eb%t5-`P}xSFlfuiq0Z=+SLV@j328QYSZ58XFIG&vcE8 zK=<^=GZbz~D17!U5GEqoH=dt!x>Nn(!`ARN<*Gnt1s@34yeDc!t^SQPT8M$3jf{-0 zW~BG(;uUBu8tiXeiJyMbZs@c7&;=Db@*)kIkd*T;Jum$Z9URjqB>5Q}nZj=08C5PRUsjo&eNpy!H*eD?m+Vtb|lD2m_ zc8aZJrDzX}ZoRZ^aui|F(Foqf*RQ?c*Jx>FWAh?nl*O#W-;2;NyA*F__Hrj3mYy-Z zaQaC|`|w7Ms1H*p%7o*#o}S)IS=oo|?0vd}V6CZa6(A;V_b_~A|BJ99A_IN`H3&$U z%v?u5hss#efF6>VkGgo#yLu`~DU?4(DdQJd)GRg~i}ub4US6t~e5kGPQJ-f$L)Fjr zNMu`PNS5 z-;_U&*`wOc1kB=)t;C9Q|DgQ7Kwu8pj$gsb#G(n>E3o*nZAwJ1Z-~S~$HC)^k5Twvsnd|wZvaCD zR!0gq7B#Lg4+DqLqD$@UL766xS{A0uh`1A)!l&{>Wbb=n2uaE{339I%33~F|&+~oY ztS!8Za4zcz7K>btUhJRO5naQD)^$(J&yz34bjZr==NBTjWM*t`HR07z0ZN;Uu6w5o zYItd+tkLA<6eTpeg`oJ90A3a$)kYebk9y7%18mB(xa}!o6G1!uxfr}QPbzoehdOZ&3(hvI z+ksk+SBi9($JO)QObj1CuPY25PJrCvsKla`=8Q)N%jHSEE=?I*$m8Ya=Z8RT9?8qE=5d_GlJu z0&i>|4NQ*EbcGetQETpk0#jxeq|KI42O1Wm`oAdb(g0yEgYQ=gsu5BtqoiEo zh3`_8pdd_s{!v@eXw|;aLTSRy9!w$r%WpO;hK=o?DBCWa&PKfj582uyUOMtDF-hpx zcMf-kK63{yBVu4Y0+7=&y>740zD)YuKLrS*fKz@naytMtbnwRNu-u9fzEjZBdsAGr zKI)-``qsbb`hpI)(RLE9Z7dsVMDci0E&}1x!+Pg-aO{E>!5m@H_*Us(f^sMD(EbH0 zK8x#N|3d3~@UfdQi2_JPfw@c_<)^V>p4U6-%|Cg?(gPX_PWOQd&>kYdIdJ9bRS^7m zXM{kM{yPw_fupRXgsZaGRr&9y96axwKG*4Gxx^%s-_$2pj2`n_@hVOY-oK zEnYDv{Z1i8<5|R4(VLQ>&Z`l6-Q=yH*M9qbbg3)u6LHvJ{)}<@YJ$BHE(5$6n))9W zNoZbmDbenjkv!8#%Tbpa7S=w#FaEHgK<8mbQ1*Q5g@uH>da?n|YZ~|%vVZT*MK!mI zFpGP6+{(~=F2IFzq0_0pnDnO52*n|yBtn_iwP_WUpvDR8bot86=K=NklQ{9#<%yPY z`IO<)?=1ieK(R1PuK!BT?`#kkenLz8a;`=pH_h!j>S?pL*{a@Y4o}~bxxSK)V*Dq? zGPS_L_`{?90s801Pjv|t$UA_1WZiBkH(!2?ZL{%V*KlUh zN&X_4-ED5CT&H``BWn&*6r^LX0UB8|DSx2oh76@NKg`n~<{0t+8sOl9TK zw#wG__U6{sU4hD}%ccJx;Q9-A+N1FSG3`2SAA+_QIqfI~F}t&z>!;t}|8G)#Q?c%8 zNItR0L~Kxs!9%*QFBN#}<`s3*CZlQMadft#TAm+1kLwYtOZu_U3Eb6ql~Uihp$KNa z+XRI{_=!dtuU#6|T@Lu3Urg#VfmXowB>Fe`Pte|vReu45Ds>d)Icv=G`ex0xvV--| z-712jE^txZoH#pR{AK^Efygga`EAZT(VGVgrm0Qj)(F=Af>EQna2as|Beu1pJ+)nmL2oCpU$xqYQSIJPpOAwww>8d`}X-J%AB|jwHeXrBQHAVLBq3e095*W>$~UF zXThwt;;^}I{nY^_Q1jWe%$f%%KlCp%VQzQO^nKbcd1~GneZq-w$KCdCMVy)53Jn=a z?5V*PQyc71FdQ0)2ac@DGwz~RyOdDn=mka;?`t=&i+E(Yz?3}W+ad!BMIM}R>CN&j zk2>z;um#b3;}7zE7}-Ht4CVX4wG4X{6ZV5Dxr)>hK@n7EyY%hbw;Bku4+pEy8llrl z(#-4$q#(uRoTv(o%iJE;5#?^o31}v5+lj8r55z7&L#QvHoS^%>_`%2c8nwvPkWxv7 zz=p$cMT&#HZxOqg64W)4*?L=L`n*9M8E{$yZXH30Y45irE)hGYr zJF0Wm#LUtrUTD{oYuqq(j$ak5iUh24om;&^;Gd$o@BzpyZe-!AOE5S_Dg)DRsd^JL zr45%nS=l@>He<22?(WI{d|vf`3@Z)1)$jh_sQ+927Z+_F*p*YlcK0*$EKW1h4Lp+OJGl=E$yZ8`dC8$?Z_?p9@{P(xn+54`1l~A&OePWE>~&$i;{wg=JVi4k zMWZ{LbiG%eoTq_V1|479_^wYIP5el#GOttFskt_;ZSCh62_Um7w`?~ff{F1sspB+} z3N1;3;Im*mv!3A7#$(llg7~(Z%s>x`G#ookcGd73OZd|H6Z*;l5~u@oTC>)Xk2ZSC z!CD!qfTisg0YcDi{eV?%_|rt5(?e@*bFO8mRQlYaQXyaOGrd2AlEt9AnX#qxKdtM& zT0VYptgGdI$yMQC6zMste!^5$6-VqMRX>F%nm9->Hsm+o+1?Ud0{!J4hDzagC~yAf zyelq15Td>v*DyDxk31FFEEWcWxHnFPL*G$M9gXa0-f+Bp{AIs@%Qfr17g-D6N$g6> zF7|goe`Ver0dQ{<2ki2{swJceY4ogJ0mPuBUssuS#>fP}$avGG%kx%r}L^WvV{}gz5M=*t!DUp(H6D+b*Lh)>cE*gsKM*O z^se#UX!)X2)fko$yzU4iJaPa0lQN6T2|?KJ6d@KomRop6%pG)0km$(#;d`j(ylO+O zdU^?3cEj)45M-$x^wEvdgsSS!=#@(#N zh1IyWjo+wvachrperx!HvUz1%x*JL`8==H)dCg$}c=u;oXT<%)c{#=jz#hs3*3#QZ zh2+6^i-S$>^Vr9q>ppm2Ej8$21j^F&DY9C;WZ}uw3=ch!>$-R`|)keG>hQ37Fn@~|#1T9*n~V)lu^5~J)os^#Z#U{MRG43hOMaDls=@9?jc-lo>3oBfb7phkhz`xvXS zCtD?Rsz%okUq87wwODwYkW8t(Q6s}PFTxXQr@Q^PrALqZl?Hy2Os{9@uf?4%k<YJ_Rzj4fqd7goa*9D8>!zYF;xHtf4yOhlYKGrYCR2F*1);tu8Ikry)G6 z@g<=NQ4D@zC&FRN6kD%3wjRQqDkF8qCd1S!{JCsqfsNJrE&$2fVzuU@ce`F2*>n&Q zh+p07?@P6RWM$*jNJXMV7NQn*MZ_kP*sJ&rv5%*FWZl;3MX*x3If#26@ePgs9+*07 zj}`dlCqlC~7tX#^AkCy&DIoMH9-tJg^kQ6d767x3EZ{7NYUR6m{r$L|zg&sqkzR*n zT=SP$?S_y3Q`Yr|y3<5Efy-GAD$|f%+2KI9f!3|9bA89s;u$?#Gqp6%IzU=&$1M`I zeI+WWbIS?IefWIrL#D#z-Kve|(lApEK05d)Q_Os> zlRn+xT7hNpP6&|gL!qxpAuw$zv!V;gpbYCTqzjJ{FQXn`^D#F&!cTHZnp_i~I#GMr|-! z!NAdh+^f%1jyE>!gg7ND>XRHRn~<(R(2o1_5MeSvVO+g^UHC&Yt^jV0`pCH6=2RoS z_^gJV8iA-xn74Ja?&Rr<1vamANszr(&fsl@mc~{E?Q+(sgmy3q-jQisk5Fv zN@v@cEO{K&p4}HZpR=C|sB^C{?{}X$AGRvk`(LJG`y)PQ9J3~}iIZ}|>jN`BM5jM( z)IDJSqenz`-;aW2Klx((pNLF+4Z^~fXJi?4jyLuO|fqV+c*AoZJ*2jMY z?6w6`R*i)Ch@CI6T$WbxDjBBb^e>5wD6w3V_9_XKUV2fvtsDbCYypyHlC?(JNgEL6 zMhjM`TQRLm)46~w_fLrx4&Hn zNSrj*atmwBydD=rCf{bGol1x)<0O41(=5z>k$WQ3Uc7DuM>bO51<d){j*>l5H_%JJ#g z_dN^Zw$$#qj}%VU%2E9#TuRdDs~{QN^RLgW@t!bqj6U*OeKvBP>Y8*``am(~51xv@ zm23zIWq+WjM0$Hg!kh@EUaXzYXyvyS_+9-vQ_o7%YNAaBvdNe2?X8(r$%zz29fm+9 z%DMN|)ZZ+p9^$V~dTgVVjdH~PLuSHs8-TL^9@lX23d_9l!5!^{nb&e!(!j2+lx36y zS0({~^whD{Z4Z%)CZEpXGrmWxlr+qPLuaeMq?zHIqS39=b>o=7 zmvjQ-EF?7U4|+Bee?0ct9$(zQv`1)YG$j+Sar{%#&59?+-lOBU_eG07>H(sEQ=9Q) z$0b}@-}aTM&rQx@T)KbzxZ;amJaXEcP(Ml>e5V346>_+#<@kS!k&_GG$d!_1VOUdj zY0)y<#`2YRZmdt&e^HGpzLq7{W?HYEY;R=>G+%_m6Z`U##3bB=CP_@=I(c;gi}{aX z+H?z%XSam3)BsQdd-JnZE2fy*x4Ge-vvvBgL-$so> zn0X|=ugrT_mlzrZRN4>2PiQ``w(zqR`F%n;W~bDT~||Kq$tRC1a>|{Fg=H z;BOuj2y9t#as4_RY9D3rivoy#z|&PD*!~xzUJ)3v)&EHfb>?o*eNw4uM5rilTyz`x zKPVRh?f-30wnmXCA2SBQ&Q2h?|99%)e@v1zF>BOg2)op9Ge zDz$%ex|{5BH(oB}8n%1wa{Tud34~eO^B<%AZhK?s>OUvtZlwV-3mWpEUolj8wJEbE zTu`5Uw9|h6wWc4}6pj>m4`D9;MdiN$d|DiZX(cZ*dTc=O8X_42btO!u6Ki`mH9=5D zabOo|keY53k^mK3r8z&OY`6Gn`rM2ezfWbwspgxYig~f^y>21;Z&BB@O#%$LQYvK3 zi~6YOd7nGH;2Kx64ABn(#23kgMW>#)(&hI;S_00_Mu6!d-N&~SkG1lICB)`Bgi}<{ zoBDd<^++d2?^UOi%p`gU5+IkpH9B5~Q6%jr-aJ@tmU~?a;~tudsWKd1*5SbEA+KCi z@wq2UTwmcUGS?=EkpSF7qe+zQpOb-`dyT#3-qjrHoU}P%$vqqvbEzA0wN>=ue+6L0 z#G`xJCTu&)@!Y3E%FWqc`~Ksn|i;EEjA?GJSA=W zyyvMs*3WoFh`V3@X2m2#RkmWR+asdt z@7ycis8ENu;jgt!?*piDpJ!oIsN*fy=IGG^OBBIZKN{WR(Q9j<9uu)-t#b_HmLSdn zPs;-sbdXkel;t}2bdUmfZ??Z$JG1LPXUW%V!sslX{M7R8kPE}&7j0lJlyc;nJpcEC z{Lb?lap>FLL(Vt*=igGDQS}O)iyPo=t3Fw;!6NUMi2I!a+zuOmHqPm37{!AF(4b15 z*%#0N3S%e!&4)UNOi4ue+%Y92)nnsdr~KVVE`>YtIwrLP+s2$R6>yV5+el-;e?x>O z{N9^LIDZ#%mT_waDS~qEWjr7A>P!-UP6YWh0+#~FF1%XAO4I_X*e>!p<++i0f8l=% zPlpx;Md6caky$}+A|<>GYc3ajijk}pX$nF$X3bCCf0Y0g!uN2wCvzt&JSa_zDnIPzq+ zbcwzqdQtGMipCFK6~m!_8c_}!VE&v#7E!$ zvVB{;w!olW$Y%NgLC5zkBJYiz+Wx4f`1clDA-M}@E}n+?Q8+b56~gzDIYh_2Vut?6 zMCb+xzX?jGl$iMo;N9?rcC0+HPozfyF>|%Eh-rC=O{c`R2{6y*YI`Cs(poa{izvoz zuU-qG4E+idU7%CG_Z7^i6Ib;7Ie*^jU1|iqxJ?SZC<(Kf8;+NR=EEM=%;eVV{`0Vh z{}=Q?0Q#ZUY_c(c+2wD^wp?5ic$sBoq9wz;frnPGddB+55U{x%u`tKypU2M&o|aTZ zaWCzjyL-x7h-Q#_%tj7w4aL~O_T^nYp71m;Vj=S((fL90j?90S5~yZ7>#7SIu$BMG zYJb4Q66oqO3Pq2kXxKz`A>0lka-xJl6=OSCz+nr4FJ~+)q@!XBoKi97#G?JKPzNka zRD0#tu?K4HcFaqKZ`%aXyyt0{se-;1>#( zZGtmbUUqp?q3YG?6|aZ5LQ|~ktkZzO7JQAF92o6CWkP8U?OgdU=Qa@t`Ue0v3+LL#JCBuHWuCp$VSgAtmRN6Jy zqldCT{g-^WXWqaG`(wJMBC2NCH2Mi85b7n<4lku4`_c1rrBfT$ zttQ_~#8l;^uk&AZI2SLE&B-75qwYz0Ve-8k)uo-yH@tTj02A-DlHZAPL$ik7YtDeo zqUsU$QMr$L>wuHq(qk~&YIBJ54{fZM99JWQBrYs@NA-Aq0=!)~yi-5i?1Th_Hmp&Q zf5yQ46!_=P09?4w%7*fmuwF*L&#LgWyMmczt?&dFAq&htEF{=pfS3iskq7uo*J@!h z-YuysIZR_v+A`Xh)t&r1YLhSdd$T_wQ)Ch!sQ(A+d)RXupb+@&ICl1O)WJBp9UX}2 zSzw}&)MtO#vbYT!GaX)b)RBC8pEY9a3}X6+Ej1xE`&(Y(P4>^$-yff9r8PG>K1+VD z?J&!u86XK9bRkDUMcjquA9fs=~aG#g$JXx`_`KM5Z1g?CuapRib@k ztoY{kmcHwmgNM*g0W8%otZ{$@-*P{^CDj?ZqjZ~I6>=|B^2k>u5PUZG67bT0v-hwI zgEu@4(4%T|-GW-`og_`kDf z{!x1xR~85-&*HzjQ;ZE2Q`%&#v^Xk?w<>hysS8d@fJB{ossRAIGwqVn382y?XD)`b z;o7uLFW%>7BbhN2uK8K)e58Xzpm{$`i^7_n5upYjXo{z@bd9U4L%)dENcpH3BoxGh zq}x-SMDX;cAaje1e{==iez~poNsG$4n$qGKfxF}`yZfiWZ;o{MtoK^46k*2$#FFko z3CR*KS-;I5^^l+&DD77%*pQ{`5`y-K*y%?mQ@-lpC zx9P4x@A`FA)zCfx@ZKR`bV}_!ZS9DYXq&pgDhZ9m1Q^+vh# z1z}$Yq0eY)&kUS$OQbykV4O) zt5Ka#4Ft8uTw6^M+z`u7eX{ebX_jp3nKs#@K_*?qLSg~O=lzfOyJ`QT9c!Tw5lOaw zezIIBYw;GfF}R1 zAmu8+JA}q}S({A*`{+suN*u}+5nLYtxc*74b5q5MHEXMS=J8s>isK0b7`Z6F$}9FK zBQpXBfOMj5SZjW-zS?uO%}8GCI%End0n(s~lxze^4&S&rEeULF?QF;R_>;C9Oy_lg zs^7ZY3UjpbV2W;o&YPECfy|@;Q#QUSMl?pcR$rJ5_6SNGT!M+OJXZmJDP7U+P@Ip- zvG5=X9pZJlM%j3@I;vXi!vLuxOaQ1k_4Vv9^cjp)Gm#y)mQdGwb^=!_c(81LV*2$# zvSATdl5`$_t1t8Gdf9?}mn zdkcTD^PsreWM%?FMJ2w7GI{7FEY`xS=1^HCEP5+T?P}Ww``vWB?%JtJ>_`!M2Dt%g z`GIAoFyW+mKEF6ZOT>{Hmn*9Mq^K(&)j&RvHE$pLlOTTfJhf*NHh?h?g{K?=zzG2( z6%~=-&NYNYH&SA01ZX^HH;)G(@Aqww|1)O_Z$IAAUIw9kM9+$=L~A8Q=BWM4?XuW_9vv zxkL6usV^zU`l}xoVE;>f<0L{JTmH!rV5o1LPxcZr){m2P>iei($u2m)i!H`l!BMHSzL%%poDXI&+wqeB zpTSnq?Pe^uSZ=KxGmapzk`h)9W2f9fZ-M;z>+@^zZHkwbdwR7`r|=nEi#z@hPF3b# z7)T*`^F~hED~ONT-+;W0gc#BvgKY)n_&HfvOSyOzT>uNKka=~UBkp13me zmDTr>H_O&zsoaX`;;IRZU#Q%3QKK3Y(pL{M=-L!~BBR}2w3==Ihm&&7%^YxMlRYUqBMfTp*HT9&k^>wa@O*yI3~bPLB!RHyz1Rk;26IVbz;8(d>c)jSZY;DYEht; zcF^kMc2UPM1JU_+_!6IIWBlD?glldFhS5w#zvx13>};e;1WH&$V9JM0+xR{#w7|H# zuYP_aseN%hBa#2%c%|DxJZWQEQ?2v3&u8ShbTb5YY){8XFs(~*j0 z1?;l}&4jd)YQ*r=zZV{=1OPTi-RkFkRm4*~Yd#fp&oi;#*tvAQ{&yp-hPV2CFBDio zXF_F#roeLnhVF0#poOsIfCqz@i$5U~!_ne!0t^b(QY%`(nZ1QRfcEdap};d?xBqPx z;{F@LuC+x6JXAqRN0W|CwSNJmZ1?z8#|tC5woNAMmY;Fa6s(j_cV3=UreFh(`U9rF zL9s{3gcfn~QK1jB4OSqbj~`;-kiE7>)>*UozNV#DJlkT(pY08L8kX3W;se|Qz4;=> z92_f6UySiY5}Mxn(%=-=`K({S_kwIDyV`^7;t9)UCPMG`kfN7=m2Kf{z+I#oBf%FJFyqYCz zGY~MK_1S0jKg2>|AqO;7AUG1}RF{gHf%DmI$cg>Y50$-21?M_-y7mYQOc!r3R&v6~?z-ZL7m@o! z4tyyL;UFZQzxb@#3TQN90aVD}izEyXmmD zC=AyL8;ICy>6C?=wFbJO)S|`V-dUf&<#wPF%q6+e7^3so{Q4~qWZc@n*Y)-eIlX#At_DX zvyYlI{-Dmup4r}|n{0j%QZhW`EWt<1=u~^PfAf&`d!FV<$ zdz(t0ba}hntcQCLWP{GHz23EDhN)gQL zXGuEg&!CupbU4T9Rs}PqGezuE&P~U!!Kn3algXgJ`o+WcCS@t7oTw91 zu0YDX?US?ot8?emU`Z8Bf@Uw@H9jg{$$;l>4?P(%D`(^TGs)Y!FhLSBw;xb_>*uBO zgNOPzy8Wc8V#@s|vl5h224<|?L$+7*%!t+*g`Wr1O05t4_`0&_aXFj8l#KijCz|nr zua6< zOu@*>of@?I2<<7e2!)`HUK6N(dr3ch6_p0XMH|s(3e&yHErfP(M-7*_i-l| zkT?d5)jkVNjBXmLflKUcst%d-!}JGhj}=xsD7=&W#G8ib1T+u_gRHqC@MdnkGF|aF zgGw{1rf4urEGqcs|r|ap5O7GH=K{Tqe zWJd76nP`7Jrt43Z22(D6w#jdbM86F)%;mLYR||+R<M7qJCsjO>>%1N=8M9PV1@gqjX(F}U44|zcJ2CN--@;kBwtQO$NVU% zC0?Iv3EvFSRaP`GUGA!bd+%-`frTVUGzlvpf-aL+aqf0fw19KDuw ze$>)n1*(6r)7^?7mcSMx2~#I6T|>|auG%L;ZTT{?jix~Ea}lfUt(d=zZ^mv7(Bs^6 zZQ^jX-X?Nz6R|dDU41q@G%W3R72dm29haiw6;y*Nm}dmp)T;oz|;H}rzxgZ zgaN=*}$3wg#RkJV_TF?2MyUb>&nLXlc=g(7?DvymEaIFzzl6OUPunM}x=0zygEAU78TYi19}lU}-@E zW8nuQscX!XQpsK?_kh;g(s-v+Lu)z6{h|Ve6lhGT;u-T)9`5FCrqFB%l_Gts(Mdye zv6V2~^)JwAX?F2tq#=`mx-zei)6iIW;<@~&N8iE7lnx`;L(-txmy#{8!s>@pr+y%4Dca*Y%>CVUyEN=7+?v6on+I*MtS} zwi$+$(#BiOO{HsokxM@QEtexSHT07bddnCPm+Z6WA_bcJFah1Rbccn}P3`3nnJzb& z%`0GPT>+XEU9^fFsaQXFw7m0vCShnLB}C@;xA6ZV9Ki&Ct=q^GW8%``1DkCFY&0wVWj! z=J_X#sEr)s+dJK=*I&PVVq^2Xn05YSum(F48`$$UyN?UFpbSCT7bf@iq}Qi{1_w`Z zRI>dF3HBUpIFgv$1n~v9N;y%c^7i&Z1DFh39V0$GhbB;iPJ7nlZq(1Mv=@J8Yv$RQ z3bOsRpi=5et8*wi(yW-Bdv)a7x!3>N4;D8IH!Pl7yt4=;CpCOMGJnCPEJR(QlUdw$ zxiFw3J+|d>3vIXHw;#cx<`jvcr#AO}6M|-StmL8uG}cZ{+o^ANhs<|qv_Owmai6p7 zj`pqLSz%wSx2nAity)~D(Pv({v_$zQ_0>GBb%Ix0n7dNrJS*s&kNvEuAESIaZfW2T z#D+#vs`TFf$Qpie#pwa+3zOp>a)p<6 z>};MV3002nNR?+3pw`U=#B^1C@Xo z-^b5!p!2Xw>(VnKYC;{OI|7`USq1gixvY$GOfNYo^UA-SP3~xs@vrThY+2+=e*bBk z<)u%H8JeE2xI#bsEd6&YhF`1sXG|&67T=NYuH(zDbH2==i1 zvXeaec0N{HR!eb&L|iE`sh7D>J2WO2omux5y>D6fZQufzt*y=TI#PTlJ_`%o!(sc) zu9mlRBOXjE!=)^AlPQ0t-;x^eJ^1t`OSbjCU{07JZnzw4V05jKd#g+FPxGfap}N{1 z-CyE7GRsDLx?=F0Qa{BvWo_@g&L`uM9(dX+e!8;5!eWCGJXGxa-kyRQX}V#61lt;KUL8TsB@NtBiJ%wE6swI!2ny;)SoU-$W!>xG-f;!8IO>jx_Lma$d&}8(muVXO#rzDWZrS*sErohL1}YG=_M` z{~dlRs8281JXM=qPSST==Di`SRh6h2R~b{7#BpxyeC*%B>Lup~Zd+gGFuCd6V#$ud zLElnLmIHbgiB%l$glawxqzM~(u}~K}mT+P6%THz=WtYw=*0*F~c`mZIQAAB*LuldT zp+DcvZC}fg_Zi}!)*TI@Wq)XHEitY%C9<8)1F|tRg-G%++@NCRw8<7`IFo`V`^$0=gBEw^wQgbjgPZLEy7mi8yq-m zI0H}3V?WeKUp#S)G7 z9nP``WrUcD%DQ&45qXr8la6G?z9am$tNSX3Cr;UewiikPFDCoL{0*D|uH_4T7hErd zh9u~r*#=bX{WpF&h6$4x0*yL>l1WK7J9lkHG>)cDbbV^^+Uk5!kuOVzmHxtOmlucp zFpdcIqiRint$Rzn%T$;;sqV6T13>`gD=8iP2P6;N%CYjkyA=+ed_xxvF1S{((l=G< zV`0yL!0^}4q>oy(Xk@Up=JCyQuv!QwSC=~qJ^D}wcGy%(e#iEf>-}Bz3$lLNr!FCH z^yg?T1}L>fO`g)9S)2_pBuLlhw3MCLJ-SL3S2tm~GPz>)@x1h{VO&qy^Ue8R4hGVW@g)Dtef{Ei?*Z9m)66Habe$&tkKc~8Jiq_7E%ALor7 zOA0tg2&OMs^LbXDnZX!6Z#n@$WktjhXbMeYVTG1!kfYDCwc>M7l22+Xu{cxj*F z*kkRv;zO$&sBRgk?W)K4`0z2TEgRV#!A|_1l&T@?5W&>syOvKC z5+B0-uOwMNsQp3mLif$S1nP-!{`Z}_;ebX}WwK8mK|n%UhL9A7?q*P6=oD!<*ZkkL z&ROexIO}|PKfQcmSTl38?|tvR?|ohO-oKy7VwRgEw9@8qu}i(%M{GtDf9Pr`C%}-yUb>GzKl5kly+Gjope!m4crD z!*Y7gW_5Nk;_Zcpl0$F@_aa$SgkfhlBn`Xfw$e3=PH>kk!sFA&8_g}m`pvkbOzY@a zWM?}1$2Fj9&L_|=Q!O^*lhYEZ-!f6S{g0=1(!(yEE|VITvhnTaCbbod8j=H;KRn>S zU!lrk=A!w&sc@jm1|d8;=DL(rmSdYz=}*tYQ&mwLb7pb6K9WQ^4CAU~e00xq!_u-o zZt0Os1GHdns5@MezkP@A*+!V|+y49OI$~qfOTODAtxeXWfdl*PWM7+$1SKMHk1+UJ z7b@o}vPj@0WQTEfqg)fuM@sf-ppZ}yE42dw%i3Y%;-sKh>cP?IyRhJh5G~z=6ywKl zb>AxS&g!YHiw?U4U(~vKs7OOT%3X4G-v>n#Hg0Tt+a^5J_$B!sTctN^)To@sHOBTy zjPrDHTaTTH=9M+$IxKBhUOxO{=Eas@jcuE41UBT}p8L z^CYs4GiK_8N743f$?6@f8PiErxc0^B6%+s9tyT5KC&K;Ja9>TggoXqAf!-LVcXaWK zW{9A}_+=!s35vQU5N?5!wcGA{?l8C8UKVVZN;kY8mc*-VGe+b7CW5g1hV)RHXY0Mi zBH5*_Gpp2gA()vw9`1VDLLffr^1P!fB8{J5;j~oCr}kGyd|FFd?aV)zL%H5TwmUOs zO;XSqH0JEL=404gG5MZ0gJ0D79vJ;vd%bc)_3$FiKB>+-h8lZytC^lbtcu_9IiF^L z1`BmeR&K#N_A6wOGfM5eR#u1uA1t>|rrU|kk9XfSo0QzyxLn_qrdR( zi7eUF;HaEJ0~KAgpT+Ja=>N3Wq+ZqDo^nnK3C|;g4lxLt@ZGyvGzh=I(aE*$nbrB1 zMH(t0B%3$uFWoU5;YvJBM8lCdVG)L4vD}@=e#O~DyptULaF|e$J={#hn;M=J(cFrj zU2k8)(A|~(#|u+h6oJxXTF zhHlGHGJy)1V;gI3oa`U(J$kRVwIAN_PJi)W?&eC9wYab+T1`bg=MTkdFQ%kFS%6-s zrgXxM>zy}|NdGF;nUwfFjfWrZX)7!J#@b`zx-Q7wmUqcymmS#?{w}_$@}xA3uA!^Z zTHLfoi|m?{WYf4M`})%~;TaA6b3G>G)|L6;B`T(3vH&tOiDOv{K1F9-n^4iBdQ6e5 zu2@Dq5H;Xuk5-NO71Ah35%K`WEF_`eB587Sw1oHFV(l?EuITQFvJ56CmIo@>45a+Stba4)DXR;=b zPbl?$l#A~NRBnB&o@zlCykndRy_9uLZPhWw_N1<5|9)rP5$&%L#umFrbOUlP-!naR zX1WezJ96yg0^IrjTwHJdXv$mxg&KNJ753rZr74tX?|uk#n1pYA$CD|1n3Dz2)F1JS zCJ3_zA*!lsr@ybI%&l+j-g1*Hd!2KeQU0D$@S28fdqf(vkb$KhznbL?#y57PUeeF? z+r57%<(Bw_5*GaZBqj60-ML*D)#}F6MXh%1Pk)oTqTx!-o&ovlI}!JO?vvk(vEp&v z5HOv5#!C4uvWB>Xw3sqPMEA5Vo)T-(JVHbfS_S>_d!a;L7GifTo<{QaC!WL`{1&c% zk*Rd4!nTd+D*HbV1L%{3&N|19w<9_P#iV(=3f}5FSr3dhSl#+EHFPtcL3+~lM6~~p z7wPU?$9rc)a($3syJr;RM8gNJuBZ&9*UFUbyRo97SlD7Onxcz5%v#8t4m_^NHht;r z>MxX{u6FJ(ou=X3)ZH3$sDQot>*gkZzuV8O!y`Z_pX5SN=-ic}BLQEHz(uZ<`#Q+l;6 z+cOnQz{MHme6=7Zj>4F)V6gP}CxeVPqu=YRII21nv&ti!4GB4dmwR?1pSej?OFGG) zd?eiJY}EHVIcK(AW>KZHCN?i|?6aA4miy|)HBIg6tIuIJTfb5Bg@f5M`V3X>EWWBy z2v+o6z4nHF&f5ekDhzwu7I9X&J;!V->RJrVO#2n(z$-K(zO)K*SLmmo6fm@ocTf9;5pW*?hQUI;vQ<)JvLe)yhi%IM+GbPrz>hyI?- zQoZ?H_fS9Ljef4OSWmt$W)27dxomjm9Z7yZR9J0VqJYd zPsr&JZAP?AxLK-e-oG zJKG+myJr8+)Omz7ojH5)#W$VFRUA5)ob5? zN+_`u%L+C8;-NPDIBm#LUWNL3V3UI6%1USBuPZYhoo(uMYK2&1J^8Wz;r9FsuKygA zwu}-f>87M*pw*RH_3KzF&Gf|dux;9Lr+Lfi-U*)&PIn2+0nf4oJizE{{NL(FKm!1R z2vrjklY~7${dDH*7D)PS^}0{ua3Zi<9)x7ZUyr;D@RmM2tFq!b7t@KfJY2Qgsn1V3 zS|8XTp<>G03nlrhEU5XJkib%ps*RdTq_*vP3A4lHRdQbmW(u;lJT6pf20FQkclP(E zM}rS?n4RKBqrYJBy_f7lea0wi-x#ba*5whS9QLN-2meFu|54Mkf z11yEZB+gROc_cCAhJV|4qs6Z9Z1EJF<**5w_cScWpsc<3xwHNa*ZRA{pYNX8^SqVN zcW`a+(Kr7{k-`ui6^ZUq%FPeLPQOY@BQoEOE9PxU#mx1DQ>t!Dpz@OYD^BiIss3vd zL9w|~JTPWxC4KQhq4x(aZ@#nandAN+cGeb;ug!3XHq_R7ni&v(3AjCd+efP~Iyktn zb@orXC|SY5GeYm6A=5o<;ZeF(Md1UPe}H-h3fAj4KJ8r9141HLm~8p zMPr|jvl*7Id4JvWz1dp7&bY9+tAvUE_X*#={fnoaT~lQ~{1?uDueZwM<3FT+o^~nv zz}83LZwiwRu?2;5%)NRQw7KPdA);xMLq|Iy&oWOsq*E=Os(3)amP)p~VF3EepN=xD z-;LrqGPB_!lYpN5=u5L0dqXHL$Pr?rqg4|Mg2A(*TeM?JNl$OIeiTrb32-Zi$?BS# zxy41j1@C7+wy^KMi$4~k?W#8{`9sX&$=dRdvLD@t$0?puS!NB_`IU36T(ad?Q~mM!n(VCgdgzW46Nhg3? zu1tn4lr9;Ssui}svXA)w{rfDacYvKh`OQ_Q#(d4Up#Aml@$lED=h54h(Y*un0xGyh<= zHdYdN?mv&iUE!7q0P4zjM8E`@E1cLq8!S*G(a(9Z6W2PKGfk72m0BsaljZFcxM1_q zyB?L0px~xNkV{I)@thFyAjBH6{jYQEd)_>d?hx)6i@^AAh_ zX}t)USN^cZ9kpS!JSs;fWpv%5Am3WDv!wD(^_x~0&E!nCY3ZPo+}mK}LzcGSl8pKE z>j8>>Wo!4SN8dW7VJoiti`IrBmG}RZaQpd{fg-3fT2yI4ktlGhn~djdaRS@!f==(p z-mK|ItFcygX!0i}WT^g*;r9>n;32I;_BL(9GTrYrrk*YWKQYs@B!2sLlyaFaLHqM} zj0mdBZn_JrYLAxPl&AZ6g!Rk(8e}9qAedT2D6f&Xz^M0%Q7$li9%bFB7K)tq9HdJg z3%u69 z^IM}bd@mzro{85xZ}Dwew~>k3UF~agahK4GS+BcEV1Ct~SkX=&@<2h7zZAsNry5aA zyZ3ojtl9)c3(@#f_{94poP~z8bI7@5;a!CtigfN%&xx|J5`(jWGeHv?emjnjs62Q7pILXpMivUkyUjpq+=VsA-0?h;Me^r^p2O zP*t3v=e^8&8yde;aw4gcblrSs^F_smM>|^vQOARq5lE)d9F=V%T%Hb|1<*}D{*%>e z-Wl*{SKEGJGV1Tyx;=I_7Jfx;#Ww0#Sazq1m2XuBqwY#RgAk}6_A|_BHPg9uBjhH} z=~lHzJRv0su~Is>Y;($h=vwcz*%x)sxW!kFw`1Mc+b%CNE%*6ZP5!g9Qce3=V79%M zA=cwY!RwSNuj^C2zdK_bErK6$urLi+JnAO3XK8`+ZE7C&G0Py<(kG;0+TC1>%}vWJ z4wdaFr#|N5E5NuuGn%wM_P=>=vH1L%U}Y|1>)`ds)cA4d=4N#CkJMzL@f-Kp=t;N{ zO7G_0?_(l2G`AG@IUnjKpuWirL=QaRRbj8|+nH_Jx%)zq+`_HJoKzZXdbzeh^O6S% zI)x=Px#C##zGQm$VDsem&qlrKmDRsaBajcc<{}L{l?Rh`pEF1k?9Ln@y=O8h|7H(z zqol|da4ZKgx4+%PA&J~c6RK>M_yS{zP z-I3hmR>A9Smh;YJmpPQ5r}EvPe7P2nA4Tz^d}D6$F`?GRW(@t@;TNA7c{prd`y=q^ zB=UtSkb{w*BeU3q|a~Dp$Yl}M;g{T^{yz;d6 z2Z5$Ojg0ADHG_$yu}M!66`GnN<3g1e8Mb3}f|6`(oy3&ugv6&{m;4FP}MjiWMd*>-c^vZhwc5AEBnNYc`Uj7mO-q~(=3Jp*o zi@?zXeF5%$%FE9#PKwjV01Xsyy>UOxOAh}s zh<0=@JP2{%sHu_Z(Coe;ku*BC<2i>4o|7q3dW0BJ0C5Zgh+i79inzk*7i4E6APlyR zn4^lyU)ZV21W~rVy*;qEKXOH>S7m@|8I(&v3kHMOnAc=X#{A;W4SBJ5Qzn`7xsTtb zp&gwLAs74f>4m1|=Kke!e~EMxyh5T34P`~eaMR0@SD{JGw*e1b-4qbE18dsP&rir}=^=>9G3n+^d_Z#rPCI6OUER%12CX1?aSznNGp#XsJ{bXjI0adcMxfd0DSg3fWvT$Z*llcm=oGY zEV{vZc(2uw7;yIL)CPh7sw+!A3n%ckY9|NkG_bP6y<-2GFIMEsaY@peF1^SWUfGYzI$ow)m4vvKhXJxeP{S})BRt-Rh z8ixXL+}8B-!c@9!U?L<`EQM$dYJK>@kA?!2Fi5%J$Q~i|OuDzY0B@(ush0e3xdH?? z-rxe(h|04-3;`9k=csnvNl!sRVHh~~MrOSH{JO!x!Nb6yZ@hRJ){^KYlH|c=iu`<) z5P^lw&0+y0%3#Oi4YaTds3C`-Ae6NAPs%-<9#B^|HUj|ZU#7t@?UOR%DeAqVvp3;b z{MEf?WOs1vCv_v?Clm^&Q3oT;PnF_0MyaS9lfw;#BQXJ4%$E9T%SOt#cE_o+sg^k;4Czt`z0v5*@OhUg(P+tmeP;pf4@w5` z50?mxIv^1^Ti#-dtwBSG?8e51SW)gUV+_4m=h3Es>32#722}tNADry=;z-8X(^JXa zk8u^#aJDCi0CCh0XgLPuc9$t+w!*w+zb@vMzTpicyAP^scX#(MJ4v9F(MxE^E8l3+ z_BZMs#woNROEXsO&JMEYk73r4-YIXb#Ig@tN67eK9=4*j4MDsfd43R?rx`*A>BLlYbciWPvWUvg~k@KZzvXu?Ko?iu$rD1cSysV&|$kh1wPa2srhI>0Zk2KRI z=hpfdi?XxtsH>^z05w_nPKJDFeQI)Y9+>0v0CNqU_CGxfi1v9Pg4Euu!7_{S@u_w9 z_WF9v|J?@_MP7Y<{R)tMhrsZ-s_E`t>do|Z(H1;;E1*Ic23Ju5pp12GyOQ6pcu9Mo z9_<0)T|zL)b6dMHKYyJOj82~a`8sp}Z}=B0{8dvQ0s;e#<)2lgr-5KXeRy?m1!&j# z%u*hTr|?QgQ&V#fD0SsU9~Vw;i9bS!-vP36Jv}vb4g~f2-~=cqUw<#gAmO6(ghw9d zrO3}uWFrD@f6j{+FY=*7tbiO~7;05Y!j9;RfMD1W2qPxH!!W%2lReI0Rvwh{D`0#4 z+M4xWMYYF+un85P6YfsP|{E= Jx@-RQe*q1ZEGGZ} literal 0 HcmV?d00001 diff --git a/pySDC/playgrounds/dedalus/mpi.py b/pySDC/playgrounds/dedalus/mpi.py new file mode 100644 index 0000000000..db01503a34 --- /dev/null +++ b/pySDC/playgrounds/dedalus/mpi.py @@ -0,0 +1,96 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +MPI space-time decomposition +""" +from mpi4py import MPI +import sys + +def initSpaceTimeCommunicators(nProcSpace=None, nProcTime=None, groupSpace=True): + + gComm = MPI.COMM_WORLD + gRank = gComm.Get_rank() + gSize = gComm.Get_size() + + def log(msg): + if gRank == 0: + print(msg) + + if (nProcTime is None) and (nProcSpace is None): + # No parallelization specified, space takes everything + nProcTime = 1 + nProcSpace = gSize + elif nProcSpace is None: + # Only time parallelization specified, space takes the rest + nProcSpace = gSize//nProcTime + elif nProcTime is None: + # Only space parallelization specified, time takes the rest + nProcTime = gSize//nProcSpace + + log("-- Starting MPI initialization ...") + log(f" - nProcTotal = {gSize}") + log(f" - nProcSpace = {nProcSpace}") + log(f" - nProcTime = {nProcTime}") + + # Check for inadequate decomposition + if (gSize != nProcSpace*nProcTime) and (gSize != 1): + if gRank == 0: + raise ValueError( + f' product of nProcSpace ({nProcSpace}) with nProcTime ({nProcTime})' + f' is not equal to the total number of MPI processes ({gSize})') + else: + sys.exit(0) + + + # Information message + if gSize == 1: + log(" - no parallelization at all") + nProcSpace = 1 + nProcSpace = 1 + else: + if nProcSpace != 1: + log(f" - space parallelization activated : {nProcSpace} mpi processes") + else: + log(" - no space parallelization") + if nProcSpace != 1: + log(f" - time parallelization activated : {nProcSpace} mpi processes") + else: + log(" - no time parallelization") + log('-- finished MPI initialization --') + + # Construction of MPI communicators + if groupSpace: + tColor = gRank % nProcSpace + tComm = gComm.Split(tColor, gRank) + gComm.Barrier() + sColor = (gRank-gRank % nProcSpace)/nProcSpace + sComm = gComm.Split(sColor, gRank) + gComm.Barrier() + else: + sColor = gRank % nProcTime + sComm = gComm.Split(sColor, gRank) + gComm.Barrier() + tColor = (gRank-gRank % nProcTime)/nProcSpace + tComm = gComm.Split(tColor, gRank) + gComm.Barrier() + + return gComm, sComm, tComm + + +if __name__ == "__main__": + + gComm, sComm, tComm = initSpaceTimeCommunicators(nProcTime=4) + + gRank, gSize = gComm.Get_rank(), gComm.Get_size() + sRank, sSize = sComm.Get_rank(), sComm.Get_size() + tRank, tSize = tComm.Get_rank(), tComm.Get_size() + + from time import sleep + sleep(0.1*gRank) + + import psutil + coreNum = psutil.Process().cpu_num() + + print(f"Global rank {gRank} ({gSize}), space rank {sRank} ({sSize})," + f" time rank {tRank} ({tRank}) running on CPU core {coreNum}") + \ No newline at end of file diff --git a/pySDC/playgrounds/dedalus/problem.py b/pySDC/playgrounds/dedalus/problem.py new file mode 100644 index 0000000000..fce924a454 --- /dev/null +++ b/pySDC/playgrounds/dedalus/problem.py @@ -0,0 +1,213 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Problem class for dedalus +""" +from pySDC.core.problem import Problem + +import numpy as np +from scipy.linalg import blas +from collections import deque + +import dedalus.public as d3 + +from dedalus.core.system import CoeffSystem +from dedalus.core.evaluator import Evaluator +from dedalus.tools.array import apply_sparse +from dedalus.core.field import Field + + +def State(cls, fields): + return [f.copy() for f in fields] + + +class DedalusProblem(Problem): + + dtype_u = State + dtype_f = State + + # Dummy class to trick Dedalus + class DedalusTimeStepper: + steps = 1 + stages = 1 + def __init__(self, solver): + self.solver = solver + + def __init__(self, problem:d3.IVP, nNodes, collUpdate=False): + + self.DedalusTimeStepper.stages = nNodes + solver = problem.build_solver(self.DedalusTimeStepper) + self.solver = solver + + # From new version + self.subproblems = [sp for sp in solver.subproblems if sp.size] + self.evaluator:Evaluator = solver.evaluator + self.F_fields = solver.F + + self.M = nNodes + + c = lambda: CoeffSystem(solver.subproblems, dtype=solver.dtype) + self.MX0, self.RHS = c(), c() + self.LX = deque([[c() for _ in range(self.M)] for _ in range(2)]) + self.F = deque([[c() for _ in range(self.M)] for _ in range(2)]) + + # Attributes + self.axpy = blas.get_blas_funcs('axpy', dtype=solver.dtype) + self.dt = None + self.firstEval = True + self.init = True + + # Instantiate M solver, needed only for collocation update + if collUpdate: + for sp in solver.subproblems: + if solver.store_expanded_matrices: + np.copyto(sp.LHS.data, sp.M_exp.data) + else: + sp.LHS = sp.M_min @ sp.pre_right + sp.M_solver = solver.matsolver(sp.LHS, solver) + self.collUpdate = collUpdate + + @property + def state(self): + return self.solver.state + + def stateCopy(self): + return [f.copy() for f in self.solver.state] + + + def computeMX0(self): + """ + Compute MX0 term used in RHS of both initStep and sweep methods + + Update the MX0 attribute of the timestepper object. + """ + MX0 = self.MX0 + state:list[Field] = self.solver.state + + # Require coefficient space + for f in state: + f.require_coeff_space() + + # Compute and store MX0 + for sp in self.subproblems: + spX = sp.gather_inputs(state) + apply_sparse(sp.M_min, spX, axis=0, out=MX0.get_subdata(sp)) + + + def updateLHS(self, dt, qI, init=False): + """Update LHS and LHS solvers for each subproblem + + Parameters + ---------- + dt : float + Time-step for the updated LHS. + qI : 2darray + QDeltaI coefficients. + init : bool, optional + Wether or not initialize the LHS_solvers attribute for each + subproblem. The default is False. + """ + # Update only if different dt + if self.dt == dt: + return + + # Attribute references + solver = self.solver + + # Update LHS and LHS solvers for each subproblems + for sp in solver.subproblems: + if self.init: + # Eventually instanciate list of solvers (ony first time step) + sp.LHS_solvers = [None] * self.M + self.init = False + for i in range(self.M): + if solver.store_expanded_matrices: + # sp.LHS.data[:] = sp.M_exp.data + k_Hii*sp.L_exp.data + np.copyto(sp.LHS.data, sp.M_exp.data) + self.axpy(a=dt*qI[i, i], x=sp.L_exp.data, y=sp.LHS.data) + else: + sp.LHS = (sp.M_min + dt*qI[i, i]*sp.L_min) # CREATES TEMPORARY + sp.LHS_solvers[i] = solver.matsolver(sp.LHS, solver) + + + def evalLX(self, LX): + """ + Evaluate LX using the current state, and store it + + Parameters + ---------- + LX : dedalus.core.system.CoeffSystem + Where to store the evaluated fields. + + Returns + ------- + None. + + """ + state:list[Field] = self.solver.state + # Require coefficient space + for f in state: + f.require_coeff_space() + + # Evaluate matrix vector product and store + for sp in self.solver.subproblems: + spX = sp.gather_inputs(self.solver.state) + apply_sparse(sp.L_min, spX, axis=0, out=LX.get_subdata(sp)) + + + def evalF(self, time, F): + """ + Evaluate the F operator from the current solver state + + Note + ---- + After evaluation, state fields are left in grid space + + Parameters + ---------- + time : float + Time of evaluation. + F : dedalus.core.system.CoeffSystem + Where to store the evaluated fields. + """ + solver = self.solver + # Evaluate non linear term on current state + t0 = solver.sim_time + solver.sim_time = time + if self.firstEval: + solver.evaluator.evaluate_scheduled( + sim_time=time, timestep=solver.dt, + iteration=0, wall_time=0) + self.firstEval = False + else: + solver.evaluator.evaluate_group('F') + # Store F evaluation + for sp in solver.subproblems: + sp.gather_outputs(solver.F, out=F.get_subdata(sp)) + # Put back initial solver simulation time + solver.sim_time = t0 + + def solveAndStoreState(self, m): + """ + Solve LHS * X = RHS using the LHS associated to a given node, + and store X into the solver state. + It uses the current RHS attribute of the object. + + Parameters + ---------- + m : int + Index of the nodes. + """ + # Attribute references + solver = self.solver + RHS = self.RHS + + for field in solver.state: + field.preset_layout('c') + + # Solve and store for each subproblem + for sp in solver.subproblems: + # Slice out valid subdata, skipping invalid components + spRHS = RHS.get_subdata(sp) + spX = sp.LHS_solvers[m].solve(spRHS) # CREATES TEMPORARY + sp.scatter_inputs(spX, solver.state) diff --git a/pySDC/playgrounds/dedalus/rayleighBenardSDC.py b/pySDC/playgrounds/dedalus/rayleighBenardSDC.py new file mode 100644 index 0000000000..ca09c5e81e --- /dev/null +++ b/pySDC/playgrounds/dedalus/rayleighBenardSDC.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Dedalus script simulating 2D horizontally-periodic Rayleigh-Benard convection +using Spectral Deferred Corrections. +This script demonstrates solving a 2D Cartesian initial value problem. It can +be ran serially or in parallel, and uses the built-in analysis framework to save +data snapshots to HDF5 files. The `plot_snapshots.py` script can be used to +produce plots from the saved data. It should take about 5 cpu-minutes to run. + +The problem is non-dimensionalized using the box height and freefall time, so +the resulting thermal diffusivity and viscosity are related to the Prandtl +and Rayleigh numbers as: + + kappa = (Rayleigh * Prandtl)**(-1/2) + nu = (Rayleigh / Prandtl)**(-1/2) + +For incompressible hydro with two boundaries, we need two tau terms for each the +velocity and buoyancy. Here we choose to use a first-order formulation, putting +one tau term each on auxiliary first-order gradient variables and the others in +the PDE, and lifting them all to the first derivative basis. This formulation puts +a tau term in the divergence constraint, as required for this geometry. + +To run using e.g. 4 processes: + $ mpiexec -n 4 python rayleighBenardSDC.py + +Then use https://github.com/DedalusProject/dedalus/blob/master/examples/ivp_2d_rayleigh_benard/plot_snapshots.py +for plotting. +""" +import numpy as np +import dedalus.public as d3 +from pySDC.playgrounds.dedalus.sdc import SpectralDeferredCorrectionIMEX +import logging +logger = logging.getLogger(__name__) + +useSDC = True +SpectralDeferredCorrectionIMEX.setParameters( + nSweeps=4, + nNodes=4, + implSweep="MIN-SR-FLEX", + explSweep="PIC") + +# Parameters +Lx, Lz = 4, 1 +Nx, Nz = 512, 128 +Rayleigh = 2e6 +Prandtl = 1 +dealias = 3/2 +stop_sim_time = 50 +timestepper = SpectralDeferredCorrectionIMEX if useSDC else d3.RK222 +timestep = 1e-2/4 +dtype = np.float64 + +# Bases +coords = d3.CartesianCoordinates('x', 'z') +dist = d3.Distributor(coords, dtype=dtype) +xbasis = d3.RealFourier(coords['x'], size=Nx, bounds=(0, Lx), dealias=dealias) +zbasis = d3.ChebyshevT(coords['z'], size=Nz, bounds=(0, Lz), dealias=dealias) + +# Fields +p = dist.Field(name='p', bases=(xbasis,zbasis)) +b = dist.Field(name='b', bases=(xbasis,zbasis)) +u = dist.VectorField(coords, name='u', bases=(xbasis,zbasis)) +tau_p = dist.Field(name='tau_p') +tau_b1 = dist.Field(name='tau_b1', bases=xbasis) +tau_b2 = dist.Field(name='tau_b2', bases=xbasis) +tau_u1 = dist.VectorField(coords, name='tau_u1', bases=xbasis) +tau_u2 = dist.VectorField(coords, name='tau_u2', bases=xbasis) + +# Substitutions +kappa = (Rayleigh * Prandtl)**(-1/2) +nu = (Rayleigh / Prandtl)**(-1/2) +x, z = dist.local_grids(xbasis, zbasis) +ex, ez = coords.unit_vector_fields(dist) +lift_basis = zbasis.derivative_basis(1) +lift = lambda A: d3.Lift(A, lift_basis, -1) +grad_u = d3.grad(u) + ez*lift(tau_u1) # First-order reduction +grad_b = d3.grad(b) + ez*lift(tau_b1) # First-order reduction + +# Problem +# First-order form: "div(f)" becomes "trace(grad_f)" +# First-order form: "lap(f)" becomes "div(grad_f)" +problem = d3.IVP([p, b, u, tau_p, tau_b1, tau_b2, tau_u1, tau_u2], namespace=locals()) +problem.add_equation("trace(grad_u) + tau_p = 0") +problem.add_equation("dt(b) - kappa*div(grad_b) + lift(tau_b2) = - u@grad(b)") +problem.add_equation("dt(u) - nu*div(grad_u) + grad(p) - b*ez + lift(tau_u2) = - u@grad(u)") +problem.add_equation("b(z=0) = Lz") +problem.add_equation("u(z=0) = 0") +problem.add_equation("b(z=Lz) = 0") +problem.add_equation("u(z=Lz) = 0") +problem.add_equation("integ(p) = 0") # Pressure gauge + +# Solver +solver = problem.build_solver(timestepper) +solver.stop_sim_time = stop_sim_time + +# Initial conditions +b.fill_random('g', seed=42, distribution='normal', scale=1e-3) # Random noise +b['g'] *= z * (Lz - z) # Damp noise at walls +b['g'] += Lz - z # Add linear background + +# Analysis +snapshots = solver.evaluator.add_file_handler('snapshots', sim_dt=0.25, max_writes=50) +snapshots.add_task(b, name='buoyancy') +snapshots.add_task(-d3.div(d3.skew(u)), name='vorticity') + +# Main loop +try: + logger.info('Starting main loop') + while solver.proceed: + solver.step(timestep) + if (solver.iteration-1) % 10 == 0: + logger.info('Iteration=%i, Time=%e, dt=%e' %(solver.iteration, solver.sim_time, timestep)) +except: + logger.error('Exception raised, triggering end of main loop.') + raise +finally: + solver.log_stats() diff --git a/pySDC/playgrounds/dedalus/sdc.py b/pySDC/playgrounds/dedalus/sdc.py new file mode 100644 index 0000000000..15895c5670 --- /dev/null +++ b/pySDC/playgrounds/dedalus/sdc.py @@ -0,0 +1,616 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +SDC time-integrator for Dedalus +""" +# Python imports +import numpy as np +from scipy.linalg import blas +from collections import deque + +# Dedalus import +from dedalus.core.system import CoeffSystem +from dedalus.tools.array import apply_sparse, csr_matvecs + +# QMat imports +from qmat.qcoeff.collocation import Collocation +from qmat import genQDeltaCoeffs + +# ----------------------------------------------------------------------------- +# Main SDC parameters +# ----------------------------------------------------------------------------- +DEFAULT = { + 'nNodes': 4, + 'nSweeps': 3, + 'nodeDistr': 'LEGENDRE', + 'quadType': 'RADAU-RIGHT', + 'implSweep': 'BE', + 'explSweep': 'FE', + 'initSweep': 'COPY', + 'initSweepQDelta': 'BE', + 'forceProl': False, + } + +PARAMS = { + ('nNodes', '-sM', '--sdcNNodes'): + dict(help='number of nodes for SDC', type=int, + default=DEFAULT['nNodes']), + ('quadType', '-sqt', '--sdcQuadType'): + dict(help='quadrature type for SDC', + default=DEFAULT['quadType']), + ('nodeDistr', '-snd', '--sdcNodeDistr'): + dict(help='node distribution for SDC', + default=DEFAULT['nodeDistr']), + ('nSweeps', '-sK', '--sdcNSweeps'): + dict(help='number of sweeps for SDC', type=int, + default=DEFAULT['nSweeps']), + ('initSweep', '-si', '--sdcInitSweep'): + dict(help='initial sweep to get initialized nodes values', + default=DEFAULT['initSweep']), + ('initSweepQDelta', '-siq', '--sdcInitSweepQDelta'): + dict(help='QDelta matrix used with initSweep=QDELTA', + default=DEFAULT['initSweepQDelta']), + ('implSweep', '-sis', '--sdcImplSweep'): + dict(help='type of QDelta matrix for implicit sweep', + default=DEFAULT['implSweep']), + ('explSweep', '-ses', '--sdcExplSweep'): + dict(help='type of QDelta matrix for explicit sweep', + default=DEFAULT['explSweep']), + ('forceProl', '-sfp', '--sdcForceProl'): + dict(help='if specified force the prolongation stage ' + '(ignored for quadType=GAUSS or RADAU-LEFT)', + action='store_true') + } + +# Printing function +def sdcInfos(nNodes, quadType, nodeDistr, nSweeps, + implSweep, explSweep, initSweep, forceProl, + **kwargs): + return f""" +-- nNodes : {nNodes} +-- quadType : {quadType} +-- nodeDistr : {nodeDistr} +-- nSweeps : {nSweeps} +-- implSweep : {implSweep} +-- explSweep : {explSweep} +-- initSweep : {initSweep} +-- forceProl : {forceProl} +""".strip() + +# ----------------------------------------------------------------------------- +# Base class implementation +# ----------------------------------------------------------------------------- +class IMEXSDCCore(object): + + # Initialize parameters with default values + nSweeps:int = DEFAULT['nSweeps'] + nodeType:str = DEFAULT['nodeDistr'] + quadType:str = DEFAULT['quadType'] + implSweep = DEFAULT['implSweep'] + explSweep = DEFAULT['explSweep'] + initSweep = DEFAULT['initSweep'] + initSweepQDelta = DEFAULT['initSweepQDelta'] + forceProl = DEFAULT['forceProl'] + + # Collocation method attributes + coll = Collocation( + nNodes=DEFAULT['nNodes'], nodeType=nodeType, quadType=quadType) + nodes, weights, Q = coll.genCoeffs() + # IMEX SDC attributes, QDelta matrices are 3D with shape (K, M, M) + QDeltaI, dtauI = genQDeltaCoeffs( + implSweep, dTau=True, nSweeps=nSweeps, + Q=Q, nodes=nodes, nNodes=DEFAULT['nNodes'], + nodeType=nodeType, quadType=quadType) + QDeltaE, dtauE = genQDeltaCoeffs( + explSweep, dTau=True, nSweeps=nSweeps, + Q=Q, nodes=nodes, nNodes=DEFAULT['nNodes'], + nodeType=nodeType, quadType=quadType) + QDelta0 = genQDeltaCoeffs( + initSweepQDelta, + Q=Q, nodes=nodes, nNodes=DEFAULT['nNodes'], + nodeType=nodeType, quadType=quadType) + + diagonal = np.all([np.diag(np.diag(qD)) == qD for qD in QDeltaI]) + diagonal *= np.all([np.diag(np.diag(qD)) == 0 for qD in QDeltaE]) + if initSweep == "QDelta": + diagonal *= np.all(np.diag(np.diag(QDelta0)) == QDelta0) + + # Default attributes, used later as instance attributes + # => should be defined in inherited class + dt = None + axpy = None + + @classmethod + def setParameters(cls, nNodes=None, nodeType=None, quadType=None, + implSweep=None, explSweep=None, initSweep=None, + initSweepQDelta=None, nSweeps=None, forceProl=None): + + # Get non-changing parameters + keepNNodes = nNodes is None + keepNodeDistr = nodeType is None + keepQuadType = quadType is None + keepImplSweep = implSweep is None + keepExplSweep = explSweep is None + keepNSweeps = nSweeps is None + keepInitSweepQDelta = initSweepQDelta is None + + # Set parameter values + nNodes = len(cls.nodes) if keepNNodes else nNodes + nodeType = cls.nodeType if keepNodeDistr else nodeType + quadType = cls.quadType if keepQuadType else quadType + implSweep = cls.implSweep if keepImplSweep else implSweep + explSweep = cls.explSweep if keepExplSweep else explSweep + nSweeps = cls.nSweeps if keepNSweeps else nSweeps + initSweepQDelta = cls.initSweepQDelta if keepInitSweepQDelta else initSweepQDelta + + # Determine updated parts + updateNode = (not keepNNodes) or (not keepNodeDistr) or (not keepQuadType) + updateQDeltaI = (not keepImplSweep) or updateNode or (not keepNSweeps) + updateQDeltaE = (not keepExplSweep) or updateNode or (not keepNSweeps) + updateQDelta0 = (not keepInitSweepQDelta) or updateNode + + # Update parameters + if updateNode: + cls.coll = Collocation( + nNodes=nNodes, nodeType=nodeType, quadType=quadType) + cls.nodes, cls.weights, cls.Q = cls.coll.genCoeffs() + cls.nodeType, cls.quadType = nodeType, quadType + if updateQDeltaI: + cls.QDeltaI, cls.dtauI = genQDeltaCoeffs( + implSweep, dTau=True, nSweeps=nSweeps, + Q=cls.Q, nodes=cls.nodes, nNodes=nNodes, + nodeType=nodeType, quadType=quadType) + cls.implSweep = implSweep + if updateQDeltaE: + cls.QDeltaE, cls.dtauE = genQDeltaCoeffs( + explSweep, dTau=True, nSweeps=nSweeps, + Q=cls.Q, nodes=cls.nodes, nNodes=nNodes, + nodeType=nodeType, quadType=quadType) + cls.explSweep = explSweep + if updateQDelta0: + cls.QDelta0 = genQDeltaCoeffs( + initSweepQDelta, + Q=cls.Q, nodes=cls.nodes, nNodes=nNodes, + nodeType=nodeType, quadType=quadType) + + # Eventually update additional parameters + if forceProl is not None: cls.forceProl = forceProl + if initSweep is not None: cls.initSweep = initSweep + if not keepNSweeps: + cls.nSweeps = nSweeps + + diagonal = np.all([np.diag(np.diag(qD)) == qD for qD in cls.QDeltaI]) + diagonal *= np.all([np.diag(np.diag(qD)) == 0 for qD in cls.QDeltaE]) + if cls.initSweep == "QDELTA": + diagonal *= np.all(np.diag(np.diag(cls.QDelta0)) == cls.QDelta0) + cls.diagonal = diagonal + + # ------------------------------------------------------------------------- + # Class properties + # ------------------------------------------------------------------------- + @classmethod + def getMaxOrder(cls): + # TODO : adapt to non-LEGENDRE node distributions + M = len(cls.nodes) + return 2*M if cls.quadType == 'GAUSS' else \ + 2*M-1 if cls.quadType.startswith('RADAU') else \ + 2*M-2 # LOBATTO + + @classmethod + def getInfos(cls): + return sdcInfos( + len(cls.nodes), cls.quadType, cls.nodeType, cls.nSweeps, + cls.implSweep, cls.explSweep, cls.initSweep, cls.forceProl) + + @classmethod + def getM(cls): + return len(cls.nodes) + + @classmethod + def rightIsNode(cls): + return np.isclose(cls.nodes[-1], 1.0) + + @classmethod + def leftIsNode(cls): + return np.isclose(cls.nodes[0], 0.0) + + @classmethod + def doProlongation(cls): + return not cls.rightIsNode or cls.forceProl + + +# ----------------------------------------------------------------------------- +# Dedalus based IMEX timeintegrator class +# ----------------------------------------------------------------------------- +class SpectralDeferredCorrectionIMEX(IMEXSDCCore): + + steps = 1 + + # ------------------------------------------------------------------------- + # Instance methods + # ------------------------------------------------------------------------- + def __init__(self, solver): + + # Store class attributes as instance attributes + self.infos = self.getInfos() + + # Store solver as attribute + self.solver = solver + self.subproblems = [sp for sp in solver.subproblems if sp.size] + self.stages = self.M # need this for solver.log_stats() + + # Create coefficient systems for steps history + c = lambda: CoeffSystem(solver.subproblems, dtype=solver.dtype) + self.MX0, self.RHS = c(), c() + self.LX = deque([[c() for _ in range(self.M)] for _ in range(2)]) + self.F = deque([[c() for _ in range(self.M)] for _ in range(2)]) + + if not self.leftIsNode and self.initSweep == "QDelta": + self.F0, self.LX0 = c(), c() + + # Attributes + self.axpy = blas.get_blas_funcs('axpy', dtype=solver.dtype) + self.dt = None + self.firstEval = True + + @property + def M(self): + return len(self.nodes) + + @property + def rightIsNode(self): + return np.allclose(self.nodes[-1], 1.0) + + @property + def leftIsNode(self): + return np.allclose(self.nodes[0], 0.0) + + @property + def doProlongation(self): + return not self.rightIsNode or self.forceProl + + def _computeMX0(self, MX0): + """ + Compute MX0 term used in RHS of both initStep and sweep methods + + Update the MX0 attribute of the timestepper object. + """ + state = self.solver.state + # Assert coefficient space + self._requireStateCoeffSpace(state) + # Compute and store MX0 + MX0.data.fill(0) + for sp in self.subproblems: + spX = sp.gather_inputs(state) + apply_sparse(sp.M_min, spX, axis=0, out=MX0.get_subdata(sp)) + + def _updateLHS(self, dt, init=False): + """Update LHS and LHS solvers for each subproblem + + Parameters + ---------- + dt : float + Time-step for the updated LHS. + init : bool, optional + Wether or not initialize the LHS_solvers attribute for each + subproblem. The default is False. + """ + # Attribute references + qI = self.QDeltaI + solver = self.solver + + # Update LHS and LHS solvers for each subproblems + for sp in solver.subproblems: + if init: + # Eventually instantiate list of solver (ony first time step) + sp.LHS_solvers = [[None for _ in range(self.M)] for _ in range(self.nSweeps)] + for k in range(self.nSweeps): + for m in range(self.M): + if solver.store_expanded_matrices: + raise NotImplementedError("code correction required") + np.copyto(sp.LHS.data, sp.M_exp.data) + self.axpy(a=dt*qI[k, m, m], x=sp.L_exp.data, y=sp.LHS.data) + else: + sp.LHS = (sp.M_min + dt*qI[k, m, m]*sp.L_min) + sp.LHS_solvers[k][m] = solver.matsolver(sp.LHS, solver) + if self.initSweep == "QDELTA": + raise NotImplementedError() + + def _evalLX(self, LX): + """ + Evaluate LX using the current state, and store it + + Parameters + ---------- + LX : dedalus.core.system.CoeffSystem + Where to store the evaluated fields. + + Returns + ------- + None. + + """ + # Attribute references + solver = self.solver + # Assert coefficient spacec + self._requireStateCoeffSpace(solver.state) + # Evaluate matrix vector product and store + for sp in solver.subproblems: + spX = sp.gather_inputs(solver.state) + apply_sparse(sp.L_min, spX, axis=0, out=LX.get_subdata(sp)) + + def _evalF(self, F, time, dt, wall_time): + """ + Evaluate the F operator from the current solver state + + Note + ---- + After evaluation, state fields are left in grid space + + Parameters + ---------- + time : float + Time of evaluation. + F : dedalus.core.system.CoeffSystem + Where to store the evaluated fields. + dt : float + Current time step. + wall_time : float + Current wall time. + """ + + solver = self.solver + # Evaluate non linear term on current state + t0 = solver.sim_time + solver.sim_time = time + if self.firstEval: + solver.evaluator.evaluate_scheduled( + wall_time=wall_time, timestep=dt, sim_time=time, + iteration=solver.iteration) + self.firstEval = False + else: + solver.evaluator.evaluate_group('F') + # Store F evaluation + for sp in solver.subproblems: + sp.gather_outputs(solver.F, out=F.get_subdata(sp)) + # Put back initial solver simulation time + solver.sim_time = t0 + + def _solveAndStoreState(self, k, m): + """ + Solve LHS * X = RHS using the LHS associated to a given node, + and store X into the solver state. + It uses the current RHS attribute of the object. + + Parameters + ---------- + k : int + Sweep index (0 for the first sweep). + m : int + Index of the nodes. + """ + # Attribute references + solver = self.solver + RHS = self.RHS + + self._presetStateCoeffSpace(solver.state) + + # Solve and store for each subproblem + for sp in solver.subproblems: + # Slice out valid subdata, skipping invalid components + spRHS = RHS.get_subdata(sp) + spX = sp.LHS_solvers[k][m].solve(spRHS) # CREATES TEMPORARY + sp.scatter_inputs(spX, solver.state) + + def _requireStateCoeffSpace(self, state): + """Transform current state fields in coefficient space. + If already in coefficient space, doesn't do anything.""" + for field in state: + field.require_coeff_space() + + def _presetStateCoeffSpace(self, state): + """Allow to write fields in coefficient space into current state + fields, without transforming current state in coefficient space.""" + for field in state: + field.preset_layout('c') + + def _initSweep(self): + """ + Initialize node terms for one given time-step + + Parameters + ---------- + iType : str, optional + Type of initialization, can be : + - iType="QDELTA" : use QDelta[I,E] for coarse time integration. + - iType="COPY" : just copy the values from the initial solution. + - iType="FNO" : use an FNO-ML model to predict values (incoming ...) + """ + # Attribute references + tau, qI, qE = self.nodes, self.QDeltaI, self.QDeltaE + solver = self.solver + t0, dt, wall_time = solver.sim_time, self.dt, self.wall_time + RHS, MX0, Fk, LXk = self.RHS, self.MX0, self.F[0], self.LX[0] + if not self.leftIsNode and self.initSweep == "QDELTA": + F0, LX0 = self.F0, self.LX0 + axpy = self.axpy + + if self.initSweep == 'QDELTA': + + # Eventual initial field evaluation + if not self.leftIsNode: + if np.any(self.dtauE): + self._evalF(F0, t0, dt, wall_time) + if np.any(self.dtauI): + self._evalLX(LX0) + + # Loop on all quadrature nodes + for m in range(self.M): + + # Build RHS + if RHS.data.size: + + # Initialize with MX0 term + np.copyto(RHS.data, MX0.data) + + # Add eventual initial field evaluation + if not self.leftIsNode: + if self.dtauE[m]: + axpy(a=dt*self.dtauE[m], x=F0.data, y=RHS.data) + if self.dtauI[m]: + axpy(a=-dt*self.dtauI[m], x=LX0.data, y=RHS.data) + + # Add F and LX terms (already computed) + for i in range(m): + axpy(a=dt*qE[m, i], x=Fk[i].data, y=RHS.data) + axpy(a=-dt*qI[m, i], x=LXk[i].data, y=RHS.data) + + # Solve system and store node solution in solver state + self._solveAndStoreState(m) + + # Evaluate and store LX with current state + self._evalLX(LXk[m]) + + # Evaluate and store F(X, t) with current state + self._evalF(Fk[m], t0+dt*tau[m], dt, wall_time) + + elif self.initSweep == 'COPY': + self._evalLX(LXk[0]) + self._evalF(Fk[0], t0, dt, wall_time) + for m in range(1, self.M): + np.copyto(LXk[m].data, LXk[0].data) + np.copyto(Fk[m].data, Fk[0].data) + + else: + raise NotImplementedError(f'initSweep={self.initSweep}') + + def _sweep(self, k): + """Perform a sweep for the current time-step""" + # Attribute references + tau, qI, qE, q = self.nodes, self.QDeltaI, self.QDeltaE, self.Q + solver = self.solver + t0, dt, wall_time = solver.sim_time, self.dt, self.wall_time + RHS, MX0 = self.RHS, self.MX0 + Fk, LXk, Fk1, LXk1 = self.F[0], self.LX[0], self.F[1], self.LX[1] + axpy = self.axpy + + # Loop on all quadrature nodes + for m in range(self.M): + + # Build RHS + if RHS.data.size: + + # Initialize with MX0 term + np.copyto(RHS.data, MX0.data) + + # Add quadrature terms + for i in range(self.M): + axpy(a=dt*q[m, i], x=Fk[i].data, y=RHS.data) + axpy(a=-dt*q[m, i], x=LXk[i].data, y=RHS.data) + + if not self.diagonal: + # Add F and LX terms from iteration k+1 and previous nodes + for i in range(m): + axpy(a=dt*qE[k, m, i], x=Fk1[i].data, y=RHS.data) + axpy(a=-dt*qI[k, m, i], x=LXk1[i].data, y=RHS.data) + # Add F and LX terms from iteration k and previous nodes + for i in range(m): + axpy(a=-dt*qE[k, m, i], x=Fk[i].data, y=RHS.data) + axpy(a=dt*qI[k, m, i], x=LXk[i].data, y=RHS.data) + + # Add LX terms from iteration k+1 and current nodes + axpy(a=dt*qI[k, m, m], x=LXk[m].data, y=RHS.data) + + # Solve system and store node solution in solver state + self._solveAndStoreState(k, m) + + # Avoid non necessary RHS evaluations work + if not self.forceProl and k == self.nSweeps-1: + if self.diagonal: + continue + elif m == self.M-1: + continue + + # Evaluate and store LX with current state + self._evalLX(LXk1[m]) + # Evaluate and store F(X, t) with current state + self._evalF(Fk1[m], t0+dt*tau[m], dt, wall_time) + + # Inverse position for iterate k and k+1 in storage + # ie making the new evaluation the old for next iteration + self.F.rotate() + self.LX.rotate() + + def _prolongation(self): + """Compute prolongation (needed if last node != 1)""" + # Attribute references + solver = self.solver + w, dt = self.weights, self.dt + RHS, MX0, Fk, LXk = self.RHS, self.MX0, self.F[0], self.LX[0] + axpy = self.axpy + + # Build RHS + if RHS.data.size: + # Initialize with MX0 term + np.copyto(RHS.data, MX0.data) + # Add quadrature terms + for i in range(self.M): + axpy(a=dt*w[i], x=Fk[i].data, y=RHS.data) + axpy(a=-dt*w[i], x=LXk[i].data, y=RHS.data) + + self._presetStateCoeffSpace(solver.state) + + # Solve and store for each subproblem + for sp in solver.subproblems: + # Slice out valid subdata, skipping invalid components + spRHS = RHS.get_subdata(sp)[:sp.LHS.shape[0]] + # Solve using LHS of the node + spX = sp.M_solver.solve(spRHS) + # Make output buffer including invalid components for scatter + spX2 = np.zeros( + (sp.pre_right.shape[0], len(sp.subsystems)), + dtype=spX.dtype) + # Store X to state_fields + csr_matvecs(sp.pre_right, spX, spX2) + sp.scatter(spX2, solver.state) + + def step(self, dt, wall_time): + """ + Compute the next time-step solution using the time-stepper method, + and modify to state field of the solver + + Note + ---- + State fields should be left in grid space after at the end of the step. + + Parameters + ---------- + dt : float + Length of the current time-step. + wall_time : float + Current wall time for the simulation. + """ + self.wall_time = wall_time + + # Initialize and/or update LHS terms, depending on dt + if dt != self.dt: + self._updateLHS(dt, init=self.dt is None) + self.dt = dt + + # Compute MX0 for the whole time step + self._computeMX0(self.MX0) + + # Initialize node values + self._initSweep() + + # Performs sweeps + for k in range(self.nSweeps): + self._sweep(k) + + # Compute prolongation if needed + if self.doProlongation: + self._prolongation() + + # Update simulation time and reset evaluation tag + self.solver.sim_time += dt + self.firstEval = True diff --git a/pySDC/playgrounds/dedalus/sweeper.py b/pySDC/playgrounds/dedalus/sweeper.py new file mode 100644 index 0000000000..13c41a337e --- /dev/null +++ b/pySDC/playgrounds/dedalus/sweeper.py @@ -0,0 +1,161 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Sweeper class for dedalus +""" +import numpy as np + +from problem import DedalusProblem +from pySDC.core.sweeper import Sweeper + + +class DedalusSweeperIMEX(Sweeper): + + def __init__(self, params): + if 'QI' not in params: params['QI'] = 'IE' + if 'QE' not in params: params['QE'] = 'EE' + # call parent's initialization routine + super().__init__(params) + # IMEX integration matrices + self.QI = self.get_Qdelta_implicit(qd_type=self.params.QI) + self.QE = self.get_Qdelta_explicit(qd_type=self.params.QE) + + def predict(self): + """Copy for now ...""" + + L = self.level + t0, dt = L.time, L.dt + qI = self.QI[1:, 1:] + + P:DedalusProblem = L.prob + assert type(P) == DedalusProblem + assert self.coll.num_nodes == P.M + + P.firstEval = True + P.solver.sim_time = t0 + P.solver.dt = dt + + # for f0, f in zip(L.u[0], P.solver.state): + # np.copyto(f.data, f0.data) + + # for f in P.solver.state: + # f.change_scales(f.domain.dealias) + # P.solver.evaluator.require_grid_space(P.solver.state) + # P.solver.evaluator.require_coeff_space(P.solver.state) + + P.updateLHS(dt, qI) + P.computeMX0() + + Fk, LXk = P.F[0], P.LX[0] + P.evalLX(LXk[0]) + P.evalF(t0, Fk[0]) + for m in range(1, P.M): + np.copyto(LXk[m].data, LXk[0].data) + np.copyto(Fk[m].data, Fk[0].data) + + # additional stuff for pySDC, which stores solution at each nodes + for m in range(1, self.coll.num_nodes + 1): + L.u[m] = P.stateCopy() + + # indicate that this level is now ready for sweeps + L.status.unlocked = True + L.status.updated = True + + + def update_nodes(self): + """ + TODO + """ + # get current level and problem description + L = self.level + # only if the level has been touched before + assert L.status.unlocked + + t0, dt = L.time, L.dt + tau, qI, qE, q = self.coll.nodes, self.QI[1:, 1:], self.QE[1:, 1:], self.coll.Qmat[1:, 1:] + + P:DedalusProblem = L.prob + assert type(P) == DedalusProblem + + # get number of collocation nodes for easier access + M = self.coll.num_nodes + assert M == P.M + + # Attribute references + RHS, MX0 = P.RHS, P.MX0 + Fk, LXk, Fk1, LXk1 = P.F[0], P.LX[0], P.F[1], P.LX[1] + axpy = P.axpy + + # Loop on all quadrature nodes + for m in range(M): + + # Initialize with MX0 term + np.copyto(RHS.data, MX0.data) + + # Add quadrature terms + for j in range(M): + axpy(a=dt*q[m, j], x=Fk[j].data, y=RHS.data) + axpy(a=-dt*q[m, j], x=LXk[j].data, y=RHS.data) + + # Add F and LX terms from iteration k+1 + for j in range(m): + axpy(a=dt*qE[m, j], x=Fk1[j].data, y=RHS.data) + axpy(a=-dt*qI[m, j], x=LXk1[j].data, y=RHS.data) + + # Add F and LX terms from iteration k + for j in range(m): + axpy(a=-dt*qE[m, j], x=Fk[j].data, y=RHS.data) + axpy(a=dt*qI[m, j], x=LXk[j].data, y=RHS.data) + axpy(a=dt*qI[m, m], x=LXk[m].data, y=RHS.data) + + # Solve system and store node solution in solver state + P.solveAndStoreState(m) + + # Evaluate and store LX with current state + P.evalLX(LXk1[m]) + # Evaluate and store F(X, t) with current state + P.evalF(t0+dt*tau[m], Fk1[m]) + + # Update u for pySDC + L.u[m+1] = P.stateCopy() + + # Inverse position for iterate k and k+1 in storage + # ie making the new evaluation the old for next iteration + P.F.rotate() + P.LX.rotate() + + # indicate presence of new values at this level + L.status.updated = True + + + def compute_end_point(self): + """ + Compute u at the right point of the interval + + The value uend computed here is a full evaluation of the Picard formulation unless do_full_update==False + + Returns: + None + """ + + # get current level and problem description + L = self.level + t0, dt = L.time, L.dt + P:DedalusProblem = L.prob + + # check if Mth node is equal to right point and do_coll_update is false, perform a simple copy + if self.coll.right_is_node and not self.params.do_coll_update: + # a copy is sufficient + L.uend = L.u[-1] + else: + raise NotImplementedError() + # start with u0 and add integral over the full interval (using coll.weights) + L.uend = P.dtype_u(L.u[0]) + for m in range(self.coll.num_nodes): + L.uend += L.dt * self.coll.weights[m] * (L.f[m + 1].impl + L.f[m + 1].expl) + # add up tau correction of the full interval (last entry) + if L.tau[-1] is not None: + L.uend += L.tau[-1] + + P.solver.sim_time = t0 + dt + P.firstEval = True