From cec5f48ece82f6a769030705aa7272007a8b9101 Mon Sep 17 00:00:00 2001
From: Nicola-Fonzi <60700515+Nicola-Fonzi@users.noreply.github.com>
Date: Mon, 7 Dec 2020 09:53:07 +0100
Subject: [PATCH 01/25] Update tutorials.yml
Added the unsteady fsi entry
---
_data/tutorials.yml | 1 +
1 file changed, 1 insertion(+)
diff --git a/_data/tutorials.yml b/_data/tutorials.yml
index e0f7787d..85806de3 100644
--- a/_data/tutorials.yml
+++ b/_data/tutorials.yml
@@ -37,6 +37,7 @@
- title: Multiphysics
tutorials:
- Static_FSI
+ - Unsteady_FSI_Python
- title: Design Features
tutorials:
From c250bd55b2171a527461a0a041bd0621e23ad012 Mon Sep 17 00:00:00 2001
From: Nicola-Fonzi <60700515+Nicola-Fonzi@users.noreply.github.com>
Date: Mon, 7 Dec 2020 09:56:15 +0100
Subject: [PATCH 02/25] Initial commit unsteady_fsi
---
_tutorials/multiphysics/Unsteady_FSI_Python/unsteady_fsi | 1 +
1 file changed, 1 insertion(+)
create mode 100644 _tutorials/multiphysics/Unsteady_FSI_Python/unsteady_fsi
diff --git a/_tutorials/multiphysics/Unsteady_FSI_Python/unsteady_fsi b/_tutorials/multiphysics/Unsteady_FSI_Python/unsteady_fsi
new file mode 100644
index 00000000..bf0297b0
--- /dev/null
+++ b/_tutorials/multiphysics/Unsteady_FSI_Python/unsteady_fsi
@@ -0,0 +1 @@
+## Test
From 2028e55fa10585137968c0e52f75457ebcb37127 Mon Sep 17 00:00:00 2001
From: Nicola-Fonzi <60700515+Nicola-Fonzi@users.noreply.github.com>
Date: Mon, 7 Dec 2020 10:02:17 +0100
Subject: [PATCH 03/25] Renamed tutorial Dynamic_FSI_Python
---
_data/tutorials.yml | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/_data/tutorials.yml b/_data/tutorials.yml
index 85806de3..e69dbc6b 100644
--- a/_data/tutorials.yml
+++ b/_data/tutorials.yml
@@ -37,7 +37,7 @@
- title: Multiphysics
tutorials:
- Static_FSI
- - Unsteady_FSI_Python
+ - Dynamic_FSI_Python
- title: Design Features
tutorials:
From 67cda58a6bf6f08629d8894eb8f4123d89780bcb Mon Sep 17 00:00:00 2001
From: Nicola-Fonzi <60700515+Nicola-Fonzi@users.noreply.github.com>
Date: Mon, 7 Dec 2020 12:12:08 +0100
Subject: [PATCH 04/25] Some modifications to .md file
---
.../Unsteady_FSI_Python/Dynamic_FSI_Python.md | 368 ++++++++++++++++++
.../Unsteady_FSI_Python/unsteady_fsi | 1 -
2 files changed, 368 insertions(+), 1 deletion(-)
create mode 100644 _tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md
delete mode 100644 _tutorials/multiphysics/Unsteady_FSI_Python/unsteady_fsi
diff --git a/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md b/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md
new file mode 100644
index 00000000..c3ad1452
--- /dev/null
+++ b/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md
@@ -0,0 +1,368 @@
+---
+title: Dynamic Fluid-Structure Interaction (FSI) using the Python wrapper and an external structural solver
+permalink: /tutorials/Dynamic_FSI_Python/
+written_by: Nicola-Fonzi
+for_version: 7.0.6
+revised_by:
+revision_date:
+revised_version:
+solver: Multiphysics
+requires: SU2_CFD, PYTHON WRAPPER
+complexity: intermediate
+follows: Static_FSI
+---
+### Goals
+
+This tutorial shows how to exploit the capabilities of the Python wrapper to couple SU2 with an external strutural solver. The problem at hand consist in a NACA 0012 airfoil,
+free to pitch and plunge, with given stiffnesses, immersed in a flow with varying Mach number. The two frequencies of the modes of the structure will vary as the speed is increased, reaching
+a point when the aeroelastic system will be unstable. The classical pitch-plunge flutter will then be visible.
+
+The considered case has been presented in $$^1$$, the reader is encouraged to refer to that reference for further details on the problem.
+
+The structual solver used in this example loosely couples SU2 with the commercial code Nastran, keeping the results of general interest. Indeed, linear strctures of arbitrary
+complexity can be analysed with the same workflow.
+
+Summarizing, this document will cover:
+- Preparing a FSI analysis to be used with the python wrapper
+- Using the strctural solver for Nastran models with SU2
+- Important keywords for the fsi and solid solver config files
+
+A sketch of the problem at hand is presented below:
+
+
+
+### Resources
+
+You can find the resources for this tutorial in the folder [fsi/steady_state](https://github.com/su2code/Tutorials/tree/master/multiphysics/steady_fsi) of the [Tutorials repository](https://github.com/su2code/Tutorials). There is a [FSI config file](https://github.com/su2code/Tutorials/tree/master/multiphysics/steady_fsi/config_fsi_steady.cfg) and two sub-config files for the [flow](https://github.com/su2code/Tutorials/tree/master/multiphysics/steady_fsi/config_channel.cfg) and [structural](https://github.com/su2code/Tutorials/tree/master/multiphysics/steady_fsi/config_cantilever.cfg) subproblems.
+
+Moreover, you will need two mesh files for the [flow domain](https://github.com/su2code/Tutorials/tree/master/multiphysics/steady_fsi/mesh_channel.su2) and the [cantilever](https://github.com/su2code/Tutorials/tree/master/multiphysics/steady_fsi/mesh_cantilever.su2).
+
+### Background
+
+The solution process will follow a very similar flow as the one explained in [this](https://su2code.github.io/tutorials/Static_FSI/) tutorial,
+which you should complete before stepping to the present one. The fluid and structural solutions will be obtained separately, iterating between them until convergence.
+Here, the difference is due to the fact that the simulation is unsteady. Thus, this inner iteration between fluid and struture will be repeated at each physical time step.
+
+The aerodynamic model is based on the compressible Reynolds-averaged Navier-Stokes equations. A central JST scheme is used for the convective fluxes, and a weighted least square
+scheme is used for the gradients. The turbulence model is the SST and a CFL number of 20, for the psuedo time step, is used.
+
+The strctural model is made by a single point, positioned at the rotation axis, with two degrees of freedom, pitch and plunge.
+Inertia and mass of the airfoil are concentrated at the center of mass of the profile, at a certain distance from the rotation axis. The equations of motions are available
+analytically and reads:
+
+
+
+The structural solver, instead of integrating
+the equations of motions for this point, which are available analytically, is intended to solve a general strctural problem. For this reason, a preprocessing step in Nastran
+will be performed, computing the mode shapes and modal frequencies of the model. Then, the structural solver will integrate a set of ODEs for the modes of the strcture.
+
+In this particular case, it may look excessively complicated, but it allows to solve an arbitrary aeroelastic problem with the same scheme.
+
+#### Mesh Description
+
+The fluid domain is discretised with 133k nodes, with refining close to the airfoil surface, in order to correctly represent the turbulent boundary layer. The first cell
+is placed at a height of $y+\approx 1$. A close up view of the mesh is pictured below:
+
+
+
+As far as the strctural mesh is concerned, this is a finite element mesh prepared for the commercial code Nastram. In the context of this example, as only a 2D problem, with
+only two degrees of freedom, is considered, the mesh is extremely simple; a set of rigid elements that connect all the nodes to the master node, positioned on the rotation axis.
+The master node only has two degrees of freedom: pitch and plunge.
+
+Another slave node is added, in the position of the center of mass, that will house the mass and inertia of the airfoil.
+
+It should always be recalled that, if interpolation is needed
+between the strctural and fluid meshes, RBF will be used. The limitation of this linear interpolation is due to the fact that, if a 2D problem is concerned, the strctural points
+cannot all lie on the same line. Equivalently, if a 3D problem is tackled, the points cannot lie all in the same plane. For this reason, the thickness is represented in the
+FEM mesh for Nastran as shown below:
+
+
+
+#### Configuration File Options
+
+We start the tutorial by definining the problem as a multiphysics case,
+
+```
+SOLVER = MULTIPHYSICS
+```
+
+We set the config files for each sub-problem using the command ```CONFIG_LIST```, and state that each sub-problem will use a different mesh file:
+
+```
+CONFIG_LIST = (config_channel.cfg, config_cantilever.cfg)
+MULTIZONE_MESH = NO
+```
+
+Now, we define the outer iteration strategy to solve the FSI problem. We use a Block Gauss-Seidel iteration as defined in the background section with a maximum of 40 outer iterations
+
+```
+MULTIZONE_SOLVER = BLOCK_GAUSS_SEIDEL
+OUTER_ITER = 40
+```
+
+Then, the convergence criteria is set to evaluate the averaged residual of the flow state vector (zone 0) (```AVG_BGS_RES[0]```) and the structural state vector (zone 1) (```AVG_BGS_RES[1]```) in two consecutive outer iterations, $$\mathbf{w}^{n+1}-\mathbf{w}^{n}$$ and $$\mathbf{u}^{n+1}-\mathbf{u}^{n}$$ respectively.
+
+```
+CONV_FIELD = AVG_BGS_RES[0], AVG_BGS_RES[1]
+CONV_RESIDUAL_MINVAL = -9
+```
+
+Finally, we define the coupling conditions. In this case, the interface between the marker ```flowbound``` in the flow field, and ```feabound``` in the structural field, is defined as
+
+```
+MARKER_ZONE_INTERFACE = (flowbound, feabound)
+```
+
+The last step is defining our desired output. In this tutorial, we will use the following configuration for the screen output
+
+```
+SCREEN_OUTPUT = (OUTER_ITER, AVG_BGS_RES[0], AVG_BGS_RES[1], DEFORM_MIN_VOLUME[0], DEFORM_ITER[0])
+WRT_ZONE_CONV = NO
+```
+
+where the convergence magnitudes are plotted alongside the minimum volume obtained in the deformed mesh, and the number of iterations required by the linear solver that updates the mesh deformation. For clarity, the command ```WRT_ZONE_CONV``` limits the output to the outer iterations, while the convergence of the flow and structural sub-problems is not written to screen.
+
+The convergence output will be set using
+
+```
+HISTORY_OUTPUT = ITER, BGS_RES[0], AERO_COEFF[0], BGS_RES[1]
+
+WRT_ZONE_HIST = NO
+CONV_FILENAME= history
+```
+
+where the individual components of the flow and structural outer-iteration residuals will be written, together with the convergence of the aerodynamic coefficients for the flow domain. In terms of result output, we use
+
+```
+OUTPUT_FILES = (RESTART, PARAVIEW)
+RESTART_FILENAME = restart_fsi_steady
+VOLUME_FILENAME = fsi_steady
+```
+
+where the volume files ```fsi_steady_*.vtu``` and restart files ```restart_fsi_steady_*.vtu``` will be appended the zone number.
+
+#### Applying coupling conditions to the individual domains
+
+Minor modifications are required for the flow and structural config files to be used in Fluid-Structure Interaction, as compared to a single-zone problem. As it is understood that the user will have a basic knowledge of the flow and structural solvers of SU2 before attempting this tutorial, only the specific commands for FSI will be discussed here.
+
+On the fluid domain (zone 0), it is necessary to indicate SU2 what is the boundary for which the flow loads will need to be computed and later applied to the structural domain, in [config_channel.cfg](https://github.com/su2code/Tutorials/blob/master/fsi/steady_state/config_channel.cfg). This is done using
+
+```
+MARKER_FLUID_LOAD = ( flowbound )
+```
+
+Next, given that we are using an Arbitrary Lagrangian-Eulerian (ALE) formulation for the flow domain, the conditions for mesh deformation need to be set. We start by defining
+
+```
+DEFORM_MESH = YES
+MARKER_DEFORM_MESH = ( flowbound )
+```
+
+where the deforming boundary is set to ```flowbound```. Next, the mesh problem defined in the background section is set. We define the stiffness of the flow mesh as a function of the distance to the deformable wall, where volumes that are farther away from the boundary will be more flexible (thus, more prone to deform largely)
+
+```
+DEFORM_STIFFNESS_TYPE = WALL_DISTANCE
+```
+
+Next, we set the properties of the linear solver. This is a pseudo-linear-elastic problem and, therefore, the properties of the linear elastic solver will be similar as those used on the [Linear Elasticity Tutorial](../Linear_Elasticity/). As a result, for this case we use
+```
+DEFORM_LINEAR_SOLVER = CONJUGATE_GRADIENT
+DEFORM_LINEAR_SOLVER_PREC = ILU
+DEFORM_LINEAR_SOLVER_ERROR = 1E-8
+DEFORM_LINEAR_SOLVER_ITER = 1000
+DEFORM_CONSOLE_OUTPUT = NO
+```
+
+As the structural side uses a Lagrangian formulation, and therefore no mesh deformation is carried out, only the wet boundary that receives the flow loads needs to be specified in [config_cantilever.cfg](https://github.com/su2code/Tutorials/blob/master/fsi/steady_state/config_cantilever.cfg)
+
+```
+MARKER_FLUID_LOAD = ( feabound )
+```
+
+### Running SU2
+
+Follow the links provided to download the [FSI config file](https://github.com/su2code/Tutorials/blob/master/fsi/steady_state/config_fsi_steady.cfg) and the [flow](https://github.com/su2code/Tutorials/blob/master/fsi/steady_state/config_channel.cfg) and [structural](https://github.com/su2code/Tutorials/blob/master/fsi/steady_state/config_cantilever.cfg) sub-config files.
+
+Also, you will need two files for the [channel mesh](https://github.com/su2code/Tutorials/blob/master/fsi/steady_state/mesh_channel.su2) and the [cantilever mesh](https://github.com/su2code/Tutorials/blob/master/fsi/steady_state/mesh_cantilever.su2). Please note that the latter is different from the structural mechanics tutorial due to a different definition of the boundary conditions.
+
+Execute the code with the standard command and using the multi-zone config file
+
+```
+$ SU2_CFD config_fsi_steady.cfg
+```
+
+which will show the following convergence history:
+
+```
++----------------------------------------------------------------+
+| Multizone Summary |
++----------------------------------------------------------------+
+| Outer_Iter| avg[bgs][0]| avg[bgs][1]|MinVolume[0]|DeformIter[0|
++----------------------------------------------------------------+
+| 0| -0.298306| -2.048554| 8.8212e-10| 0|
+| 1| -1.221569| -2.661506| 8.7196e-10| 44|
+| 2| -2.176547| -3.088096| 8.7369e-10| 44|
+| 3| -2.822207| -3.526246| 8.7331e-10| 44|
+| 4| -3.479962| -3.962556| 8.7339e-10| 44|
+| 5| -4.134404| -4.399315| 8.7337e-10| 44|
+| 6| -4.789556| -4.835974| 8.7338e-10| 44|
+| 7| -5.444539| -5.272643| 8.7338e-10| 44|
+| 8| -6.099727| -5.709490| 8.7338e-10| 44|
+| 9| -6.756832| -6.147703| 8.7338e-10| 44|
+| 10| -7.393475| -6.578762| 8.7338e-10| 44|
+| 11| -8.032563| -7.003230| 8.7338e-10| 44|
+| 12| -8.670813| -7.426572| 8.7338e-10| 44|
+| 13| -9.309448| -7.851164| 8.7338e-10| 44|
+| 14| -9.949022| -8.277094| 8.7338e-10| 44|
+| 15| -10.588983| -8.703770| 8.7338e-10| 44|
+| 16| -11.228939| -9.130581| 8.7338e-10| 44|
+```
+
+The code is stopped as soon as the values of ```avg[bgs][0]``` and ```avg[bgs][1]``` are below the convergence criteria set in the config file.
+
+The displacement field on the structural domain and the velocity field on the flow domain obtained in ```fsi_steady_1.vtu```_and ```fsi_steady_0.vtu``` respectively are shown below:
+
+
+
+
+
+#### Relaxing the computation
+
+In order to increase the speed of convergence, it is normally a good option to apply a relaxation
+
+
+to the sub-iterations of the FSI problem. This is done by adding the following commands
+
+```
+BGS_RELAXATION= FIXED_PARAMETER
+STAT_RELAX_PARAMETER= 0.8
+```
+
+to the FSI master config file. Running SU2 now would result in the following convergence history
+
+```
++----------------------------------------------------------------+
+| Multizone Summary |
++----------------------------------------------------------------+
+| Outer_Iter| avg[bgs][0]| avg[bgs][1]|MinVolume[0]|DeformIter[0|
++----------------------------------------------------------------+
+| 0| -0.298306| -2.048554| 8.8212e-10| 0|
+| 1| -1.316281| -2.740561| 8.4063e-10| 44|
+| 2| -2.215592| -3.528252| 8.6618e-10| 44|
+| 3| -3.014938| -4.639836| 8.7192e-10| 44|
+| 4| -3.732471| -5.402349| 8.7309e-10| 44|
+| 5| -4.436227| -5.790750| 8.7332e-10| 44|
+| 6| -5.137723| -6.263485| 8.7337e-10| 44|
+| 7| -5.840368| -6.711814| 8.7338e-10| 44|
+| 8| -6.536232| -7.339662| 8.7338e-10| 44|
+| 9| -7.247336| -7.555240| 8.7338e-10| 44|
+| 10| -7.946736| -8.297281| 8.7338e-10| 44|
+| 11| -8.648240| -8.586354| 8.7338e-10| 44|
+| 12| -9.353917| -9.043231| 8.7338e-10| 44|
+```
+
+where the -9 convergence parameter is reached in almost half the iterations as per the non-relaxed case.
+
+#### Printing the inner convergence
+
+Sometimes, especially when things don't work as expected, it's interesting to visualize the convergence of each of the subproblems. This can be achieved enabling
+
+```
+WRT_ZONE_CONV = YES
+```
+
+in the main config file of the multizone problem, ```config_fsi_steady.cfg```. In most cases we will not be interested on every single inner iteration of the fluid domain; we can use
+
+```
+SCREEN_WRT_FREQ_INNER = 10
+```
+
+in ```config_channel.cfg``` to limit the printout to every 10 iterations. The resulting screen output is the following:
+
+
+```
++----------------------------------------------------------------+
+| Zone 0 (Incomp. Fluid) |
++----------------------------------------------------------------+
+| Outer_Iter| Inner_Iter| rms[P]| rms[U]| rms[V]|
++----------------------------------------------------------------+
+| 0| 0| -4.799717| -19.870540| -32.000000|
+| 0| 10| -5.421646| -4.859176| -5.558343|
+| 0| 20| -6.200459| -5.584293| -6.120761|
+| 0| 30| -6.852188| -6.205093| -6.958953|
+| 0| 40| -7.502751| -6.854100| -8.034484|
+| 0| 50| -8.316190| -7.662838| -8.492639|
+| 0| 60| -9.153994| -8.509430| -9.696622|
+| 0| 70| -10.119472| -9.498203| -10.117162|
+| 0| 80| -11.035958| -10.448171| -11.143732|
+| 0| 90| -11.613994| -10.982943| -11.842820|
+| 0| 100| -12.264740| -11.616077| -12.556728|
+| 0| 104| -12.684964| -12.046580| -12.778620|
++-----------------------------------------------------------------------------+
+| Zone 1 (Structure) |
++-----------------------------------------------------------------------------+
+| Outer_Iter| Inner_Iter| rms[U]| rms[R]| rms[E]| VonMises|
++-----------------------------------------------------------------------------+
+| 0| 0| -1.121561| -2.706864| -4.419338| 2.2876e+03|
+| 0| 1| -1.708123| 1.269382| -1.716845| 2.2429e+03|
+| 0| 2| -2.583705| 0.631041| -3.175890| 2.2237e+03|
+| 0| 3| -2.508726| -0.795411| -5.686295| 2.0701e+03|
+| 0| 4| -2.751613| -1.306797| -6.550284| 2.0544e+03|
+| 0| 5| -3.147799| -1.271416| -6.914834| 2.0364e+03|
+| 0| 6| -2.962716| -2.664273| -7.800357| 2.0171e+03|
+| 0| 7| -4.083257| -2.080066| -8.553973| 2.0155e+03|
+| 0| 8| -4.522566| -4.400710| -11.005428| 2.0149e+03|
+| 0| 9| -7.252524| -5.240218| -14.876499| 2.0149e+03|
+| 0| 10| -10.834476| -10.711888| -23.634504| 2.0149e+03|
++----------------------------------------------------------------+
+| Multizone Summary |
++----------------------------------------------------------------+
+| Outer_Iter| avg[bgs][0]| avg[bgs][1]|MinVolume[0]|DeformIter[0|
++----------------------------------------------------------------+
+| 0| -0.298306| -2.048554| 8.8212e-10| 0|
+...
+
++----------------------------------------------------------------+
+| Zone 0 (Incomp. Fluid) |
++----------------------------------------------------------------+
+| Outer_Iter| Inner_Iter| rms[P]| rms[U]| rms[V]|
++----------------------------------------------------------------+
+| 16| 0| -15.808168| -15.358431| -15.769282|
+| 16| 10| -16.714759| -16.074145| -17.174991|
++-----------------------------------------------------------------------------+
+| Zone 1 (Structure) |
++-----------------------------------------------------------------------------+
+| Outer_Iter| Inner_Iter| rms[U]| rms[R]| rms[E]| VonMises|
++-----------------------------------------------------------------------------+
+| 16| 0| -11.937375| -11.637725| -25.999458| 1.8531e+03|
+| 16| 1| -15.794016| -11.636749| -28.541279| 1.8531e+03|
+| 16| 2| -16.165350| -11.621837| -28.528692| 1.8531e+03|
+| 16| 3| -16.056047| -11.608302| -28.544055| 1.8531e+03|
+| 16| 4| -16.156344| -11.627675| -28.565975| 1.8531e+03|
+| 16| 5| -16.204699| -11.633976| -28.581277| 1.8531e+03|
++----------------------------------------------------------------+
+| Multizone Summary |
++----------------------------------------------------------------+
+| Outer_Iter| avg[bgs][0]| avg[bgs][1]|MinVolume[0]|DeformIter[0|
++----------------------------------------------------------------+
+| 16| -11.228939| -9.130581| 8.7338e-10| 44|
+```
+
+where it can be observed that the multizone summary corresponds to the non-relaxed case presented first, and each of the subproblems are converged to a very high accuracy level.
+
+### References
+$$^1$$ Sanchez, R. (2018), A coupled adjoint method for optimal design in fluid-structure interaction problems with large displacements, _PhD thesis, Imperial College London_, DOI: [10.25560/58882](https://doi.org/10.25560/58882)
+
+$$^2$$ Matthies H. G. and Steindorf J. (2002), Partitioned but strongly coupled iteration schemes fornonlinear fluid–structure interaction, _Computers & Structures, 80(27–30):199–1999_.
+
+### Attribution
+
+If you are using this content for your research, please kindly cite the following reference in your derived works:
+
+Sanchez, R. _et al._ (2018), [Coupled Adjoint-Based Sensitivities in Large-Displacement Fluid-Structure Interaction using Algorithmic Differentiation](https://spiral.imperial.ac.uk/handle/10044/1/51023), _Int J Numer Meth Engng, Vol 111, Issue 7, pp 1081-1107_. DOI: [10.1002/nme.5700](https://doi.org/10.1002/nme.5700)
+
+
+This work is licensed under a Creative Commons Attribution 4.0 International License
+
+
+
diff --git a/_tutorials/multiphysics/Unsteady_FSI_Python/unsteady_fsi b/_tutorials/multiphysics/Unsteady_FSI_Python/unsteady_fsi
deleted file mode 100644
index bf0297b0..00000000
--- a/_tutorials/multiphysics/Unsteady_FSI_Python/unsteady_fsi
+++ /dev/null
@@ -1 +0,0 @@
-## Test
From 5c81c63115af2ad831daf2492cda4de8ed5adc1d Mon Sep 17 00:00:00 2001
From: Nicola-Fonzi <60700515+Nicola-Fonzi@users.noreply.github.com>
Date: Mon, 7 Dec 2020 13:08:09 +0100
Subject: [PATCH 05/25] Description of cfg files
---
.../Unsteady_FSI_Python/Dynamic_FSI_Python.md | 237 +++++++++++++-----
1 file changed, 172 insertions(+), 65 deletions(-)
diff --git a/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md b/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md
index c3ad1452..ca97c562 100644
--- a/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md
+++ b/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md
@@ -1,5 +1,5 @@
---
-title: Dynamic Fluid-Structure Interaction (FSI) using the Python wrapper and an external structural solver
+title: Dynamic Fluid-Structure Interaction (FSI) using the Python wrapper and a Nastran structural model
permalink: /tutorials/Dynamic_FSI_Python/
written_by: Nicola-Fonzi
for_version: 7.0.6
@@ -46,15 +46,60 @@ Here, the difference is due to the fact that the simulation is unsteady. Thus, t
The aerodynamic model is based on the compressible Reynolds-averaged Navier-Stokes equations. A central JST scheme is used for the convective fluxes, and a weighted least square
scheme is used for the gradients. The turbulence model is the SST and a CFL number of 20, for the psuedo time step, is used.
+Different Mach numbers will be considered, namely $M=[0.1, 0.2, 0.3, 0.357, 0.364]$. The Reynolds number is fixed at 4 millions, and the temperature is equal to 273K.
+
The strctural model is made by a single point, positioned at the rotation axis, with two degrees of freedom, pitch and plunge.
Inertia and mass of the airfoil are concentrated at the center of mass of the profile, at a certain distance from the rotation axis. The equations of motions are available
analytically and reads:
+$m\ddot{h} + S\ddot{\alpha} + C_{h}\dot{h} + K_{h}h = -L$
+$S\ddot{h} + I\ddot{\alpha} + C_{\alpha}\dot{\alpha} + K_{\alpha}\alpha = M$
+
+Where $m$ is the mass of the airfoil, $I$ the inertia around the center of mass, $S$ the static moment of inertia at the rotation axis, $C$ and $K$ the dampings and stiffnesses respectively. $L$ and $M$ are the lift and pitching up moment.
+
+These equations are usually adimensionalised to obtain results independent from the free-stream density of the flow.
+Indeed, we can define the following parameters:
+
+$\Csi=\frac{S}{mb}$, $r_{\alpha}^2=\frac{I_f}{mb^2}$, $\bar{\omega}=\frac{\omega_h}{\omega_{\alpha}}$, $\mu=\frac{m}{\pi \rho_{\inf} b^2}$
+
+Where $b$ is the semi chord of the airfoil, $\omega_h = \sqrt{\frac{K_h}{m}}$ $\omega_{\alpha} = \sqrt{\frac{K_{\alpha}}{I_f}}. If we fix them, the structure will behave always the same regardless of $\rho_{\inf}$.
+
+In this context $\Csi=0.25$, $r_{\alpha}^2=0.5$, $\bar{\omega}=0.3185$ and $\mu=100$.
+
+Note that, as we will vary the Mach number, the density will also change accordingly. Thus, with given nondimensional parameters, the inertias and stiffnesses must be
+varied accordingly.
+
+No strunctural damping is included and a time step of 1ms is used.
+
+The structural solver, instead of integrating the equations of motions for this point, which are available analytically, is intended to solve a general strctural problem.
+For this reason, a preprocessing step in Nastran will be performed, computing the mode shapes and modal frequencies of the model. Then, the structural solver will
+integrate a set of ODEs for the modes of the structure.
+
+To perform the preprocessing step, in the case control section of the Nastran model (i.e. at the very beginning of
+the file), the following lines must be added:
+
+ECHO = SORT
+DISPLACEMENT(PRINT,PUNCH)=ALL
+A real egeinvalue analysis will then be performed.
+This will produce, in the f06 file, an equivalent, ordered, model that will
+eventually be read by the python script to create the interface. Further, it will
+be created a punch file where all the mode shapes, together with their stiffness,
+are stored.
-The structural solver, instead of integrating
-the equations of motions for this point, which are available analytically, is intended to solve a general strctural problem. For this reason, a preprocessing step in Nastran
-will be performed, computing the mode shapes and modal frequencies of the model. Then, the structural solver will integrate a set of ODEs for the modes of the strcture.
+IMPORTANT: The modes should be normalised to unit mass.
+
+Further, in the Nastran model, a SET1 card must be added that contains all the
+nodes to be interfaced with aerodynamics. Note that, as of now, only one SET1 card
+is allowed. However, this should be sufficiently general not to create issues.
+
+The input and output reference systems of the interface nodes can be defined as
+local, but these reference systems must be directly defined with respect to the
+global one. Thus, if node x is referring to reference system y, y must be defined
+with respect to the global one.
+
+In the structural input file the keyword NMODES must then be defined to select which,
+of all the modes in the punch file, to be used.
In this particular case, it may look excessively complicated, but it allows to solve an arbitrary aeroelastic problem with the same scheme.
@@ -78,107 +123,169 @@ FEM mesh for Nastran as shown below:

-#### Configuration File Options
+#### Configuration File for the fluid zone
-We start the tutorial by definining the problem as a multiphysics case,
+First of all, the users should know that three configuarions file are required for this case: one for the fluid zone, one for the solid zone and one for the interface.
+The configuration file for the fluid zone is very similar to a configuration file for a simple single zone simulation. Indeed, SU2 does not know about the external strctural
+solver, it will only see the points on the aerodynamic mesh changing positions.
+
+For this reason, the solver keyword is set as:
```
-SOLVER = MULTIPHYSICS
+SOLVER = RANS
```
-We set the config files for each sub-problem using the command ```CONFIG_LIST```, and state that each sub-problem will use a different mesh file:
+A new marker is introduced, MARKER_MATCH_DEFORM_MESH. This marker is effectively
+a symmetry marker for the mesh deformation only. It may be useful in cases where
+symmetry in the mesh is required, but not in the fluid simulation. An example may
+be the simulation of a plane half-model, in wind tunnel, where the effect of boundary
+layer on the tunnel walls must be studied, but a pitch movement of the model is
+also allowed. Fluid symmetry cannot be used, but at the same time the mesh should
+move on the tunnel wall to match the deformation given by the pitch motion. However, in the context of this example, this marker is not required.
+
+The only difference from a common single zone configuration file is the addition of the following lines:
```
-CONFIG_LIST = (config_channel.cfg, config_cantilever.cfg)
-MULTIZONE_MESH = NO
+%-------------- Coupling conditions -------------------------------------------%
+DEFORM_MESH = YES
+MARKER_DEFORM_MESH = ( airfoil )
+DEFORM_STIFFNESS_TYPE = WALL_DISTANCE
+DEFORM_LINEAR_SOLVER_ITER= 200
+MARKER_FLUID_LOAD = ( airfoil )
```
-Now, we define the outer iteration strategy to solve the FSI problem. We use a Block Gauss-Seidel iteration as defined in the background section with a maximum of 40 outer iterations
+Where we selected the airfoil as our marker for coupling.
+
+The simulation is unsteady with a physical time step of 1ms:
```
-MULTIZONE_SOLVER = BLOCK_GAUSS_SEIDEL
-OUTER_ITER = 40
+TIME_DOMAIN = YES
+TIME_STEP= 1e-3
```
-Then, the convergence criteria is set to evaluate the averaged residual of the flow state vector (zone 0) (```AVG_BGS_RES[0]```) and the structural state vector (zone 1) (```AVG_BGS_RES[1]```) in two consecutive outer iterations, $$\mathbf{w}^{n+1}-\mathbf{w}^{n}$$ and $$\mathbf{u}^{n+1}-\mathbf{u}^{n}$$ respectively.
+#### Configuration File for the solid zone
-```
-CONV_FIELD = AVG_BGS_RES[0], AVG_BGS_RES[1]
-CONV_RESIDUAL_MINVAL = -9
-```
+This configuration file will be read by the structural python solver included in SU2, that will read the preprocessed Nastran model.
-Finally, we define the coupling conditions. In this case, the interface between the marker ```flowbound``` in the flow field, and ```feabound``` in the structural field, is defined as
+The solver can work in two ways:
-```
-MARKER_ZONE_INTERFACE = (flowbound, feabound)
-```
+1) It can impose the movement of a mode, with prescribed law, to provide forced
+response analysis
-The last step is defining our desired output. In this tutorial, we will use the following configuration for the screen output
+2) It can integrate in time the modal equation of motions to study the linearised
+structural deformations when the body is surrounded by the flow
-```
-SCREEN_OUTPUT = (OUTER_ITER, AVG_BGS_RES[0], AVG_BGS_RES[1], DEFORM_MIN_VOLUME[0], DEFORM_ITER[0])
-WRT_ZONE_CONV = NO
-```
+Available keyword for the config file:
-where the convergence magnitudes are plotted alongside the minimum volume obtained in the deformed mesh, and the number of iterations required by the linear solver that updates the mesh deformation. For clarity, the command ```WRT_ZONE_CONV``` limits the output to the outer iterations, while the convergence of the flow and structural sub-problems is not written to screen.
+NMODES (int): number of modes to use in the analysis -> if n modes are available in
+ the punch file, but only the first m
Date: Mon, 7 Dec 2020 15:42:06 +0100
Subject: [PATCH 06/25] Update Dynamic_FSI_Python.md
---
.../Unsteady_FSI_Python/Dynamic_FSI_Python.md | 221 ++++--------------
1 file changed, 48 insertions(+), 173 deletions(-)
diff --git a/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md b/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md
index ca97c562..fa5cde44 100644
--- a/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md
+++ b/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md
@@ -33,9 +33,9 @@ A sketch of the problem at hand is presented below:
### Resources
-You can find the resources for this tutorial in the folder [fsi/steady_state](https://github.com/su2code/Tutorials/tree/master/multiphysics/steady_fsi) of the [Tutorials repository](https://github.com/su2code/Tutorials). There is a [FSI config file](https://github.com/su2code/Tutorials/tree/master/multiphysics/steady_fsi/config_fsi_steady.cfg) and two sub-config files for the [flow](https://github.com/su2code/Tutorials/tree/master/multiphysics/steady_fsi/config_channel.cfg) and [structural](https://github.com/su2code/Tutorials/tree/master/multiphysics/steady_fsi/config_cantilever.cfg) subproblems.
+You can find the resources for this tutorial in [this folder] (https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python) of the [Tutorials repository](https://github.com/su2code/Tutorials). There is a [matlab file](https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python/FlatPlateModel.m) that can be used to produce validation data with Theodorsen theory and the [mesh file](https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python/airfoil.su2).
-Moreover, you will need two mesh files for the [flow domain](https://github.com/su2code/Tutorials/tree/master/multiphysics/steady_fsi/mesh_channel.su2) and the [cantilever](https://github.com/su2code/Tutorials/tree/master/multiphysics/steady_fsi/mesh_cantilever.su2).
+In the [main directory] (https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python), there are other 5 subdirectories containing the configuration files and structural models for the different Mach numbers. Please do not mix those files as the structural models and configurations are different at the different aerodynamic conditions.
### Background
@@ -163,6 +163,18 @@ TIME_DOMAIN = YES
TIME_STEP= 1e-3
```
+It is important to set an appropriate convergence criteria for the fluid zone. Indeed, while the structural part is linear and requires no interations, it is important that
+the fluid zone correcly converges for accurate results.
+
+The relative residual cannot be used as this is reset at each inner iteration. Thus, aver the first inner loop, it would be difficult to obtain convergence as the absolute
+residuals are already quite low. For this reason we will use:
+
+```
+CONV_CRITERIA = RESIDUAL
+CONV_FIELD= RMS_DENSITY
+CONV_RESIDUAL_MINVAL= -9.0
+```
+
#### Configuration File for the solid zone
This configuration file will be read by the structural python solver included in SU2, that will read the preprocessed Nastran model.
@@ -177,7 +189,7 @@ structural deformations when the body is surrounded by the flow
Available keyword for the config file:
-NMODES (int): number of modes to use in the analysis -> if n modes are available in
+NMODES (int): number of modes to use in the analysis. If n modes are available in
the punch file, but only the first m
This work is licensed under a Creative Commons Attribution 4.0 International License
From e8496480f8aa44f193d485fb78a6e4a940d8c6e3 Mon Sep 17 00:00:00 2001
From: Nicola-Fonzi <60700515+Nicola-Fonzi@users.noreply.github.com>
Date: Mon, 7 Dec 2020 16:20:00 +0100
Subject: [PATCH 07/25] Update Dynamic_FSI_Python.md
---
.../Unsteady_FSI_Python/Dynamic_FSI_Python.md | 24 ++++++++++---------
1 file changed, 13 insertions(+), 11 deletions(-)
diff --git a/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md b/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md
index fa5cde44..7d348c0d 100644
--- a/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md
+++ b/_tutorials/multiphysics/Unsteady_FSI_Python/Dynamic_FSI_Python.md
@@ -6,7 +6,7 @@ for_version: 7.0.6
revised_by:
revision_date:
revised_version:
-solver: Multiphysics
+solver: RANS
requires: SU2_CFD, PYTHON WRAPPER
complexity: intermediate
follows: Static_FSI
@@ -33,9 +33,9 @@ A sketch of the problem at hand is presented below:
### Resources
-You can find the resources for this tutorial in [this folder] (https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python) of the [Tutorials repository](https://github.com/su2code/Tutorials). There is a [matlab file](https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python/FlatPlateModel.m) that can be used to produce validation data with Theodorsen theory and the [mesh file](https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python/airfoil.su2).
+You can find the resources for this tutorial in [this folder](https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python) of the [Tutorials repository](https://github.com/su2code/Tutorials). There is a [matlab file](https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python/FlatPlateModel.m) that can be used to produce validation data with Theodorsen theory and the [mesh file](https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python/airfoil.su2).
-In the [main directory] (https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python), there are other 5 subdirectories containing the configuration files and structural models for the different Mach numbers. Please do not mix those files as the structural models and configurations are different at the different aerodynamic conditions.
+In the [main directory](https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python), there are other 5 subdirectories containing the configuration files and structural models for the different Mach numbers. Please do not mix those files as the structural models and configurations are different at the different aerodynamic conditions.
### Background
@@ -52,19 +52,19 @@ The strctural model is made by a single point, positioned at the rotation axis,
Inertia and mass of the airfoil are concentrated at the center of mass of the profile, at a certain distance from the rotation axis. The equations of motions are available
analytically and reads:
-$m\ddot{h} + S\ddot{\alpha} + C_{h}\dot{h} + K_{h}h = -L$
-$S\ddot{h} + I\ddot{\alpha} + C_{\alpha}\dot{\alpha} + K_{\alpha}\alpha = M$
+$$m\ddot{h} + S\ddot{\alpha} + C_{h}\dot{h} + K_{h}h = -L$$
+$$S\ddot{h} + I\ddot{\alpha} + C_{\alpha}\dot{\alpha} + K_{\alpha}\alpha = M$$
-Where $m$ is the mass of the airfoil, $I$ the inertia around the center of mass, $S$ the static moment of inertia at the rotation axis, $C$ and $K$ the dampings and stiffnesses respectively. $L$ and $M$ are the lift and pitching up moment.
+Where $$m$$ is the mass of the airfoil, $$I$$ the inertia around the center of mass, $$S$$ the static moment of inertia at the rotation axis, $$C$$ and $$K$$ the dampings and stiffnesses respectively. $$L$$ and $$M$$ are the lift and pitching up moment.
These equations are usually adimensionalised to obtain results independent from the free-stream density of the flow.
Indeed, we can define the following parameters:
-$\Csi=\frac{S}{mb}$, $r_{\alpha}^2=\frac{I_f}{mb^2}$, $\bar{\omega}=\frac{\omega_h}{\omega_{\alpha}}$, $\mu=\frac{m}{\pi \rho_{\inf} b^2}$
+$$\Csi=\frac{S}{mb}$$, $$r_{\alpha}^2=\frac{I_f}{mb^2}$$, $$\bar{\omega}=\frac{\omega_h}{\omega_{\alpha}}$$, $$\mu=\frac{m}{\pi \rho_{\inf} b^2}$$
-Where $b$ is the semi chord of the airfoil, $\omega_h = \sqrt{\frac{K_h}{m}}$ $\omega_{\alpha} = \sqrt{\frac{K_{\alpha}}{I_f}}. If we fix them, the structure will behave always the same regardless of $\rho_{\inf}$.
+Where $$b$$ is the semi chord of the airfoil, $$\omega_h = \sqrt{\frac{K_h}{m}}$$ $$\omega_{\alpha} = \sqrt{\frac{K_{\alpha}}{I_f}}$$. If we fix them, the structure will behave always the same regardless of $$\rho_{\inf}$$.
-In this context $\Csi=0.25$, $r_{\alpha}^2=0.5$, $\bar{\omega}=0.3185$ and $\mu=100$.
+In this context $$\Csi=0.25$$, $$r_{\alpha}^2=0.5$$, $$\bar{\omega}=0.3185$$ and $$\mu=100$$.
Note that, as we will vary the Mach number, the density will also change accordingly. Thus, with given nondimensional parameters, the inertias and stiffnesses must be
varied accordingly.
@@ -106,7 +106,7 @@ In this particular case, it may look excessively complicated, but it allows to s
#### Mesh Description
The fluid domain is discretised with 133k nodes, with refining close to the airfoil surface, in order to correctly represent the turbulent boundary layer. The first cell
-is placed at a height of $y+\approx 1$. A close up view of the mesh is pictured below:
+is placed at a height of $$y+\approx 1$$. A close up view of the mesh is pictured below:

@@ -317,7 +317,9 @@ TIME_MARCHING = YES
```
### Running SU2
-Follow the links provided to download the [main directory] (https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python).
+Follow the links provided to download the [main directory](https://github.com/su2code/Tutorials/tree/master/multiphysics/unsteady_fsi_python).
+
+The preprocessing step in Nastran has already been performed, thus you will directly find the required structural model.
Copy the mesh file in each subdirectory, then run the following command from inside each subdirectory:
From 8926caca57bd9ac18ff5db62c3d3ad049707f45b Mon Sep 17 00:00:00 2001
From: Nicola-Fonzi <60700515+Nicola-Fonzi@users.noreply.github.com>
Date: Mon, 7 Dec 2020 17:13:51 +0100
Subject: [PATCH 08/25] Create dummy
---
tutorials_files/multiphysics/unsteady_fsi_python/images/dummy | 1 +
1 file changed, 1 insertion(+)
create mode 100644 tutorials_files/multiphysics/unsteady_fsi_python/images/dummy
diff --git a/tutorials_files/multiphysics/unsteady_fsi_python/images/dummy b/tutorials_files/multiphysics/unsteady_fsi_python/images/dummy
new file mode 100644
index 00000000..8b137891
--- /dev/null
+++ b/tutorials_files/multiphysics/unsteady_fsi_python/images/dummy
@@ -0,0 +1 @@
+
From 48ab04b6b8cff715e3a62caadeeb6b71f6f401dc Mon Sep 17 00:00:00 2001
From: Nicola-Fonzi <60700515+Nicola-Fonzi@users.noreply.github.com>
Date: Mon, 7 Dec 2020 17:14:23 +0100
Subject: [PATCH 09/25] Add files via upload
---
.../unsteady_fsi_python/images/Setup.png | Bin 0 -> 188315 bytes
1 file changed, 0 insertions(+), 0 deletions(-)
create mode 100644 tutorials_files/multiphysics/unsteady_fsi_python/images/Setup.png
diff --git a/tutorials_files/multiphysics/unsteady_fsi_python/images/Setup.png b/tutorials_files/multiphysics/unsteady_fsi_python/images/Setup.png
new file mode 100644
index 0000000000000000000000000000000000000000..842e82ef6aae9038612a9127ebe9972520941f7b
GIT binary patch
literal 188315
zcmd3NWmp{9wl*~GlHkFDySrQP2<{Tx3GR}F5F|iwx1d3qH0}`GT^kGT?)nvT&b?^-HUdu~KQM`6?Ft@TbgMpC?iAzAzQtiPH*oard)u8~XNIGQ8VNlR~#Fr5NPQgft
zLLEzp`#zLZU!%1$j9$YF=-9(Z>6j)AAU_5?^)3U(?5Q5ew$52ff@uhg9q>w0w
zRNloWClv8Z4M34xaF=VMNMI<&E`2=lVP;ug4>CSAl@OYZq^MU2!exLFM2v!if7^uZ
z1X{)M^i)+B#)80&CrB(HB&8k24HP%e3=u+xHY_lH2$lv)Vt%DJt-U3Bu&};&@99&=
z95hB^Qf7|)tlEL)?7?MWofB#mRi8^V!tq^C%8w4``w1qD0v;0;?!>Rou=G=$dGP`YJ7o!{PeWDjz#-Z0)07eF-C4RT4<}Mx@)9
zfYo0#>hNZr+KKzzAylS8;szYnj?xoV2xBxZDX+q9YS)bzntCZ$iD|A*z5W7gyF2b3
z`To6J~PsfBDF&4h}b*-IK!4!E&p;f3C5}~LY!s(o4r?v+
zTf?^UOE|G|_Z;;w*U4x_-w!f4F=M5zsCy;qnbRbc-_mn@B$Pm>5?rSBfY}e0ZMDEI
zw&y3gAUZt{6aB@Jm3A25_+?n4hDSz)`Jpr6T_)Tt)|hnak!yPmX~k9u(#rj?U`06B
zyC+NUy<0r_Jt$IVo+#83H>>|BzgJF`!tsf7{zBFf_g>}GoA>mXixJ~)qvy{c*}a=n
zS8r~TIbzlKUk_EEI5+o}V^@7$w0;{__F`A`vHt
zmu%m`-LHxZMf>-6P4yS{=KQSetn54cmWTE9a~?qfFQNN)+fRDYG^|GNmNdoNm*o+t
zWXtCC=JjA-`>KwE3pM)#9T0vBKgh%{!Uw!!cxvwZ^{51v|CSJ@Glh+FC&AIX#rRlb7G
zt8G3}vY&R4nw8R?dY_6H2O5GE(nyWllensMOxF>0nFCTmu_Cp?+T!2hEQ~pgD3+05Ha?VJ4Mi=yH;1}jA?AuHty
z-3m!TAU>=~uY~!wI8BTXb6&n+Fv{vs9IyvyywYvvISovk8?MJFi`AUHJuLTti)Y
zT`{-y#mH;osQN}XrTG;!z4p(#
z=e>)$#v&qYOKvqoCPX8I66F)7Axgt}!=+19NhGxzeBk{+Cgm%I)kfX6^(pR?EQU}b
z;aEhE$vS5O2NpY+OAMWhq?3r2xXMFoudJ_aff$1il@pgJj(Cr&&iu4B(d?OZs!7Ao
znuW_!6#HdU&IQxt9@-dH+W9t&Pd69@=xn5ZMwbqG1%)wIsar~8c6lXvjppm-zsxzy
zUCrqR>IRb5CDv!wiv}W-+_l-6LCnF}@3d-*Q8mztx<L~B%g<4=cB$O=-!1(_#NR$
z;+_fFweHo)UG=p)d>LpIt_t1-UI+0Qt_Va9uZ@KA!3|#I7rzrkz?*;QgZ_uL*5LPx
zPQJC}waS!QavpM4a=Hb6mTxz&hrNw8jJ1t|Ke&(d5(g&j=*iIM4BTd+du2(-pOBBx5X_L?_I2vCoQj<7oa!FR9;O%}+VEs{
z@1H?WW1h+@NHH?-Fu2NEAG8mL*wDJ=dFI9C=x3OUSzYHi$0*Zoyc89bu|W0jZFQ0K
za=P9>jl3oF!)xaIknn@@n`{eSF`n=fL=;W5dg=S@xNN$l`y?i=4D&0ClHP4xq|m;O
zdE8Rd##I?N>GSp7Y5@)bA?Iu780U+%uP37?n2Kk)0;BPT>#F&$>=Fd}6ii94nyze)
zy;r{;lcGv<$d9}>g26f{Q9xRA@4X3v4T=jo$fHpf;@5YmbMe|i&nF7A5N1GSx_(D1
zs>)vO8NXJ1RlYy8-oGG-loZ2}?^${_7vwCo5i|tl$q*E?xCeKXb})q*vdI|?*9_bK
zY6!w=MtdmDldhLGa9Z(hKHy{fNg$73$>L_%V1KlsV!U<5IR!a}L@L574nDA72Zxeh
zCz!L^>8@AWjyWwmZ5{70fi!n@3acTOYi;T)8WSbvC9L{)UIVc_8@(dEeBj(t_0qk9
z&)+SzEXwV#>LPbKcEIY0n)^_rI^&~H)g(bAY!0ZiYV|X=&NJ&+y41S%_85ygnr52n
zOjbn$Scdi`n&OPb75v|Ci`rhIRR1d{_ZZj#cKYujtPhB9ZL=-<)pgnUBz7H2N0-_#%h{B>
zr?y&Hp5OK|L%2=A`L_6BhHArSh*We))X=|mM^eAYVkgT`#9w55eL};1Nk5UJ>j4rgK0$j_2hqkj@0|{ZpysprboXG#{L|9?
z`D5?Ai|>B!x4wqKnq4C|loz_>?P4`|uMZDrs5-tyrlFCKktc}#@KCuM-`o0HT2`v)
zuXOwEBxzUM(7h+a%irE(|L}IB`Tf21UefT`mE?W)!&98BjsCK@fO;p22WekpK#l6K#4??+@fNY0OHRF<8T1k*5w0waQu
zPRiP9^HcgP%^{Ow+j|=(@S%O%`^LvX}<+54HyjBfYwO7cXg8G}+
z6HjIG>QiR%F%~?K53X+vlVHWI;g@0>nwuR={2%B%WIOC(hhxD_Hs4Wk7BH#^3k%LZ
zv0hdc9Ap`GT-u?jk$D6-K>cAk)q<<|dnG1-q$A!;OWs^r8HN!!MutIv#f5y+RiXA_%x5du<|O@
zhcGZ8S1UCw7cFHaK@$f%R%26#w`Qy!c8-ttff4o)1P<-YT#P9^>}>6w1wBNb|1m-k
zIDS0M_MGC6Aucu|&$X0aQ%E^DnNjeta0s?I8oNSz&EWijBXHR<LU??ubixal^Q>CQItcN?Z3|bb-aq1vxBYc
z<8XC*D;H7De^2?#>3`1>W_#S2KW)^%q~MRUfEb9P3bXyMNQt7tt95$9z=*@hOG~JE
z!0yhwyO2)#Zr+f49gtGUVZ(US$;ih_^L(F$*eX_;YX5i&l}{C>ySQ4QzGz%asEodm
z)!6SnR(9I+CQYyriZ|`63S}u$9e7G}{`?#SqImOp!Q20a?4e;{_0l57Ivk2-HT!woIaVR!{l2Tmde?>s(^~J==f*
z)B|C+}5r%eC^
z#{Hkxiv9nYz?2D2?V0zSXGf!BDduw8~~%>ukOX&Av$M0Q#l#k6!wO#jN@~u=
zZVGxnPb)V%CPY8q2sgTNq`gyD;%ZJ(URpn7YZC9R9AIYaa=w8PBw}JhMZ>|WgS-5e
znnQFpc^T{%!EvU!jk_}7bt15!6Y+JLQ1f;h54HNU&0po3!U$A)X|uThs~i;1@Y1-!yduy>lz)mDZhgEiznO&J?{hGkCT(n*0vPT*gGYBo|AY
z{ihq?IhhCS8o9DRhoUCL=J3lahw`YVo(ZrByc~5^P{!VIR%pO0NE!*8A^b0CdQZEw
zEv_gpA3_yVh_q)vcoH&{Dq!OM32I+&H%bw=
z;eA;wE*yY4KDFBGx%(db=;$aL7W|sQZLtzi`y@S>UnhPcg2Dz;QXi%~K@%yRJMC%|
z0phRWU%colCZiv)n;$hL>*m{E50=MaSKR#{#(J)_ug@i+4+nc3>f&D=?A!$>*Lyx!7l1W){@)
zZeGz>sYbpp<-^Ywf{Fjj`Xt5$tiZ8UR#s{_xegVLz)lt{GKREAywSNi5h-b0dOF$G{(e=dI@ptS
z{e)$SOg9$3FX76ujx4LwPi;fyv)WonM*Oj+=N+nFoLQ#RpZ
zNEOKtRoVGl080BNwzFVbk?nEa`<*#7PNN^*YOF_vICqF+u5*FDr!;1CfY2Tg65-v=5E9@P>$2nb}S_ycc7voJI
z@4`kWCR$$ELaUTDN_2S4_7g)eQal=;>{wjukHJJORGV}~Xn!IlCDjm@gA$XFFcD+-
z4h+QU+Sjfn8&Ba;XO|U&NaeOSOG8Xud$Hi1o0)0!Q=>_2f2LgaWh(!>0YqY@
z5%DmxiwRA`%GqTi+hNl6-V%El)qw8uui$(DFuO@RH>qhn&gD7Ie?o`fGx{0JNjgH7&fA>GvljO(+hJBAnf(4G&e
z=EVaLq_DEEur%==X@f<+_+7S<_ljK1QC-ZtFUB_N?v{Gii6|r|C$DbhXHv0bO&|x30b}EU@|-x~;(gwu)D*v>O#;e-=>g
zj0&I1Fg(FUk>5*PY}F8igYyg0ed@u$s+15dlQik`#{Ro>vKQ)KGZ-|PmXSdGy-MRL
zXKLZO+c#py9^bE;w1ct8<1NxXJ9*z+DGilKS&7{2mLKmq_;g^1*{2kx@BwLPu$8uc1#3D7^3^bS#?!_E=$GVnQ)AtFn^`=IQrb-kYA
zpd;lr|E`Nv6k0mz7RWTT^Ft_L-EAZw;yf=&dF0^Iz3GG$s=wjIW{0If148yd;Z8DO
zh6w0?*NeY+<$C8%r?gG8KIdkQMJ+^9%7cxP*jRv6vNlS{xdCWY{q<+l8
z%)7H`+LC{X)hVRmpOIP~O^`^Q=qkml&1GD?@P1}#BuzdI#UlkV2?^<{sH)01J6mF6
z`dWeo^n4JIPh0Csr{-9iLiGbAt0yu8!ByepHl0xD;tMA7sA!^u1a(s
zN4-qoZm2`^-VV2hua8ed+>SIsMOT=*q~9z$;qj+s5I&3x>R{WB6HJZ-V7{K=)->u^Z?W4zTG?TuumRyxK>udrvNAAZB60hb?f~fjIsIc8eY)d+08kYGh6|u)z5ez
z%Y^nnVUVn51(V@6TBOtTqsyx3=rAtcuJe_ovL<77m8S}Oj-rx%F`Q~YbF&-h?|-_-
zvr^P>_!Dh(a&o;=2kNv_Qmq3>?GGhpFQv55(9qSJ#fK(43TvO6p4852;|Lhg6l#V0
zSyF&PSgblv&y}LGiq&xuI;Q6ZVZ9q$JP`8}f=Kn@oG#T$c($(E+S+m`_+;W!TKXU~
zNt3ZDNv1cFitMv}0yuw(a}Q}AB?xSQNcS=oJMH!juifk_*D%EWzyPdj>G&&zE|~1)cU2KUTWO5
z!3WV8Ji9i-y9vV{vHKOVSp)F?PM;NI1%eqEM4kEIW=t-btgh%nSZATuXxAWuyg!
z=V+TL<7~hA;pmFBi7ant*w9Vt^}O?t$QUNX5{iVh7sK#sk|ZxHtCdNPeLqu#`Tm&g
zL5162mLgQJqN>U(&00_bre3_yv{9P7dlU-YaV;H4;d{l@&Dmt4q7rR+@aR~>HTBcn
zU0pqaKIbO@;e9K3bCr2l~^8EVoXEi2pSI;GFo#9vaY
zKSbjDIXxrez}eq$sLwJB6?;@SoXT6S`Q9Xjz>sL+?c$xMk8Q`!B>?-p{zl$(cQaj
z%XAryO&|hoD<}Vru1C-b7uZNMNLIPIYc8+Q7Q1Q~J_|&fN`vf6^G$toxjx^#g0fGk
zr-5*osPt6k7x41Vto!aXC6h++VI{O-99k&@C(|_7Ny0H^@I1L-q0|yHc?SL&1OnN`
zw9vJ1|I*Y%lI+~Y>XFx0YvmB1gpX0taDTmHP_K;FeLNe83NI-x{$wg4CFd-^k7RW`
z$curiXQu~@3*w7c6gMJ6aH-)K=*Sr>svnd~ln<7E_*#UG3dgdE-Xs~(28OawQ
zMfE$Vp#_S-sJPrV%Y;SUoB0))-AHDBg6$7+jtmdi@)RFT6`(q3I2jE8ltV05I$vqy
zXmWGDw~_*HYd;lV2$GrQufMBUU=4V!sYz!sls+*r^>n}lVmkw|+0LtH7Z$dyFVuJP
zuCKC={MIdZ;7=wSNme;fm@U=Y(b+kh?hH|A1iN(%_Sg$;V>^}yock0N^0oFsyqLaGN>gn%ITrID@L>RaU$9;+o9VPPc&e2ebcQDNW
zDTTq0zbAp%)}Vz~(WH(+#Ki@Lg+J^<&N<~@zO2$G_^jZ|ULd0GGCD_NXJV7J#L
zeDC}lCTu<$keLH|mX^xzK+ENDa5@;ddTNvD3pgP%*u5haw;wafemCbRAzwmb)y)%s
zoG~4n-W?0~^r0>s#O}Yz$^DC)Rz>0G9L^SbM!WuksQL5w_;@C~wo47Y?~*y9#!8CBi`5
z=2m2`pTC$DA~Tc~C^Gcc{mZtoyMAtNZgf=ClILeD78|3%a)4Ypd;3xgo<8In6{B|H
z)p21a^0Y=~^0XqJo?WOkZKO>7v7i+<2)RZ|+&8y3%So>pxu&%VAtWC`S;+>rqafO
zg6@;oCr!IWEUNRtNvF6#I*K!}f%`ezQ*C*9xl~|6sq{@SMmZA+pM6D9x<;CQW^P&8
zutVF5auJU0$AAn3uWQ7$>Ds=@5B(;zEu-^7n;Tq?{6KQtA_f{LS$A)2G@E=%L}dS-
zMF8SC@qDJDP8>37et{e0uhR*{Tfc>hbJR=SK)h6@8frb&XVg^2({eni>(JBE(o+2L
zYg(EeiFo2RWgVbP4k6*;;o6lZtU@!DGR~_?6VSJ`jncinz3ZdxM0-Ud*1}=)Jl5KgXsA~JP)jF6VcHr7C?l;u<@+ATQgj-r}MYrW!uXQI83+RqmdaaMx
znZ?IGV%ie8gk)rKD;wLzTVj6huTf7_O#8MaXKdfI(gYDTz~kiH{IvSrTk|2;Fb;F7
zze4p7Y(DBF$sTi2(9%leNlikOJe9c5mI*5kh{Akx1_bm_Cc{28tdTvkTcqte5%5$g
z3o-q1T#K)dfBl_H+Pv$gqh27RmxFMazOFAqn+vkQP;^ofL%OG-aP?7T@s&2KUd>YF
z>|l=#6(cq>Bw*c<5p2)V8b>VPQuo4IS;ZR_)CA8>S^E-)HedbYpfcz2=g%Rozi{p2
zORkHczK}1pY*P>3#l9sbK`9T09q4!bj8x_|%B0zY$&umOZ4Q_zUKr8xR-Sgo`i3r}QM2d9Sh6JaMY^-!eLj7cS2gSa@m2-ejs57y^uk;HNr08r`
z%TP)y+|LKblK6{b;9dB0FQ@54xKw@gD$WE%C;I_BG8%rzxuWumUj&yYt4m5vWc7%6
z$XH*&l~_I4U(Rb_+~QO7&fe5kuaiT@Ab9*NaY(c3#)HIWKJM(phXsM=RKNKO2QSd|
z$^N`AL9zTvO|juux~UR|Smh$r;)93Nh@TpNVOSR9QjNHxt83+6+eT{Yvm>9^`v)o6
zA)AIkMh9&(H#6&b@=BO;7HTZGw%id`m7;u(J7b`t&ly{axDA4kR*3ao}I5kbX&WWiFnUJwS=lITXoTr}B`0mS9
zNaN3ehd0F3qw?lSDJGVd(t_UCxkVP@f^|0o8?alWo);}$XLnDwk%|Mp^d(xB;tXem
zAw_8F#+?7u(T`@)ox<}~RejW5R-;Aq*$W*4YRNE=5-
z7rX}=YQM?OD(Sl)EV*EI@Jiqg|Jrz&3a-w}NNQls4dYt7b`4wedg9Te@Csd}h@a*6
zSV=q+;TeHt(jMG>?a!^@*8nPTE*MVhKb6kRwWU8YJxxC=)LA|)(vng01u%vC2{Y?3
zNglgNepNPMh@pAJgMId5BGq7{uXWz|^;y$@m<*I1ETFPedF`BWlgiT%ZO3!7n-EXb
zXTz&s;P+y#H{PdBZ=2d!yg#qjh=3V7Q}vYW4r?2E0;!i|JS&b=zrGH#gQX*MM-}b*iOkN@i%`m;wFM3*`3179|4>2+x&u<$1a#
zM{8nUl|2*T*qr&N#_=@`uDODw>Lkh%ZgLEoOOOdIt^OY#`{)=)yN-SywQjD=@5ib+
zEkbi(#c%EbTij>?Uhj$J(OKq#A?GZd%t@O^Jc&cYyUs>@9}F?x37+Y*vP5h0q4R5-
znNpT9BHjUeF}oNr%zH%py83#s)j5#-PC(nIM5k{g{{u_rJ25fTTDPc4(h`h7li3>gflk?R?UE+OJ2x=6
z60_ze`G%57L;^C9a(K<#Fi7V&NyZc#M3cD~@;>#c4W!;_?@d-dX4Hr
zO+AeT$0e$!1K`d3qGlt`q8D?{RZuH1k>N)5e4qlnIb~+sdn?YH#xB<6bfz<-Z
zDgDz{&zCQvSIv{;HqCK37BHL!rR)-&6(rDXc7l
zj?Wn;Jw5%Q*XwH0`vEsbQJ3AYK
zDU)7icI=bAtjYjVuF!M#*&ild75nuI-yVShC!oLpmSWMVkh5<*(*S$byKTQ{n*AOq
zig9Tr-c{akK5KIv2lnZx>q_6DWwnWt`5@s2OE%9cB;@@r-{u)M`l%_P
za;@flk#So2L^C;r_>dMv^_MOIc$Jnnf$>I&Bt8ggZu@f#Kr|Y&@MQ06mH;0kBcqQ!
z--MJ)z0K6G@7%AkH{VXhZ=F^pTCw3YmYVg&_Xr&{thk0GudH{8-T&|b660#nGXoMI
zy*EIJjGzr9oDIAq4T6zk_chZPH)if%>e9f!Mk-{(5gU8^`WgqnwFPWpTV`DKTQp35
z(Tt6Y>lu3`#*&tnwxVv}Y$!nn#wZjwA_fVj4s*nG`wfjK0LrC-h!~vzgCM2<1%)W9
z4yvV_nVz
zNK!}N0zC6`8@_iW`yKV-Xsf5P@HsF>m+au5oU}#uKSTYkJcs--fCf8tVkgmRGk==u>IhY8ciiaN4$EyX%IqpYW`s`fr~m}+~@#`Y@v
z8{9+YU@40rw}#%_s?_ZCUjCHm%$;tnQ!ZXh8B{aa=;jxeKRvYCp&ImZrO^;507=6!
zCGV`Z1j1D{I*jj;EM9&p$kP9TRnBxvsnb0ZSdJil
zbgLCqqs7L`4Q|XkcVyKl^HP?Amsce|8EK$nvP^MH{~82_AB+%ML;*xjoMh
zKj*%>@PBaM+DVxqT)n2S2MUG9K);F((ADB|-cRIEWHJ1+luj-a?Au@mz+yA&M|)ix
zO{AlY
z&0!yd*HFrK9CMZlR?G$}_HziFDo>>Ica%$USgy;9^0L}ipIgOq#>A{jx|WYmx+gx!
zhF}j#*SDT(`Rosv9N+48nYzkH%(>|oaS$?mQC;``9
zm&3Q!JY^ZDn@QXvx`lD?^53!5fEFL*k4e?&o@v8U8s3c0M17{Wld{{Cw#v+B6(z~<
zcl5|1V4X^Z3Zx*>%xb2(y;|oD=_l1Y7NiOaF|`S?v3gHpe?Z-{rnaz+>0YH*yjMRv
zJDYrMp8P~Mcxu#u(hJMQz%+W$W+65wCnqy%>(Hbtk~~2w7*l%@$L<@VI`>j+YAW;0
z%nZcu&WnkVFWr!q9%C)1^eIr>kE61%UTcISSiP(jT)0hd)mM5M{uNJV2I6BCpJbtp
z(nv-1)4~y+R|Y0f9BmW)EdhYpZ*ysF_}fq7lvNgOA8=o0Lf7L@SbAQ?2szNlr6`Lg
z0I6&>cZ%Ed-#|-BO7Tn@YI~DS!;kA#@7jiSt}{QUVrD(~rVcULTJLX
z+7h;1*(St9rlozQU`jC89Yf_B>QiuYa<&4RGL=*_zZDq{YVaDQeczcy5q#CgDkvyf
zWvBOC$X?5QEj785)%N}Ga3Z@b(B0?^Z|V*pLq??_;>?@PIWtcvwXmq$IVTowKFq*?
zLbQ5^fM#B_ZI4_?##<~6ssI3jK0U`uT44R+)VcayP=A4*z()7VQ=sFEzDg9WyMrz(
zU&^nJx7nmzZx5fAnDhMRg_cBtC+^oL~@7w3kL5`FKp&T5_P|jWs2iZ>dd2$8mZE3E!`IM{a4JI>3@D*xtmyt+j!2
zGIYxyEd>0f{UcBS1?bFtLN4@(5pxmLyq;io*b`v$j@0KA@bzHaWJn|=w?ojFWbKQo
z>~o@_qvH2W6oz(if{^)jGyCAdui~cqLQc+{AjJ_dRr7nOPo7i^eyw*fpw!i1i~j7W
zDk1S(LV%rE4s^&0JaK|~WT!Pyo5W!%J6=zBR?(gpJi|n@40AlmbeM@&(ykVM#arB+
zuq?B+DySoVk)bGK)C5nQT)#SfB`9G3$z0g(6Uh!aVLDi~_PBR7vSts*ZU
z?z+Y^Yzf@|KTmF90Gdlvo#-r)K*%3+#1&=1J&jHa3XY32BG}PW{3&4y$`NsAv;~lN
zW`V_AjT*Si*%r=(=Z|a8qhS~3WVB!+DG}C^9h&T2`!r~^SB=dZ8jhY|N0QLp$X61z
z{4+eru}6>U_ic0sw$YDnEdam<0=loll#e@+Hdb3g7^8W12!4jYET?8u=1s}6<^{86
zPN6=?0cDp0&=;0@r&Xmd&qVd9ESQo_+u`5RKkp0WTx-(R=smMRqj~=P8SZbC!30`8
zw${=_%xJZpu9)u9x@|(pYtlCf-JOJP&ORPw8b5*7ie1~UU?r`5nWk>9Y5fE9S;3jY20
zVG4L`$CLnCBR>nnq+M>5eW=jQ+@w!)6`v@t0&i(AUu&Es7qeDVdx!MZ+2u2IAYe^=
z1-m!n)QQgFcA0bTRv5otwQwbzwh{^SyHG;;2>R}znsR<_p({{DxZhod5g?){@xmKG
zZaIYpbiJ@QY7X}7x8ZSGM$fa2N5?;-OHl^Lm6#VKP0$G7_Be8(ofWmXQemrrOiWCi
zwwzL)aZ4n8pexGBYof|2b$CHw^2vXE=_=n|pT^gU9@i)>w=@2E`?RzGxx}BmTerp$2l}F4J{yo}ckdxqrWNhVo%TjM&iv1T7prWQ#&%dJfylmD(XoU1daZyOs#=rSeDdNB`C
zaht0NH-1jAEe2mo^5+Dh6p)Zwg;t>CC*GoK)PjyuwPRGK+$f4>ZG#z
z=`3vR#mV~d^OK}METChFIM70f=Qv_?D-g42=Hg=P>N+x~Csf5!{+?AZmKq3?{-j|*
zmnlHVvi?RgK!`%%{sJsh?cE^Oki5l3ODr^a6!B`NN5d#+Nr8S|y%d-VqGTs?nepDu
zPHmp9$n>g{e(|aAoJDah5UQbQpy5POBLHaXpL)TU2mxh^sKw)khlNk27`FcO!E^7f
z{CNMKbBmnNO)>^F`l%QwFH~j8+0WKU;{%c>-YYql*|oY^W+6I5mH7A@pGey7R#x;Mqp-#5#4W7YAwaKxA#|;bClA^GZJSE6G&K
zQg0p&b;X7uEHH{a5_Ij~pZ+bubV!yUw<^noY-qJu3apz>mjV6pLIZDuzsR@(@i$f-
z0u~To?1K|fYH&;l5Iww}uctFRkWIjEuetug@zrTQS+ceR{Q#23t@-?C+nj5gCamt-eLyqiJn`VB`YAt
zvYEZYo7a0nbgL~r+fn3FVnA~vcR{wvh67JqZM(GWHPY8_fAG6s$H>K2>7`*o!aVJ!
z)1ma)d%{vk4ZD=6|CG*5!*}tYKao#|uKU^|e%l`Yz)S?4y1KfYrly5=6x_&xM#zte
zpc7|^37oP2k1BNokDHIrMAfM{@jr7R9y~*+WpLU7+>+Ueio9Y3nt>EQ(A5S=cMb^d
zLmT8lI&N1~jWVB!Sn?hr062-g
z02)H_tXHv*@c}s_fOnA*7_^?IGyrX2WKG^8unAi+KY%Wh|41XhT7Y<%lmq&J&ZGLS
zANcKH8xN(x^ku!?Z+-j1P+K+HX9h9BGu~GMS!k{eRoxxYCNN>!)_TOQjA89pFe-fOeVLV
z3)$<@{^8V-|LxSZ*g<9q$*XC;W{(Q?(V8s*uS5HI)0f0-ZFCpG*|7HYpT*I!8W4N;
zMKK-n=Z{SDqk`67uU6Cc+<7DQnt6D*>?65MJq9Qw8z|XJ%&{5S+fiz6u7Ru&_JWaQ>cs{p`n84OwTYIb7C65*6d1$;Y%6?Dd
zayE&qD!`~f7!*s^1XY6VwGt6Cr%M;*&9ilG&+aNat>=af?jKfbV6}J!m*M-XY4Gs={Lv5)g9A2QJPzzw~?9XV}^&eFnx1
z@vM+H_h0Fxaf(DK+w?el`Y)35sD?)uOg`b&8+p0M%gau+xVXI=;1yI`C!I>@-Z(Iz
zvyy*TYQ=}~XXrIB2#A->$PP|&6%&9r6S~gPfwwh)zDSe$2nBfTLl$=Hg}J*W$6Voy
zbq*`5gR{mgv?dZMDciSQ55CeN$CYAqUjZ_VC}md-n(d85^53
z!vmd#O1ev%GdiB<6W>%+7diHP=gNLo;R>
zF~GG4m$^~B#__3yRDgAz^pc<)J$1|ak6u0=wLZQggmry)cc-GQ%^+xAWr2c)m9V`0
zcI1{B;AdjGxQuJq?Q!Kdr|b8=hfZFy?3jZ8~ao&Bd(1SI~k;GCvu#^CmM|vY4G$YLbvWXdn|fobw4G(FM=e2LpYJ?
zfa1&Jy9jeEnQftacbm}-_pREsYS8ulu%zMZ=nkS|*^HUzOra3k@nf7nTS08(@JFvA
zwBu+jXq%`%0uzR)1B46gwc-vBQ4Od2@}6zN@BoUG=3|mj%%T%Zn=#v@U_VzTO}N}n
z!}rDfG1;PIb*#&lFz|l}eeXVE
z0?>C;Q@B-u_eg~iEXp{KYwu3|5P;G!nL0YQpJ|d25wSidN*T>ICDff=HTPaMk2%Bt
z8ZO7CM}p_Sjf}iuU}NKRSj3tk`vLS8FR!fFO}zok=yUQU>5HKTCuvqbm%4N+Mn*mR
z89nBVCfD3m3of>fM{f9Y*^o(UnPwjW*fm9IpxqhHeu>c**+W<
z{oV!3{_2rzcwM0T#m?@0MbE<2RPH6FL-|cea@)s5ab>As;k%nQr)8hC9k3EPkSi`^
zHUcso9sP9YdgmAD4+`+}T84)|Pt>oea}{g(swV8<97AGOkC3Xm*VWN~?OdnyApA+H
zmtH?65;Yw2^Ybd|>So?KKp9!P&ZY2UP|!?iK5k$TW2VD^=)>eYPYZ{n@MOxiIXzC5
zl{LQ@|FuC91sY!3(SZa>C}7F}R@X4H18)KW;7?zlKR?Tw`ndoy*namt)=vCJ&l4>B
zukW9ErPM9CjY$8A)p|65DT%1Ifl#(_*JMlJV^jf+8%v=9alsw%n(2$4N_E~y`UBA!
zG0a#CO)u;h)Tq(rm$!B0SN?cp9yY&gmVw{H0>cZtKjZxYHjl3ENJE*|RLjk+%GBI^@Cbc3LVOc>Qz|*5*+tHJTS(9MjUw9ETDf~^
zT!{_-i0;6Idvxa?F?_#TW_F#70T2^O<$nj>*2(i0vbf3lNO-mEOH#a0@at6Aeeg>b
z&7Yik4dBfCp4scYdKIZx=hAbI{4X=rE~TYQ%1DqbuBsb(dmIPEBp#JTY`!w@;(h$O
z1hpOhA6Z`+Q03NjE1QsRL8Oro>5%RaX+)$;TIudGKtQ^sK?J3{ySqd>C8fLJu8l|G
z{qBF*&x#pijv32`MMTC$P^5hOXooV~=2XCN?P=a6L
z;GD}>te=xj-)p_`q~Ge#Z!L)?r=)iwm%CYgQS|ZE8gFU&qF8M>r0BfvEgKev%SkFj
zgmd#*KqP2F+Z8RXfF&+29;JadY|mFqE!231M8vOQiiY!ZbaG>SZ9Jbo&ENb*r;mRrEpK3ro>UL0owkm7&e^$X9;2@-Pfk%4|=qrQ5M^
zM`^m^yym(_N(ykFd$wRcHyMryA0cSsEq-ku>9*3_H$dxD`51nS7n+pY60>0#P#C#hW5)-?ZXpC7i%%&*qEF<
z&3(tW09qI{p9jrW>iOVI)H;82I~j_jqVk*B&~m&oVZOCH7(C!!M7MwbRQ{HiS2jh*
z#`eN)lo=U2D$a6!!b|HW?D@3C=X<=3zzYoMPpxtW9Ntn1Ofj%awtCVB55jhLZT2Km
zk2IVI1_l6S+(@d+wr%j@2$Yo_busC1Q&h=9G_7ZfaA+C3^C9$W(aRNAMpp2DE{98X
z)<)l#uWs-CIq
z`g>foO4zKMC9E+%1XX;SGqr;cjM%Ap2JJU88gAqUI4e(JS--xXn3`gE`SN9{M#h&f
z&j5A3!Sfgk@rn1i5Bf<*z}*Ez{W65SusO4k*D!^jyePTZ!r4yp%o}1^f)+`
z0JU`7Txc;(YHUN;3#J(6k%4q47oO;!*8Oh)w$L*u_tm=bqwL7OJ-Tes;TZv(FueRp);@vj=z;`k}=?G4oV6t#U#0I>`(M+MDfNk3Yzu-Wd9pdEGKj-y6w?5PP+p0So
z0^ruBH-6S&Zx~0_vA5Y#IhzrXWJ=>IYK*LHt9572J3e7-(D+0*oCY15Y5
zZ7nU(Zl??`NEVh=p3qf&{z24gFaNSZ_gRvcK9fwFLft{7o+;`9C?IJo4fpv^AnR!$
zu%ZulTQ`e#cU?2SlzFa%-(x&3ws6ip?`YzDrtFf6qArbd>&rs*Gtg_r%*n}#3_?H4
zOh8<`^WQ9h6jc9!04k2C6v$UHgA!?Yxu{{V_Zg<19RA$Ow}Dx^iv{xLVU3eYDkX8P
z|EYcbJ@4E^80Vc0oZa49Ruwdr2ROvE!P9lZn>EhxH5$jPx2T-_df-ZAY?@qMU9U@y
zG*+PK=svu+0!}~YGtZsQSI%2P740#uxfYVbpU($i;|yhA=ejPT)ykxk!~3%n(}$xy5rVqs|+-4)F=kvp77Kc4pCtxzVS31gHFKt0oA
zI!m89pX%#VKrwn(hNT7r%hD@^znUVVqHAU*b4XG}S?99x@U?|)U;bA*|DuIH(lKff
zO{Tcomb3ozT405k5C#q%h^(6@=%!<%-p?Ny2d&?4l
z>dU~+Q(c~vjp7_Z*82+(Ax@G=6aGF%i1|&-7OCaSmzRAYQ@jK}2TFkOM
zU%7Y`W2cqGyVrl?_3i*64V9$x(}Q4xqdDY}bCRN0hKU_hg@O3uON)dJ!nXt!F&6+z
z49|_0)Ec#7l>lQ&DXFl=Mv=ZZQCL<9^)2JdCzpAL&t%BbXd@tA_$U5X~lJy9qB4Xlb-C?yuGtKkRcGL%)`9AhLPqZ
zC3^gm!{kk=X=yb-pL=_0z&-6bUSg_qRM&P|D<^o+#~ZU!yw)#P_}O=n%`uxCBVfn^
zed7B$-uc~Bh#b+q2?Qmsu)VBF^f
z(9qGPoSiFHNKF%|`1!SLY;F0o)SlKqfs)Y5On7D0JI+m~A}71B8kBT_DW{ryRB0)H
z><`>?U0o|xkE32m=ivn_>n0Sc)o<3j%fZu>rqU;b4?pS<-WF9a9jF@WW`p4Yy}i9i
zPgDdznO|3-VHYI@MWZQb&?*Sa0CzpQ9yu
z)Zag8h(_bEFOvhS;bdyV(OuO6DW)aom{MwtF-1vXe8+8EOox~J<|1JfpwrZ>%K75!C
zoJp|?3Oo!a5-3k-cl0Qg9_L9ARIUhmaam44O%K#Vw)syG1W+OoBQrZY1yhZql=fv4
zEx`jD2an43b}JBZ#w}cK&}E;8pl75z|BufERCPeZ%Gz2us3At8Pjc(->IyGrpQb+(
zdoM4C-%ljUbT#)0IQ^JkZ8Uv-avXlTWt|UNsU1&8{Qtukvv0h8q8LwtXNbe}jhT2a
zd44kVmDrSRL;QbP-a~`N-~XNL?029bqSmW}RI95yB#NN+qVsnTbZ)@CPZ-csO;1zt
z*v{g&`7QD;JdSq22z9q_szjNj+`L4LRVF5F>t3Mf^^cK7z@f)y7W&|rnmalK&p^@n
z5YzOdksLK#sdG(~TfUNL>Lcj@l3ft?*aW(H_^-yAQM
zfx64=XeA?eOovx|0eF7Il;^H5=eqoq}R>@-E`#edkCDgbD
zo-rSJ!xNi)`~E#_bJK!p?Up&GL8`7fs^;6<#UfzmSC3cUwGCU8ud_-60gcTYoY((4
zoyd9sol!JH9Fdk`N}OXoH7>quOzkMG=k_UWLh{!weFS(Mg`C`9FOajHzm39h6xv!-<=21`Bh5|lVH|M>2^>?t3CI#q=)pbE5g%bb76Z6~K%+AM~
z2fG!-ruU_~3mpA!eTkmZ8=k;~fQTpp1)p(rding^O;VO>;0ptKq%h8XwqL$==e}&a
z@dScy1_taV8!_w61JAK`Ke&!`L6c8&uU>_lck=wzulgMbNW?`ZRa*(EkGm32{Cv`7
zvk(_NMgL9GI>`JX+8^z?Rv}$-MEZ6iX)UdU;*t{QsTojPNz1^nn8VCO@c1ol|Lf`R
z=SMmlp+KMS77>Iey0#3`n^tdEXD>*W>Hl?&WV6)QIlvY=UqDDNQt0=ZIZG*ix@V?V
zMc)^T_T}X*38lZ|4Uca^Ktd8#Sjaj+0!gIepvD$PrBMHLx)?qgUm<+Z?zYb;lkr7C5!JB
znt0VH;sPGU?Q6LLG&-E^>rQS^JfnOBs9F0bCJ*K%Z7f9L3D}psb=&J*KNbsH*pW
z3npjkADChibbeK5(%<4oHCc_i2bC=2E!C#&T)ArDVx)hjn`{9Nd@jAZS`z`GdWKJ0
zS~|R%FHFK`vTVW-zj&`{Z}*dADd1&(mq)FY&w!vqiloN!omgx&lr#OyxrH8Kyln5&
z;h}WD5Di)wE&eWaIWIi{2B5TtZlpcGxjV%7{`GH>iud>T+o8nZmIfFa(%ih=<_ol9
zkFj5n*$hpvU^(;R$ab2Lxdn-JXnwJev_?`nn)$pLS1B9He}Zp%!$0abaG|(5%R6iz
z^M5L#?QZ1q)@?1gRNDHL&%+6cw14Y@!$gGlhQ~4amYW-uo6DHz8Wj`s2$UEe0%+IR#KIADgkm-9{nvB
z=y6XK^PeKT4G0Q~8dRMAz56yrPa7VLD_I!HBaH-U3@n1UkrzF`XS1Yj7nn$P@5A;yrUFUZ2aedNuI<#)>9)nc))p@gdSY+##3(ARWXUR?AbgV+Q0Lw5Q
zsN-6Mis5?npuM9b^b>09#6%|i?xK{roPn}=P}9MXkZ?t@v)rZ
zaVW1VSgxfI2M%;_EK1ak(Ma_Dkju1hZG+lyz{c#xPN?jK1+rNn<@^b;TrX~=k!wlCJuTyy+!)hyZ?$|7O#BxM21Ipd1^QRY7z0myC_2syC{7_3z2gI1KL
z3W(h1mX@)cRCJg+B+wkz$;Lz)$|=aVT&91zCiHpxQl6w%*g~~#RL6IL>*jL;aNmEA
zFSyp7D*$-lW`EFId+U-JH>;ZE$uhzo~
z@wnITMt`-pI%V^9()mPHnZ75?@^eZftK{3?Bd)y-|AeOl4x}m$F=DUI(wbEd;
z&$<1xIwFa-`AgvE7gMT6B;TLN+#C>Mb-g#v%iUkgN7C|nM>CJTEqiE;f6Uw?J=>=~
z+0V%m__l!aSG=GQPlgt=&^u~+TD!t&Cly}^bgt0@TS`GeQ4Nix1mO;e^TJ3kUGh*-
z(a(NadGJ1Uk%~Xp(WzAf8HnSJ7@`{y)OEN{(A0UHct>cNB~HOCS?6tNM5&8*KNa8M
z9O7fOy{Rhl#>7S(6Y5&*icYVonOak(qM_*mopRN`K{Vc(r8M22J(oBv#7~hIN>7
zGu0FQAqk+$LiPnt1T=kORQDzD(g5%0nS+amC#8rjNh11NW8EuYV47~`YM2CysbF$W
zz92!p)yf57_8e(3C;s2)9~6Q(L)!XAI*yVq6nLuc5)E^0VC6JK%{QV-lmxWiSShku
zcJC0>Vmk*1DOuUT=n{1NuT4HX^WLJekK|2}
zgHE9fPphDtNPbquyEIqeWCn<Gh+H}uqhMV%A5P{)JB8>g1O&?bCU-7ub
z7}Gv}9Cm!{1daS-BI^@_89ReLisXR~7?oUAuNzzbx$CK=*}uU$S$5H%ymiUC2wS
z)=~`9BOYda@a{i+0RCQx8M-q?EiDEx5KsH-pB(V7m-6rsHmv&(NJ2^Z%j5<4Y^