{
"cells": [
{
"cell_type": "markdown",
"id": "253688c6",
"metadata": {},
"source": [
".. index:: Equations; Darcy flow\n",
"\n",
".. index:: Methods; hp adaptive Discontinuous Galerkin\n",
"\n",
"# Problem Description\n",
"This example considers two-phase porous-media flow.\n",
"In all that follows, we assume that the flow is immiscible and\n",
"incompressible with no mass transfer between phases.\n",
"A detailed description of the mathematical model and parameters\n",
"can be found ."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "8d17b3e7",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:29:15.529671Z",
"iopub.status.busy": "2024-02-29T12:29:15.529327Z",
"iopub.status.idle": "2024-02-29T12:30:54.268616Z",
"shell.execute_reply": "2024-02-29T12:30:54.267667Z"
}
},
"outputs": [],
"source": [
"try:\n",
" import dune.femdg\n",
"except ImportError:\n",
" print(\"This example needs 'dune.femdg' - skipping\")\n",
" import sys\n",
" sys.exit(0)\n",
"\n",
"from matplotlib import pyplot\n",
"\n",
"import numpy\n",
"from ufl import *\n",
"from dune.ufl import Space, Constant\n",
"\n",
"import dune.fem as fem\n",
"import dune.create as create\n",
"from dune.generator import algorithm\n",
"from dune.common import FieldVector\n",
"from dune.grid import cartesianDomain, Marker\n",
"from dune.fem.function import levelFunction, gridFunction\n",
"from dune.fem import integrate\n",
"from dune.plotting import plotPointData as plot\n",
"\n",
"from limit import createOrderRedcution, createLimiter\n",
"\n",
"fem.threading.use = 4"
]
},
{
"cell_type": "markdown",
"id": "ed1b9124",
"metadata": {},
"source": [
"Some parameters"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "12d92d6a",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:30:54.271980Z",
"iopub.status.busy": "2024-02-29T12:30:54.271667Z",
"iopub.status.idle": "2024-02-29T12:30:54.276010Z",
"shell.execute_reply": "2024-02-29T12:30:54.275347Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"maxLevel = 3\n",
"maxOrder = 3\n",
"dt = 5.\n",
"endTime = 800.\n",
"coupled = False\n",
"tolerance = 3e-2\n",
"penalty = 5 * (maxOrder * ( maxOrder + 1 ))\n",
"newtonParameters = {\"tolerance\": tolerance,\n",
" \"verbose\": \"false\", \"linear.verbose\": \"false\",\n",
" \"linabstol\": 1e-8, \"reduction\": 1e-8}"
]
},
{
"cell_type": "markdown",
"id": "6e448e33",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"# Defining the model\n",
"using a Brooks Corey pressure law"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "57c71d3d",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:30:54.278526Z",
"iopub.status.busy": "2024-02-29T12:30:54.278297Z",
"iopub.status.idle": "2024-02-29T12:30:54.284397Z",
"shell.execute_reply": "2024-02-29T12:30:54.283554Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"def brooksCorey(P,s_n):\n",
" s_w = 1-s_n\n",
" s_we = (s_w-P.s_wr)/(1.-P.s_wr-P.s_nr)\n",
" s_ne = (s_n-P.s_nr)/(1.-P.s_wr-P.s_nr)\n",
" if P.useCutOff:\n",
" cutOff = lambda a: min_value(max_value(a,0.00001),0.99999)\n",
" s_we = cutOff(s_we)\n",
" s_ne = cutOff(s_ne)\n",
" kr_w = s_we**((2.+3.*P.theta)/P.theta)\n",
" kr_n = s_ne**2*(1.-s_we**((2.+P.theta)/P.theta))\n",
" p_c = P.pd*s_we**(-1./P.theta)\n",
" dp_c = P.pd * (-1./P.theta) * s_we**(-1./P.theta-1.) * (-1./(1.-P.s_wr-P.s_nr))\n",
" l_n = kr_n / P.mu_n\n",
" l_w = kr_w / P.mu_w\n",
" return p_c,dp_c,l_n,l_w"
]
},
{
"cell_type": "markdown",
"id": "c30ab8ec",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Constants and domain description for anisotropic lens test"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "feabf6e4",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:30:54.287334Z",
"iopub.status.busy": "2024-02-29T12:30:54.286898Z",
"iopub.status.idle": "2024-02-29T12:30:54.296310Z",
"shell.execute_reply": "2024-02-29T12:30:54.295486Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"class AnisotropicLens:\n",
" dimWorld = 2\n",
" domain = cartesianDomain([0,0.39],[0.9,0.65],[15,4])\n",
" x = SpatialCoordinate(triangle)\n",
"\n",
" g = [0,]*dimWorld ; g[dimWorld-1] = -9.810 # [m/s^2]\n",
" g = as_vector(g)\n",
" r_w = 1000. # [Kg/m^3]\n",
" mu_w = 1.e-3 # [Kg/m s]\n",
" r_n = 1460. # [Kg/m^3]\n",
" mu_n = 9.e-4 # [Kg/m s]\n",
"\n",
" lens = lambda x,a,b: (a-b)* (conditional(abs(x[1]-0.49)<0.03,1.,0.)* conditional(abs(x[0]-0.45)<0.11,1.,0.)) + b\n",
"\n",
" p_c = brooksCorey\n",
"\n",
" Kdiag = lens(x, 6.64*1e-14, 1e-10) # [m^2]\n",
" Koff = lens(x, 0,-5e-11) # [m^2]\n",
" K = as_matrix( [[Kdiag,Koff],[Koff,Kdiag]] )\n",
"\n",
" Phi = lens(x, 0.39, 0.40) # [-]\n",
" s_wr = lens(x, 0.10, 0.12) # [-]\n",
" s_nr = lens(x, 0.00, 0.00) # [-]\n",
" theta = lens(x, 2.00, 2.70) # [-]\n",
" pd = lens(x, 5000., 755.) # [Pa]\n",
"\n",
" #### initial conditions\n",
" p_w0 = (0.65-x[1])*9810. # hydrostatic pressure\n",
" s_n0 = 0 # fully saturated\n",
" # boundary conditions\n",
" inflow = conditional(abs(x[0]-0.45)<0.06,1.,0.)* conditional(abs(x[1]-0.65)<1e-8,1.,0.)\n",
" J_n = -5.137*1e-5\n",
" J_w = 1e-20\n",
" dirichlet = conditional(abs(x[0])<1e-8,1.,0.) + conditional(abs(x[0]-0.9)<1e-8,1.,0.)\n",
" p_wD = p_w0\n",
" s_nD = s_n0\n",
"\n",
" q_n = 0\n",
" q_w = 0\n",
"\n",
" useCutOff = False\n",
"\n",
"P = AnisotropicLens()"
]
},
{
"cell_type": "markdown",
"id": "060216e1",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Setup grid, discrete spaces and functions"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "aeaeb1b0",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:30:54.298912Z",
"iopub.status.busy": "2024-02-29T12:30:54.298667Z",
"iopub.status.idle": "2024-02-29T12:36:54.660084Z",
"shell.execute_reply": "2024-02-29T12:36:54.658672Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"grid = create.view(\"adaptive\", \"ALUCube\", P.domain, dimgrid=2)\n",
"\n",
"if coupled:\n",
" spc = create.space(\"dglegendrehp\", grid, dimRange=2, order=maxOrder)\n",
"else:\n",
" spc1 = create.space(\"dglegendrehp\", grid, dimRange=1, order=maxOrder)\n",
" spc = create.space(\"product\", spc1,spc1, components=[\"p\",\"s\"] )\n",
"\n",
"solution = spc.interpolate([0,0], name=\"solution\")\n",
"solution_old = spc.interpolate([0,0], name=\"solution_old\")\n",
"sol_pm1 = spc.interpolate([0,0], name=\"sol_pm1\")\n",
"intermediate = spc.interpolate([0,0], name=\"iterate\")\n",
"persistentDF = [solution,solution_old,intermediate]\n",
"\n",
"fvspc = create.space(\"finitevolume\", grid, dimRange=1, storage=\"numpy\")\n",
"estimate = fvspc.interpolate([0], name=\"estimate\")\n",
"estimate_pm1 = fvspc.interpolate([0], name=\"estimate-pm1\")"
]
},
{
"cell_type": "markdown",
"id": "c49adf74",
"metadata": {
"lines_to_next_cell": 2
},
"source": []
},
{
"cell_type": "code",
"execution_count": 6,
"id": "89b7eeaa",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:36:54.664926Z",
"iopub.status.busy": "2024-02-29T12:36:54.664396Z",
"iopub.status.idle": "2024-02-29T12:36:54.675381Z",
"shell.execute_reply": "2024-02-29T12:36:54.674140Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"uflSpace = Space((P.dimWorld,P.dimWorld),2)\n",
"u = TrialFunction(uflSpace)\n",
"v = TestFunction(uflSpace)\n",
"cell = uflSpace.cell()\n",
"x = SpatialCoordinate(cell)\n",
"tau = Constant(dt, name=\"timeStep\")\n",
"beta = Constant(penalty, name=\"penalty\")\n",
"\n",
"p_w = u[0]\n",
"s_n = u[1]\n",
"\n",
"p_c,dp_c,l_n,l_w = P.p_c(s_n=intermediate[1])"
]
},
{
"cell_type": "markdown",
"id": "1145820e",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Bulk terms"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "ced270bb",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:36:54.679404Z",
"iopub.status.busy": "2024-02-29T12:36:54.678871Z",
"iopub.status.idle": "2024-02-29T12:36:54.689687Z",
"shell.execute_reply": "2024-02-29T12:36:54.688456Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"dBulk_p = P.K*( (l_n+l_w)*grad(p_w) + l_n*dp_c*grad(s_n) )\n",
"dBulk_p += -P.K*( (P.r_n*l_n+P.r_w*l_w)*P.g )\n",
"bulk_p = -(P.q_w+P.q_n)\n",
"dBulk_s = P.K*l_n*dp_c*grad(s_n)\n",
"dBulk_s += P.K*l_n*(grad(p_w)-P.r_n*P.g)\n",
"bulk_s = -P.q_n"
]
},
{
"cell_type": "markdown",
"id": "c5e5cc07",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Boundary and initial conditions"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "df36a387",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:36:54.693814Z",
"iopub.status.busy": "2024-02-29T12:36:54.693280Z",
"iopub.status.idle": "2024-02-29T12:36:54.699548Z",
"shell.execute_reply": "2024-02-29T12:36:54.698302Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"p_D, s_D = P.p_wD, P.s_nD,\n",
"p_N, s_N = P.J_w+P.J_n, P.J_n\n",
"p_0, s_0 = P.p_w0, P.s_n0"
]
},
{
"cell_type": "markdown",
"id": "b862ba13",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Bulk integrals"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "a8a97cc3",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:36:54.703807Z",
"iopub.status.busy": "2024-02-29T12:36:54.703042Z",
"iopub.status.idle": "2024-02-29T12:36:54.714187Z",
"shell.execute_reply": "2024-02-29T12:36:54.712932Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"form_p = ( inner(dBulk_p,grad(v[0])) + bulk_p*v[0] ) * dx\n",
"form_s = ( inner(dBulk_s,grad(v[1])) + bulk_s*v[1] ) * dx"
]
},
{
"cell_type": "markdown",
"id": "e87d8879",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Boundary fluxes"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "68ada3bb",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:36:54.720206Z",
"iopub.status.busy": "2024-02-29T12:36:54.719666Z",
"iopub.status.idle": "2024-02-29T12:36:54.728081Z",
"shell.execute_reply": "2024-02-29T12:36:54.726798Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"form_p += p_N * v[0] * P.inflow * ds\n",
"form_s += s_N * v[1] * P.inflow * ds"
]
},
{
"cell_type": "markdown",
"id": "cd2a58e7",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"DG terms"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "6c70e378",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:36:54.733924Z",
"iopub.status.busy": "2024-02-29T12:36:54.733374Z",
"iopub.status.idle": "2024-02-29T12:36:54.743921Z",
"shell.execute_reply": "2024-02-29T12:36:54.742634Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"def sMax(a): return max_value(a('+'), a('-'))\n",
"n = FacetNormal(cell)\n",
"hT = MaxCellEdgeLength(cell)\n",
"he = avg( CellVolume(cell) ) / FacetArea(cell)\n",
"heBnd = CellVolume(cell) / FacetArea(cell)\n",
"k = dot(P.K*n,n)\n",
"lambdaMax = k('+')*k('-')/avg(k)\n",
"def wavg(z): return (k('-')*z('+')+k('+')*z('-'))/(k('+')+k('-'))"
]
},
{
"cell_type": "markdown",
"id": "af48685a",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Penalty terms (including dirichlet boundary treatment)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "24e1289e",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:36:54.749749Z",
"iopub.status.busy": "2024-02-29T12:36:54.749214Z",
"iopub.status.idle": "2024-02-29T12:36:54.770255Z",
"shell.execute_reply": "2024-02-29T12:36:54.768955Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"p_c0,dp_c0,l_n0,l_w0 = P.p_c(0.5)\n",
"penalty_p = [beta*lambdaMax*sMax(l_n0+l_w0),\n",
" beta*k*(l_n0+l_w0)]\n",
"penalty_s = [beta*lambdaMax*sMax(l_n0*dp_c0),\n",
" beta*k*(l_n0*dp_c0)]\n",
"form_p += penalty_p[0]/he * jump(u[0])*jump(v[0]) * dS\n",
"form_s += penalty_s[0]/he * jump(u[1])*jump(v[1]) * dS\n",
"form_p += penalty_p[1]/heBnd * (u[0]-p_D) * v[0] * P.dirichlet * ds\n",
"form_s += penalty_s[1]/heBnd * (u[1]-s_D) * v[1] * P.dirichlet * ds"
]
},
{
"cell_type": "markdown",
"id": "c9776562",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Consistency terms"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "aad0b5ef",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:36:54.776399Z",
"iopub.status.busy": "2024-02-29T12:36:54.775841Z",
"iopub.status.idle": "2024-02-29T12:36:54.799046Z",
"shell.execute_reply": "2024-02-29T12:36:54.797782Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"form_p -= inner(wavg(dBulk_p),n('+')) * jump(v[0]) * dS\n",
"form_s -= inner(wavg(dBulk_s),n('+')) * jump(v[1]) * dS\n",
"form_p -= inner(dBulk_p,n) * v[0] * P.dirichlet * ds\n",
"form_s -= inner(dBulk_s,n) * v[1] * P.dirichlet * ds"
]
},
{
"cell_type": "markdown",
"id": "6f19956b",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"# Time discretization"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "114ea6e5",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:36:54.805217Z",
"iopub.status.busy": "2024-02-29T12:36:54.804657Z",
"iopub.status.idle": "2024-02-29T12:36:54.816831Z",
"shell.execute_reply": "2024-02-29T12:36:54.815531Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"form_s = P.Phi*(u[1]-solution_old[1])*v[1] * dx + tau*form_s"
]
},
{
"cell_type": "markdown",
"id": "b1981ee5",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"# Stabilization (Limiter)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "5a1a9e68",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:36:54.825271Z",
"iopub.status.busy": "2024-02-29T12:36:54.824713Z",
"iopub.status.idle": "2024-02-29T12:37:25.011600Z",
"shell.execute_reply": "2024-02-29T12:37:25.009975Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"limiter = createLimiter( spc, limiter=\"scaling\" )\n",
"tmp = solution.copy()\n",
"def limit(target):\n",
" tmp.assign(target)\n",
" limiter(tmp,target)"
]
},
{
"cell_type": "markdown",
"id": "e16ec73a",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Time stepping\n",
"Converting UFL forms to scheme"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "2a9b4783",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:37:25.018336Z",
"iopub.status.busy": "2024-02-29T12:37:25.017803Z",
"iopub.status.idle": "2024-02-29T12:38:37.954927Z",
"shell.execute_reply": "2024-02-29T12:38:37.953116Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"if coupled:\n",
" form = form_s + form_p\n",
" tpModel = create.model( \"integrands\", grid, form == 0)\n",
" # tpModel.penalty = penalty\n",
" # tpModel.timeStep = dt\n",
" scheme = create.scheme(\"galerkin\", tpModel, spc, (\"suitesparse\",\"umfpack\"),\n",
" parameters={\"newton.\" + k: v for k, v in newtonParameters.items()})\n",
"else:\n",
" uflSpace1 = Space((P.dimWorld,P.dimWorld),1)\n",
" u1 = TrialFunction(uflSpace1)\n",
" v1 = TestFunction(uflSpace1)\n",
" form_p = replace(form_p, { u:as_vector([u1[0],intermediate.s[0]]),\n",
" v:as_vector([v1[0],0.]) } )\n",
" form_s = replace(form_s, { u:as_vector([solution[0],u1[0]]),\n",
" intermediate:as_vector([solution[0],intermediate[1]]),\n",
" v:as_vector([0.,v1[0]]) } )\n",
" form = [form_p,form_s]\n",
" tpModel = [create.model( \"integrands\", grid, form[0] == 0),\n",
" create.model( \"integrands\", grid, form[1] == 0)]\n",
" # tpModel[0].penalty = penalty\n",
" # tpModel[1].penalty = penalty\n",
" # tpModel[1].timeStep = dt\n",
" scheme = [create.scheme(\"galerkin\", m, s, (\"suitesparse\",\"umfpack\"),\n",
" parameters={\"newton.\" + k: v for k, v in newtonParameters.items()})\n",
" for m,s in zip(tpModel,spc.components)]"
]
},
{
"cell_type": "markdown",
"id": "704ed3d6",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Stopping condition for iterative approaches"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "b95fb555",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:38:37.961427Z",
"iopub.status.busy": "2024-02-29T12:38:37.960904Z",
"iopub.status.idle": "2024-02-29T12:38:37.969286Z",
"shell.execute_reply": "2024-02-29T12:38:37.967745Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"def errorMeasure(w,dw):\n",
" rel = integrate([w[1]**2,dw[1]**2], order=5, gridView=grid)\n",
" return numpy.sqrt(rel[1]) < tolerance * numpy.sqrt(rel[0])"
]
},
{
"cell_type": "markdown",
"id": "07287e6c",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"# Iterative schemes (iterative or impes-iterative)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "4aac5bf2",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:38:37.974997Z",
"iopub.status.busy": "2024-02-29T12:38:37.974493Z",
"iopub.status.idle": "2024-02-29T12:38:37.983886Z",
"shell.execute_reply": "2024-02-29T12:38:37.982444Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"def step():\n",
" n = 0\n",
" solution_old.assign(solution)\n",
" while True:\n",
" intermediate.assign(solution)\n",
" if coupled:\n",
" scheme.solve(target=solution)\n",
" else:\n",
" scheme[0].solve(target=solution.p)\n",
" scheme[1].solve(target=solution.s)\n",
" limit(solution)\n",
" n += 1\n",
" # print(\"step\",n,flush=True)\n",
" if errorMeasure(solution,solution-intermediate) or n==20:\n",
" break"
]
},
{
"cell_type": "markdown",
"id": "cbb1f35e",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
".. index:: Adaptation; hp-Adaptation\n",
"\n",
"# HP Adpativity\n",
"\n",
"Setting up residual indicator"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "cb4a4c48",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:38:37.989418Z",
"iopub.status.busy": "2024-02-29T12:38:37.988921Z",
"iopub.status.idle": "2024-02-29T12:39:24.971512Z",
"shell.execute_reply": "2024-02-29T12:39:24.970007Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"uflSpace0 = Space((P.dimWorld,P.dimWorld),1)\n",
"v0 = TestFunction(uflSpace0)\n",
"\n",
"Rvol = P.Phi*(u[1]-solution_old[1])/tau - div(dBulk_s) - bulk_s\n",
"estimator = hT**2 * Rvol**2 * v0[0] * dx + he * inner(jump(dBulk_s), n('+'))**2 * avg(v0[0]) * dS + heBnd * (s_N + inner(dBulk_s,n))**2 * v0[0] * P.inflow * ds + penalty_s[0]**2/he * jump(u[1])**2 * avg(v0[0]) * dS + penalty_s[1]**2/heBnd * (s_D - u[1])**2 * v0[0] * P.dirichlet * ds\n",
"estimator = replace(estimator, {intermediate:u})\n",
"\n",
"estimatorModel = create.model(\"integrands\", grid, estimator == 0)\n",
"# estimatorModel.timeStep = dt\n",
"# estimatorModel.penalty = penalty\n",
"estimator = create.operator(\"galerkin\", estimatorModel, spc, fvspc)"
]
},
{
"cell_type": "markdown",
"id": "afd5c6c5",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
".. index:: Adaptation; Marking for h-Adaptation\n",
"\n",
"# Marker for grid adaptivity (h)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "166c8fd0",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:39:24.977148Z",
"iopub.status.busy": "2024-02-29T12:39:24.976897Z",
"iopub.status.idle": "2024-02-29T12:39:24.982782Z",
"shell.execute_reply": "2024-02-29T12:39:24.981780Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"hTol = 1e-16 # changed later\n",
"def markh(element):\n",
" center = element.geometry.referenceElement.center\n",
" eta = estimate.localFunction(element).evaluate(center)[0]\n",
" if eta > hTol and element.level < maxLevel:\n",
" return Marker.refine\n",
" elif eta < 0.01*hTol:\n",
" return Marker.coarsen\n",
" else:\n",
" return Marker.keep"
]
},
{
"cell_type": "markdown",
"id": "6364be88",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
".. index:: Adaptation; Marking for p-Adaptation\n",
"\n",
"# Marker for space adaptivity (p)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "2c5ed943",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:39:24.985866Z",
"iopub.status.busy": "2024-02-29T12:39:24.985626Z",
"iopub.status.idle": "2024-02-29T12:39:24.991627Z",
"shell.execute_reply": "2024-02-29T12:39:24.990597Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"pTol = 1e-16\n",
"def markp(element):\n",
" center = element.geometry.referenceElement.center\n",
" r = estimate.localFunction(element).evaluate(center)[0]\n",
" r_p1 = estimate_pm1.localFunction(element).evaluate(center)[0]\n",
" eta = abs(r-r_p1)\n",
" polorder = spc.localOrder(element)\n",
" if eta < pTol:\n",
" return polorder-1 if polorder > 1 else polorder\n",
" elif eta > 100.*pTol:\n",
" return polorder+1 if polorder < maxOrder else polorder\n",
" else:\n",
" return polorder"
]
},
{
"cell_type": "markdown",
"id": "d35261a1",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Operator for projecting into space with a reduced order on every element"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "45788c28",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:39:24.995896Z",
"iopub.status.busy": "2024-02-29T12:39:24.995655Z",
"iopub.status.idle": "2024-02-29T12:39:41.086120Z",
"shell.execute_reply": "2024-02-29T12:39:41.084823Z"
},
"lines_to_next_cell": 2
},
"outputs": [],
"source": [
"orderreduce = createOrderRedcution( spc )"
]
},
{
"cell_type": "markdown",
"id": "16cab1ee",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"# Main program\n",
"\n",
"Pre adapt the grid"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "434b6bac",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:39:41.091411Z",
"iopub.status.busy": "2024-02-29T12:39:41.091169Z",
"iopub.status.idle": "2024-02-29T12:41:29.930960Z",
"shell.execute_reply": "2024-02-29T12:41:29.929884Z"
},
"lines_to_next_cell": 2
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"pre adaptive ( 0 ): 240\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"pre adaptive ( 1 ): 231\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"pre adaptive ( 2 ): 363\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"final pre adaptive ( 2 ): 5.0 273\n"
]
}
],
"source": [
"hgrid = grid.hierarchicalGrid\n",
"hgrid.globalRefine(1)\n",
"for i in range(maxLevel):\n",
" print(\"pre adaptive (\",i,\"): \",grid.size(0),end=\"\\n\")\n",
" solution.interpolate( as_vector([p_0,s_0]) )\n",
" limit(solution)\n",
" step()\n",
" estimator(solution, estimate)\n",
" hgrid.mark(markh)\n",
" fem.adapt(persistentDF)\n",
"\n",
"print(\"final pre adaptive (\",i,\"): \",dt,grid.size(0),end=\"\\n\")"
]
},
{
"cell_type": "markdown",
"id": "57d6988c",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Define the constant for the h adaptivity"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "cf52e198",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:41:29.934992Z",
"iopub.status.busy": "2024-02-29T12:41:29.934776Z",
"iopub.status.idle": "2024-02-29T12:41:29.974409Z",
"shell.execute_reply": "2024-02-29T12:41:29.973075Z"
},
"lines_to_next_cell": 2
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Using timeTol = 3.216131221875019e-15\n"
]
}
],
"source": [
"solution.interpolate( as_vector([p_0,s_0]) )\n",
"limit(solution)\n",
"estimator(solution, estimate)\n",
"timeTol = sum(estimate.dofVector) / endTime\n",
"print('Using timeTol = ',timeTol, end='\\n')"
]
},
{
"cell_type": "markdown",
"id": "fe576240",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Time loop"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "2d209ccc",
"metadata": {
"execution": {
"iopub.execute_input": "2024-02-29T12:41:29.978432Z",
"iopub.status.busy": "2024-02-29T12:41:29.978239Z",
"iopub.status.idle": "2024-02-29T12:46:36.529262Z",
"shell.execute_reply": "2024-02-29T12:46:36.528108Z"
},
"lines_to_next_cell": 2
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"5.0 240 0.0 5.890350223214321e-17 # timestep\n"
]
},
{
"data": {
"image/jpeg": "/9j/4AAQSkZJRgABAQEAMgAyAAD/2wBDAAgGBgcGBQgHBwcJCQgKDBQNDAsLDBkSEw8UHRofHh0aHBwgJC4nICIsIxwcKDcpLDAxNDQ0Hyc5PTgyPC4zNDL/2wBDAQkJCQwLDBgNDRgyIRwhMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjL/wAARCACxAkUDASIAAhEBAxEB/8QAHwAAAQUBAQEBAQEAAAAAAAAAAAECAwQFBgcICQoL/8QAtRAAAgEDAwIEAwUFBAQAAAF9AQIDAAQRBRIhMUEGE1FhByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3ODk6Q0RFRkdISUpTVFVWV1hZWmNkZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6/8QAHwEAAwEBAQEBAQEBAQAAAAAAAAECAwQFBgcICQoL/8QAtREAAgECBAQDBAcFBAQAAQJ3AAECAxEEBSExBhJBUQdhcRMiMoEIFEKRobHBCSMzUvAVYnLRChYkNOEl8RcYGRomJygpKjU2Nzg5OkNERUZHSElKU1RVVldYWVpjZGVmZ2hpanN0dXZ3eHl6goOEhYaHiImKkpOUlZaXmJmaoqOkpaanqKmqsrO0tba3uLm6wsPExcbHyMnK0tPU1dbX2Nna4uPk5ebn6Onq8vP09fb3+Pn6/9oADAMBAAIRAxEAPwD3+iiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKK+Y7D4p+Op7m1WTX2KyOgYfZIOQSM/8s6idRQtfqdOHwlTEKTh0PpS+vbbTbKa9vJkhtoVLySOeFFV9K1qw1lJmspXZoWCSxywvDJGSMjcjgMMg5GRzWZ4u0+7ufCjRW4lvZoLi2uGTC751inSRlwABkqpAGBk4qro1w1x4q1bWvsd9BZXcVpZwedaSI8joZWZyhXcq/vFXcwA4PbFWcxtWviDTbzWrjSIJ3N9AheSNoJFG0HaSGKhW5IHBNIPEOmtqrabHLNLcq/luYraV442xna0iqUU4xwSDzWC+ooPifG32PU/KFg1oZxp05i80yqwHmbNuMD72dvvWFpenarpWqG2h/tr+0jrkkzsfM+xtZyTM7En/AFWdjH/bDYHQUAd/Prmm22t2ujS3QGoXSNJFAFYkquckkDA6HGSM4OM4NQWPibTNS1CSytPtkkkcskDyfYJxCHjJVx5pTZwVI+9yRxXGx6X4nt/G2j3t3pNnM8t7cSXF3DduwWMxlVBHlYQKmAoz8zZ5BYmrttAlv4mtB4eg1yEyX80uppdJOtqUfezkCX5MmQqR5f48ZoA7yiiigClq8d1LpVxHZlhOV+UI21iM8hW4wSMgHIwT1Fc3cWfiGWxtUhF9GYZZyENyu5suGg8xt+WjCkqwyWJHGeDXY0UAYmiWl7a6lqjTNfNazSCSH7ZOshDFnLBApIWMAptBweue1ed+LvjfN4S8U32hy+Go7h7UoDKmoFQwZFcceVxwwr1+vkr4y/8AJWtd+sH/AKTx13ZfQhWrck9rEydlod5/w0i3/QpD/wAGX/2qj/hpFv8AoUh/4Mv/ALVXhFFe1/ZeG7fiZ87Pd/8AhpFv+hSH/gy/+1Ug/aSYjP8AwiQ/8GX/ANqrwmkX7opf2Zhua1vxDnZ7x/w0i3/QpD/wZf8A2qj/AIaRb/oUh/4Mv/tVeEUU/wCy8N2/EOdnu/8Aw0kwH/IpD/wZf/aqP+GkmI/5FIf+DL/7VXg56GgdBS/szDc1rfiHOz3j/hpFv+hSH/gy/wDtVH/DSLf9CkP/AAZf/aq8Iop/2Xhu34hzs92P7STAZ/4RIf8Agy/+1Uv/AA0i3/QpD/wZf/aq8Hb7ppaX9mYbmtb8Q52e7/8ADSLf9CkP/Bl/9qo/4aRb/oUh/wCDL/7VXhFFP+y8N2/EOdnux/aSYf8AMpD/AMGX/wBqpf8AhpFv+hSH/gy/+1V4O3b60tJZZhrtW/EOdnu//DSLf9CkP/Bl/wDaqP8AhpFv+hSH/gy/+1V4RRT/ALLw3b8Q52e7f8NJNkD/AIRIf+DL/wC1Uv8Aw0i3/QpD/wAGX/2qvBz94UtJZZhm3p+Ic7Pdv+Gkmzj/AIRIf+DL/wC1UH9pMg/8ikP/AAZf/aq8HOd3GOlIQS3bpWUsuoLaLHzs94/4aUP/AEKY/wDBl/8AaqP+GlDz/wAUmP8AwZf/AGqvBWGP7vSg/wAX3ayeBop6r+tQ5mfTuifGS51zSLu/h0GxiFu4QxTasVZyQOn7n3q9P8TNXgNxnw7preRt3bdbHO7pj91Xy3FeXVujpDcyRqxBKo5AP5VM2rai2/OoXB3Y3Zmbn9a8+vlVWVSUqVXlj2cb2+d/L8RxrSjo4J/Nn06/xO1hJzCfDenbg6Jka2uMsMj/AJZUf8LO1j/oW9O/1/kf8htfvf8AfrpXzCdW1EuWN/cE7gc+c3UdO9H9raj/AM/9x9/d/rm6+vWsv7IxP/P9f+A/8Ev6w/8An3H75H09/wALO1j/AKFvTv8AX+R/yG1+9/366Uf8LO1j/oW9O/1/kf8AIbX73/frpXzD/a2o/wDP/cff3f65uvr1o/tbUf8An/uPv7v9c3X160f2Rif+f6/8B/4IfWH/AM+4/fI+nk+J2sPOIR4c04MXdMnWxjKjJ/5ZURfE7WJRblfDenfv923/AInY429c/uq+Yhq2ohwwv7gHcTnzm6nr3oXVtSXZt1C5G3O3EzcfrR/ZGKf/AC/X/gP/AAfMPrD/AOfcfvke53P7Rb2txJA/hRC6OyNt1PIyDj/nlUP/AA0oeP8Aikxz/wBRL/7VXg7O0jbnYMxJJJOSTTR/D92vW+pUHbTt1fkZ8zPev+GlD/0KY/8ABl/9qpf+Gkmxn/hEh/4Mv/tVeChTtzgdKd823t0qoZfRa96L2DmZ7z/w0kxH/IpD/wAGX/2qj/hpJgP+RSH/AIMv/tVeDjoKD0Nb/wBmYblvb8Rc7PeP+GkmI/5FIf8Agy/+1Uf8NIt/0KQ/8GX/ANqrwcdBS01leGtt+Ic7Pd/+GkW/6FIf+DL/AO1Uh/aSYDP/AAiQ/wDBl/8Aaq8JpG+6aUsswyi3b8Q52e8f8NIt/wBCkP8AwZf/AGqj/hpFv+hSH/gy/wDtVeEUU/7Lw3b8Q52e7/8ADSLf9CkP/Bl/9qpD+0kw/wCZSH/gy/8AtVeE0jdvrSllmGSvb8Q52e8f8NIt/wBCkP8AwZf/AGqj/hpFv+hSH/gy/wDtVeEUU/7Lw3b8Q52e7/8ADSLf9CkP/Bl/9qpP+GkmyB/wiQ/8GX/2qvCaQ/eFKWWYZLb8Q52e8f8ADSLf9CkP/Bl/9qo/4aRb/oUh/wCDL/7VXhFFP+y8N2/EOdnu/wDw0i3/AEKQ/wDBl/8AaqT/AIaSbOP+ESH/AIMv/tVeE0n8f4UpZZhlbT8Q52e8f8NIt/0KQ/8ABl/9qo/4aRb/AKFIf+DL/wC1V4RRT/svDdvxDnZ7v/w0i3/QpD/wZf8A2qvQvhx8Qj8QLO/uDpQsBaSKm37R5u/IJz91cdPevkavoT9nL/kDa7/18R/+gmuLH4GjRo88FrcqMm2e2UUUV4hoFfH2lw4lsmx/Eh/lX2DXyrp9pi3tXAwdqEHHsK4MbPlcPU+n4dhGUK9/5f8AM+qJJI4YnlldUjQFmdjgKB1JPYVXsNSsNVtvtOnXtteQZK+bbyrIuR1GVJGawfG8M7+C5VmAnCT20l2IoyA8Kzo0uFyTjYGyMnjNVNA1LTrnxhreo2N3bHTLmOzt450kHlz3I83cEPRm2mMHGemO1d58wdPDqunXN/NYQX9rLeQjdLbpMrSRjOMsoOR+NSy31pDdwWkt1BHc3AYwwvIA8m0ZbavU4HXHSuM/tnQoPi0tsupafHcNpzQNEJ0DmYzKdhGc7yOcdax5JtS/4WJod/qWg6gl3NfzxRuWgKJAIZFRUxJnABMjZAOSQM4UUAeivrOlx6mmmPqVmuoOMratOolYYzwmcnj2ptrrek3t7LZWmqWVxdxZ8yCK4R3THXKg5FecaBe32l6kLU6msurT67Mt1pJgTcYWlc+eTjzOE2sHztwAuK17XXNC8SeNbGGy1HTYIdHuJRDEsyLPdTlGRgiZ3CNQzEnHzEA9FyQDv6KKKACiqWr3klhpVxdRBS8a5BcZVecFm9h1PTgVzdx4uuo7G1lhNjITLOs0oJ2uI3AVY13f6yRSGVcnI9Qc0AdjXyV8Zf8AkrWu/WD/ANJ46+mNE1m41LUL6CbyCsJOFiUhoj5kibXyTk4QHt16dDXzH8YJVn+K2uSIHALQjDoUPEEY6EA16eU/7x8mRPY4iiiivpjEKRfuilpF+6Kl/Ev67ALRRRVAIehoHQUHoaB0FT9oBaKKKoBG+6aWkb7ppalfE/67gFFFFUAjdvrS0jdvrS1K+JgFFFFUA1uowKXJ/u/rQfvClrNJtvX+rDG5O7p29aOd2cfrS/x/hS0owbvr18gGNknG3tSH+L5acc7uBnimkFifl/WsKkd7b/8AAfkMD/F8tB/i+Wjbkkbf1o25JG39azal/Xz8gA/e+73FH/Af4qNuWPy/rQFz/D39aLSv/X+QB/wH+Kj/AID/ABUBc/w9/WgLn+Hv60Wl/X/DAA+993uaB/D8tAXn7vf1oC8A7f1oSlf+vLyAB/D8tABwDtFAXgHb+tKpO0fL+tVTjdq/5enkDAZ2Yx+tLk7enb1pV+6KD0Nbxh7id+ghATgfL+tBJweP1pR0FB6Gq5Xy7/kADoKWkHQUtaR2QgpG+6aWkb7ppT+FgLRRRVAFI3b60tI3b61M/hAWiiiqAKQ/eFLSH7wqZ7fd+YC0UUVQBSfx/hS0n8f4VMugC0UUVQBX0J+zl/yBtd/6+I//AEE18919Cfs5f8gbXf8Ar4j/APQTXmZt/u/zRcNz2yiiivmTYK+brG2xpls2P+WKn9BX0jXgtrbEeH4XA5FqpGf9yvEzifLKl6/5Hv5HU5I1vNf5nvVFcv43muE8FymYi3Dz20d2YZCQkLTosuGwDjYWycDjNQ+Gra103xhr+naVFHDpkUFpIIIQBHFO3m7woHAJVYiQPUHvXtngHXUVx9vp9lZfFaSW2t44pbnSXkmdRzI3nLyTVdNM0688a28mgWUMD2Fy8uqajEuDKzIw+zlushywZgchdo7kYAO4orhJvD+i3XxLsVstIsIH0yNtRu7iG2RXeaTckSswGT/y1c5PUKaxWsr7SPEGiD+zFGrTa1J52qrOhN3bsZCy4B3kKmMqwCrsGD0oA9VooooAKKKKACvkn4zE/wDC2tdwO8H/AKTx19bV8lfGX/krWu/WD/0njr0sqV8R8mRPY4PJ/u/rRk/3f1paK+l5X3/IxEyf7v60ik7Rx+tOpF+6KhxfMtfy8hhk/wB39aMn+7+tLRV8r7/kIaScH5f1oBOB8v60p6GgdBUcr5t/yGGT/d/WjJ/u/rS0VfK+/wCQhrE7Tx+tLk/3f1ob7ppahRfM9fy8xiZP939aMn+7+tLRV8r7/kIaxPHHf1pcn+7+tDdvrS1Ci+Z6/kMTJ/u/rRk/3f1paKvlff8AIQ0k7hx+tLk/3f1oP3hS1EYu71/LsMbk7unb1pcn+7+tH8f4UtEYvXX8gG5O7p29aBnJOP1pf4/wpaIwb69fIBoJ3Hj9aBnJOP1pR940tEYNrfq+wDQTuPH60DIzx+tKPvGlojB7379gGqTzx39aBkZ4/WlXv9aWiEHZO/5ANUnnjv60DIGMfrSr3+tLRCDsnf8AIBqk7Rx+tAyBjH60q/dFLRCD5U79PIBqk7Rx+tBJwfl/WlX7ooPQ0oxfItenkHUQE4Hy/rQScH5f1pR0FB6Gnyvl3/IBATgfL+tLk/3f1oHQUtOMXZa/kAmT/d/WkYnaeP1p1I33TSnF8r1/IAyf7v60ZP8Ad/Wloq+V9/yEJk/3f1pGJ447+tOpG7fWonF8u/5DDJ/u/rRk/wB39aWir5X3/IQmT/d/WkJO4cfrTqQ/eFROLtv27dxhk/3f1oyf7v60tFXyvv8AkITJ/u/rSZO7p29adSfx/hUSi9NfyGGT/d/WjJ/u/rS0VfK+/wCQhMn+7+tfQn7OP/IF13I/5eY//QTXz5X0J+zl/wAgbXf+viP/ANBNebmqaw+/VFw3PbKKKK+aNgrxy3tseFYmx/y4g/8Ajlex15nHb7fBCPtzjTgcD/rnXzHEU+WdD/F/kehgKnIp+aPSpI45onilRXjdSrIwyGB6gjuKpW2iaTZ20dta6XZQQRS+dHFFboqpJ/fAAwG9+tZfi7ULu28KNLbmWynnuLa2ZyV3wLLOkbNkEjIViQcnBxUXh/zbDxVrOirdXVxZQW9rcxfaZ2meNpDKrLvcliP3atgk43HtX0555pt4Y0BtT/tNtD0w6hvEn2o2kfm7x0bfjOffNJF4W8PQX4v4dB0uO9DmT7QlnGJN56tuAznnrWVbxS2nxLliW9vZILjTHnaGW4Z41fzVAKoTheOOBWFFeX32e28RHUbw3c3iFrB7Uzt5Ig+1Nb+X5WdoIUB92N2R1xxQB6HHa28M808UESTTkGWRUAaQgYG498Djmq1vo2l2d/NfWum2cF5Pky3EUCrJJk5O5gMnn1rz6TUtV03WNOkuP7a/tWfWzbXCv5n2NrV3cJsB/d8JsI2/PkHdxmpNF1m/uNQj1rVbe8+zz6vLYROmpOqxESvEim2ACFflALEltxJxjFAHpdFFFAFe+vI7CzkuZVZlQD5UHzMScADPckgfjWbL4lght7WdrK+Mc8zwuyxqVt2V9jeY27CjcDyCRxWrc20N5bSW86b4pBhhkj8iOQfcVR/sDTvs4t9k5hEhlMZupSrMSCdw3fMCRkg5BJJ6k5AHaZrEepz3cAtbq2ltZNjJcoELjJAdRkkodpwe+PrXy38Zf+Sta79YP/SeOvqu1062s55p4hIZZsb3kleQkDJABYnAG5sAccmvlD4wQxW/xW1yOGJI0DQkKigDJgjJ4HuSa9PKf94+TInscRRRRX0xiFIv3RS0i/dFS/iX9dgFoooqgEPQ0DoKD0NA6Cp+0AtFFFUAjfdNLSN900tSvif9dwCiiiqARu31paRu31palfEwCiiiqAaxwRS7vY0H7wpazSd3Z/1YY3d83Q9KXd7Gj+P8KWlFS11Abu+boelLu9jR/H+FLRFS11AaG+Y8Gl3exoH3jS0QUrb9/wAwGhvmPBpd3saB940tEFK2/f8AMBqt14PWl3exoXv9aWiCly7gNVuvB60u72NC9/rS0QUuXcBqt8o4NLu9jQv3RS0QUuVagNVvlHBoLcHg0q/dFB6GlFS5Fr0DqIG4HBoLcHg0o6Cg9DTtLl3AB0FLSDoKWtI7IQUjfdNLSN900p/CwFoooqgCkbt9aWkbt9amfwgLRRRVAFIfvClpD94VM9vu/MBaKKKoApP4/wAKWk/j/Cpl0AWiiiqAK+hP2cv+QNrv/XxH/wCgmvnuvoT9nL/kDa7/ANfEf/oJrzM2/wB3+aLhue2UUUV8ybBXn0ZH/CBr/wBgwf8AoqvQa84jfHgVf+waOn/XKvnM/o+0lQ8pf5ESq+zt5nf31lbalZTWV5Ck1tMpSSNxwwqlY+HdN01GW1S4RmlWZ5GupXkkYDA3uzFnAHG1iR7VW8TaxPY+Gjd6f8lxPNBbQtNER5bSypEGKtg8b84PpUWhXd/D4i1XQ76+e/FtBb3UNzLGiSFZTIpVtgVeDESCAOG9s19GWWW8KaU2ujWj9u+3jo41G4C4znbs37NuQDtxj2p6+F9HXV/7UFoftXmmb/Wv5YkI2mQR52B8cbsZ965j+3dZ2f299vP2L+3P7N/s7yk2eV9o+zbt2N+/d8/3sY4x3qm3jG8tdYtjcat/pcutf2fLonkp+6gaUxxy527+RsfeW2ndgCgDto/Dumx6t/aZimkugzOhmuZZEiZsglEZiqEgkfKBwTUa+FtGXVf7SFoftHmmfHnP5QlIwZBFu2B/9rbn3qjaXGpQ/EC60+41KS5spLD7VHA0UaiE+aVwCBuPHqTVE+IdSvPHOlx2k4TQ5ZLi227FJuXjQlnDEZCqw2jGMlWPIxQB2lFFFABRTJZY4InlmkWONFLO7nAUDqST0FVzqunKlu5v7ULcnEBMy4lPovPzfhQBbr5J+MxA+LWu59YP/SeOvq+C8tbqSWO3uYZXhbbKsbhijejAdDwetfKPxl/5K1rv1g/9J469LKr/AFjTsyJ7HB7h60bh60tFfS2l3/r7zETcPWkVhtHNOpF+6KhqXMtf608xhuHrRuHrS0Vdpd/6+8Q0sMHmgMMDmlPQ0DoKi0ubf+vvGG4etG4etLRV2l3/AK+8Q1mG080u4etDfdNLUJS5nr/WvmMTcPWjcPWloq7S7/194hrMOOe9LuHrQ3b60tQlLmev9feMTcPWjcPWloq7S7/194hpYbhzS7h60H7wpaiKld6/1b1GN3Dd+FLuHrS0VSjJdfw/4IDdw3de1LuHrR/H+FLUxUtdf6+8BoYbjS7h60tFUoyXX+vvAaGG480u4etA+8aWpgpW37/n6gNDDn60u4etLRVKMkrX/D/ggNVhzz3pdw9aF7/WlqYKXLv/AF94DVYBaXcPWloqlGSVr/h/wQGqw2jmgsMHmlX7opamKk4LX+vvAaGGBzQWGDzTqQ9DTako7/194CBhgc0u4etA6CloipWWv9feAm4etIzDaeadSN900pqXK9f6+8A3D1o3D1paKu0u/wDX3iE3D1pGYcc96dSN2+tRNS5d/wCvvGG4etG4etLRV2l3/r7xCbh60hYbhzTqQ/eFRNStv2/P1GG4etG4etLRV2l3/r7xCbh60m4buvanUn8f4VElLTX+vvGG4etG4etLRV2l3/r7xCbh619Cfs4nOi67j/n5j/8AQTXz5X0J+zl/yBtd/wCviP8A9BNebmvN9X17ouG57ZRRRXzRsFeZRN/xRSD/AKhw/wDRdem15ZC3/FGxj/qHj/0XXDjaPtHDyZ5WZVfZun5v/I9H1bS7bWdMm0+7D+TKBko21lIIZWU9iCAQfUVS07w6unSzXA1K+nvLiRGuLmbyy8qoCFjICBQoyfugHknPJqxrurpoejzX7wvMVZI44UIBkkdwiKCemWYDPbNV9G1q4vr++03ULJLPULNY5HSKfzo2jk3bWViqk8o4IKjBXvXceqV/+EPsf7T+1farz7N9r+3fYN6+R9o6+Z93dnd82N23dzjNWJPDqXOqx3t5qN9dRwzefBaSmMQxPzhhtQM2MnG5mx16gVm/8JjL9s8z+zB/Y/8AaP8AZn237R8/nb/Lz5e3GzzPkzuznnGKlv8AX9a0/XtOsJNIsJIb+7aCJ4tQcyiNVZmkZDCAMKvI3dSBk5zQBP8A8Iw//CT/ANu/25qfm7PK+z4g8ryt27Z/qt2M99273qFvAPhj+2LHVIdHs4LmzlaZDDAih3I4LfLkkHkcjB5qK78Xtp/im10a6h08tdT+VFHb3/mXQBBId4Ng2pxydxxV201y8l8WXOiXWmxwIlv9phnW53mRN+zldo2+vU0AbtFFFAFTUrL+0NPltfM8svgq+3IDAgjI7jIGR3rIk8NTTWUdvLdWjMtzJcCf7GfMjLvvbyzv+Q7i2Dzxjg4JPRUUAZOkaK2l3NxK1yJVdRHGoj2lUDyOMnJ3HMh546dOtfLnxgiWD4ra5GhcgNCcu5c8wRnqSTX1zXyV8Zf+Sta79YP/AEnjr08p/wB4+TInscJRRRX0xiFIv3RS0i/dFS/iX9dgFoooqgEPQ0DoKD0NA6Cp+0AtFFFUAjfdNLSN900tSvif9dwCiiiqARu31paRu31palfEwCiiiqAQ/eFLSH7wpamO7/roAUUUVQCfx/hS0n8f4UtTHqAUUUVQCD7xpaQfeNLUw2+/8wCiiiqARe/1paRe/wBaWph8IBRRRVAIv3RS0i/dFLUw+FAFIehpaQ9DTlsADoKWkHQUtEdkAUjfdNLSN900p/CwFoooqgCkbt9aWkbt9amfwgLRRRVAFIfvClpD94VM9vu/MBaKKKoApP4/wpaT+P8ACpl0AWiiiqAK+hP2cv8AkDa7/wBfEf8A6Ca+e6+hP2cv+QNrv/XxH/6Ca8zNv93+aLhue2UUUV8ybBXksJB8HxqeQbAAg/8AXOvWq8YtroHwvEnHNkB/45W9Cg6t/I+Z4hm4ToW/m/yPTtf0Ial4ebTrExWskckU9t8mESSKRZEyB/DuUA47E1V0rTdYh1i91m/isftd4ILcwQ3DskMEe85DlAXYtIxxtUdBnjJ1tW1S20bTJ9QuywhiAyEXczEkBVUdySQAPU1X0jXYtWlurdrS6sry1K+dbXQXeoYZVsozKQcHkE9CO1YH0xgf8Irqvmf2Z5tn/Yn9rf2p5u9vP/13n+Vs27ceZ/Fu+7xjvW0ukXEvjF9ZuXiMEFmLayjUkspZt0rNkYBO2MDGeFPrUH/CYWP9p/Zfst59n+1/Yft+xfI+0dPLzu3Z3fLnbt3cZzUth4kbUtTuLS20bUTBb3L20l6WgEQdevHmbyOn8PegCtfWHiHVNStElXTrSztb1bgXMNw7zSIpJCbCgC7hgMdx4J45pn2DxB/wnZ1b7Jpn9nm2+yZ+2yebs37t+3ycZ7bd341fvPE1jZ+JtP0BlmkvL0MQ0agpFhHcbznjcI3xgH7p6VWTxjZPqIt/sd8LZrxrBb4onkmcEgp97f8AeBXO3bkdaAOiooooAKKpavHdS6VcR2ZYTlflCNtYjPIVuMEjIByME9RXPz2viSTS7OK2jkjaG+8yRZrzbI8ImVkXeu7cojJDAnJK/wAQPzAHW18lfGX/AJK1rv1g/wDSeOvprQrS9tL7UvOa+a1lkEkX22dZHDFnLBdpIWMDZtBweue1fMPxgSQfFbXBO6PJuhyyKVGPIjxwSe2O9ellV/rGnZkT2OIopNo9KNo9K+lvLt/X3GItIv3RRtHpSKo2jioblzLT+tPIY6ik2j0o2j0q7y7f19wgPQ0DoKQqMHigKMDiovLm2/r7hjqKTaPSjaPSrvLt/X3CBvumlprKNp4pdo9KhOXM9P618hi0Um0elG0elXeXb+vuEDdvrS01lHHHel2j0qE5cz0/r7hi0Um0elG0elXeXb+vuEB+8KWmlRuHFLtHpURcrvT+regxaKTaPSjaPSrvLt/X3CD+P8KWm7Ru6dqXaPSoi5a6f19wxaKTaPSjaPSrvLt/X3CAfeNLTQo3Hil2j0qIOVtu/wCfoMWik2j0o2j0q7y7f19wgXv9aWmqo5470u0elRBy5dv6+4YtFJtHpRtHpV3l2/r7hAv3RS01VG0cUu0elRBy5Vp/X3DFpD0NG0elIVGDxTblbb+vuAUdBS00KMDil2j0oi5WWn9fcAtI33TRtHpSMo2nilNy5Xp/X3AOopNo9KNo9Ku8u39fcIWkbt9aNo9KRlHHHeom5cu39fcMdRSbR6UbR6Vd5dv6+4QtIfvCjaPSkKjcOKiblbbt+foMdRSbR6UbR6Vd5dv6+4QtJ/H+FG0elJtG7p2qJOWmn9fcMdRSbR6UbR6Vd5dv6+4QtfQn7OX/ACBtd/6+I/8A0E189bR6V9Cfs4jGi67j/n5j/wDQTXm5rzfV9e6Lhue20UUV80bBXz9aXO7Q4EyebZR1/wBmvoGvmawus2Fqmf8Alkg/QV7WUQ541fT/ADPn88pc8qPk/wDI978TaPPfeGjaaeN9xBNBcwpNKT5jRSpKFLNk87MZPrVXSIdSbxFqOuXelXFsLyO2s47YyxNIiRmRjI+1yuMykYBJwOnOK6C+vbbTbKa9vJkhtoVLySOeFFV9K1qw1lJmspXZoWCSxywvDJGSMjcjgMMg5GRzXin0ByP9haxs/sH7AfsX9uf2l/aPmps8r7T9p27c79+75Pu4xzntU11oklz4ltbnTvDP9mXMWoedcatvhXz4gTuHyMXfeOMOoAznqK6H/hKNH/tf+y/tZ+1eb5P+qfy/Mxu8vzMbN+Oduc+1L/wkmmDVF053uIp3kMSNLaSxxO4z8qyMoRjweATnFAHKyeFfEtv4n0q+jvtPuYRqk13cTG0ZZFVoXQBiZfmAUrGu0DHBIODmvF4a1hfEEVw1hcm7XWGu3vmu0NmYC54EG7iTyiF3bM7ud5rtB4h01tVbTY5ZpblX8tzFbSvHG2M7WkVSinGOCQeaZZ+JtM1DUpLG1+2SSxyyQO4sJxCroSGHmlNnBBH3uvAoA16KKKACiiigAr5K+Mv/ACVrXfrB/wCk8dfWtfJXxl/5K1rv1g/9J469PKf94+TInscJRRRX0xiFIv3RS0i/dFS/iX9dgFoooqgEPQ0DoKD0NA6Cp+0AtFFFUAjfdNLSN900tSvif9dwCiiiqARu31paRu31palfEwCiiiqAQ/eFLSH7wpamO7/roAUUUVQCfx/hS0n8f4UtTHqAUUUVQCD7xpaQfeNLUw2+/wDMAoooqgEXv9aWkXv9aWph8IBRRRVAIv3RS0i/dFLUw+FAFIehpaQ9DTlsADoKWkHQUtEdkAUjfdNLSN900p/CwFoooqgCkbt9aWkbt9amfwgLRRRVAFIfvClpD94VM9vu/MBaKKKoApP4/wAKWk/j/Cpl0AWiiiqAK+hP2cv+QNrv/XxH/wCgmvnuvoT9nL/kDa7/ANfEf/oJrzM2/wB3+aLhue2UUUV8ybBXyZp1zuS0TdjIQcfhX1nXxnpN1m7sUz/HGP1FfT8OU+eFfyj/AJnn46lzuHkz6p8Xafd3PhRorcS3s0FxbXDJhd86xTpIy4AAyVUgDAycVV0a4a48VatrX2O+gsruK0s4POtJEeR0MrM5QruVf3iruYAcHtiurkkjhieWV1SNAWZ2OAoHUk9hVew1Kw1W2+06de215Bkr5tvKsi5HUZUkZr5g9A4L7HfeR/wj39nXn2z/AISL7f8AavIbyfI+1faPM83G3Oz5Nud2e2OavayINU8T6ZLZ6drE1/a36eYLmO5S0jjXIaQbsREhclSuSTiut/tfTf7T/sz+0LT+0Nu77L5y+bjGc7M5xj2oh1XTrm/msIL+1lvIRult0mVpIxnGWUHI/GgDzzS9O1XStUNtD/bX9pHXJJnY+Z9jazkmZ2JP+qzsY/7YbA6CtOKBLfxNajw9BrkMsmoSS6ktwk62hjbeXI8z93ksQVMfP4ZrrhrOltqZ0xdSszqAGTaidfNAxn7mc9Pakh13SLjUn02DVbGW+jJD2yXCNKpHXKg5GPpQBfooooApaveSWGlXF1EFLxrkFxlV5wWb2HU9OBXPz+LLiPS7OeFI7iR777PNLDC8kaoJlTOFJKs6srKCcYOeeAetooAwtC1m41O9vYZvIxCAQIlIaIl5F8t8k5YBAe33unTPzH8YJGk+K2uM8TxHdCNjkEj9xH6Ejnr1719ZXN5a2SI91cwwK7BFMrhQzHsM9T7VwGu/C7wn4m8W397qtzI2pXWyQQQ3IVhGsaIDt69VPPvXZgcRHD1eeW1iZK6PlWivqL/hQvgr+5qH/gT/APWo/wCFC+Cv7mof+BP/ANavX/tij2f4f5mfs2fLtIv3RX1H/wAKF8Ff3NQ/8Cf/AK1H/ChfBQ/g1D/wJ/8ArVP9r0b3s/w/zH7Nny7RX1F/woXwV/c1D/wJ/wDrVRk+DHgAh5Y72QQ2shF2xvgRGACME/wndjr6Gq/tij2f4f5i9mz5qPQ0DoK+oU+BHgiWNXT7e6MMqy3WQQe4OKd/woXwV/c1D/wJ/wDrVP8Aa9G97P8AD/Mfs2fLtFfUX/ChfBX9zUP/AAJ/+tR/woXwV/c1D/wJ/wDrVX9sUez/AA/zF7Nny433TS19Rf8AChfBR/g1D/wJ/wDrVBc/BHwDZIj3U11ArsEUy3oUMx7DI5PtU/2vRvez/D/Mfs2fMlFfS8HwS8Cm5eyluZ3vVLOYUvBvCZ+UlcZ6Fcn1q3/woXwV/c1D/wACf/rVX9sUez/D/MXs2fLjdvrS19Rf8KF8Ff3NQ/8AAn/61H/ChfBX9zUP/An/AOtUrN6N27P8P8x+zZ8u0V9Rf8KF8Ff3NQ/8Cf8A61H/AAobwV/d1D/wJ/8ArVX9sUez/D/MXs2fLh+8KWvpGX4OfDw20V9HqDCyWQiWZr9dhG08BugOdv4ZrQX4D+CWUMovypGQRddf0qVm9FN6P8P8x+zZ8vUV9Rf8KF8Ff3NQ/wDAn/61H/ChfBX9zUP/AAJ/+tVf2xR7P8P8xezZ8ufx/hS19Rf8KF8FZzs1D/wJ/wDrUf8AChfBX9zUP/An/wCtUrN6K6P8P8x+zZ8u0V9NzfBHwDbzQwzTXUcsxIiR70BpCOu0Ec9R0pLb4G+CWItZZbqS9iiRp0ju+RnIDbcZAJVsZ9D6VX9sUez/AA/zF7NnzIPvGlr6i/4UL4K/uah/4E//AFqP+FC+Cv7mof8AgT/9apjm9FLZ/h/mP2bPl2ivqL/hQvgr+5qH/gT/APWo/wCFC+Cv7mof+BP/ANaq/tij2f4f5i9mz5cXv9aWvqH/AIUL4KH8Gof+BP8A9as9vg78PBHHejUWGnndG0pv12mTjaA3Tpu4+lTHN6KVrP8AD/Mfs2fN1FfUI+A3gojIXUCD/wBPP/1qX/hQvgr+5qH/AIE//Wqv7Yo9n+H+YvZs+XF+6KWvqL/hQvgofwah/wCBP/1qP+FC+Cv7mof+BP8A9apjm9FJKz/D/Mfs2fLtIehr6j/4UL4K/uah/wCBP/1qgb4I+ABdizaa6Fyy7xCb0byvrtxnHB5pvOKLWz/D/MPZs+Yx0FLX0xY/A/wPcQBRPcTzxYScwXYIWQDkYxxz2NWv+FC+Cv7mof8AgT/9ahZvRS2f4f5h7Nny7SN9019R/wDChfBX9zUP/An/AOtR/wAKF8FH+DUP/An/AOtSlm9Fpqz/AA/zD2bPl2ivqL/hQvgr+5qH/gT/APWpkvwK8DQRPLMb6ONFLO73eAoHUkkcCq/tij2f4f5i9mz5gpG7fWvpJ/g58PYNtzPqDJZzKogdr9QHbJ3YY8Efd6e9aP8AwoXwV/c1D/wJ/wDrVMs3otWs/wAP8x+zZ8u0V9Rf8KF8Ff3NQ/8AAn/61H/ChfBX9zUP/An/AOtVf2xR7P8AD/MXs2fLtIfvCvqP/hQvgr+5qH/gT/8AWo/4UL4K/uah/wCBP/1qmWb0Wtn+H+Y/Zs+XaK+mLj4IeBvMeyhuLhNQeItFFJeAt3w23GSMj9KdbfBP4f3vmfZLi5n8ttj+VehtrehwODVf2xR7P8P8xezZ8y0n8f4V9R/8KF8Ff3NQ/wDAn/61H/ChfBWc7NQ/8Cf/AK1S83ovo/w/zH7Nny7RX1F/woXwV/c1D/wJ/wDrUf8AChfBX9zUP/An/wCtVf2xR7P8P8xezZ8u19Cfs5f8gbXf+viP/wBBNbNz8D/AdnA091JeQQrjdJLeBVGTgZJHrXReEfC2g+B7me00u+UR6isckcE06s7su7LL3IIK9PSuTG5hTr0uSKdyowaZ2NFFFeOaBXxJpJP9o2IBGfOjxn6ivtuvh/R2/wCJrYD/AKbx/wDoQr6bh2t7OFfzj/mY1lex9ceN4Z38FyrMBOEntpLsRRkB4VnRpcLknGwNkZPGaqaBqWnXPjDW9Rsbu2OmXMdnbxzpIPLnuR5u4IejNtMYOM9Mdq7OivmTY8s8y32f2Plf+Ej/AOEp+1eVj975X2nf5uOuz7P8u7pj5fatn+2dCg+LS2y6lp8dw2nNA0QnQOZjMp2EZzvI5x1ruqKAPI4bgHXoLYz2gP8Awksk/wDY6oftqEyMPOL5/wBXz5mNuNpxvxxW7rN9pmoeLdF+x629/PBqIDaRbtFiFgrK0rbU8xduSTubB6dxXf0UAFFFFABRRRQBk67oraxCES5EB8uWFy0e/KSLtbAyMN0wecehzTItEmGqpey3aMPNFy6JCVJm8nySQdxwm3nbgnP8VbNFABRRRQAUUUUAFYFp4eurMXaxX8JWa2FtHutiSiqzlCfn+YjzGz0zx0wc79FAFbTraSy022tZZVleGJYy6psDYGM4ycfmas0UUAFFFFABWTruitrEIRLkQHy5YXLR78pIu1sDIw3TB5x6HNa1FAGXFp16dZjvru9gmjjh8uOFLcpsYgb2B3nqR3BwOAepOpRRQAUUUUAFRXMC3VrNbuWCSoyMVOCARjipaKAOfj0C9j069tTqMBa7ZGZ1t5I9uEVDjbKCMhE6Ec7uoOBtWkH2Wzgt9+/yo1TdtAzgYzgcD8KmooAKKKKACiiigDG1bRJtQu0uILtIGAj3CSEyAmOQSJjDDHzDn1Hp1p2n6F/Z+uX+pLezy/bFUNFIq/KQzHO4AEgBgoB6AdTxjXooAKKKKACiiigCK5gW6tZrdywSVGRipwQCMcVgf8IzcPp91az3NhcG4dJMzWBZUZYxHkL5nXCrjng7vUAdJRQBHbxeRbRQmR5DGgXe5yzYGMn3qSiigAooooAKxp9Emk1uPUIrtERZlmaJ4SxLBGjOG3DAKMeMcHnnkHZooAxtA0AaFHJGk6vGURFVY9nC5+ZuTuc55bjO0cDFbNFFABRRRQAVU1Ky/tDT5bXzPLL4KvtyAwIIyO4yBkd6t0UAYsOhzCxjtLi7SSMXpu5NkJXflzLt5Y4HmHOf7o2+praoooAKKKKACiiigDGn0SaTW49Qiu0RFmWZonhLEsEaM4bcMAox4xweeeQU8P6A+iefvu/tHmhFHyEbQufVmx1+6u1R2UZOdqigAooooAKKKKAKWpWMl7DD5MyxTQyrNGzpvXIyOVyMjBPcc4NZKeFsR6XG90pWxWJGZI3RpVibdGDh9vGBnKnq2NueOjooAKKKKACvnu0/Z81u0u4JxrWnv5Lq4Uq43YOcZxxX0JRW1HEVKN/Zu19xNJ7nL+N5rhPBcpmItw89tHdmGQkJC06LLhsA42FsnA4zUPhq2tdN8Ya/p2lRRw6ZFBaSCCEARxTt5u8KBwCVWIkD1B711kkcc0TxSorxupVkYZDA9QR3FUrbRNJs7aO2tdLsoIIpfOjiit0VUk/vgAYDe/WsRnnflQbP7Ywv/CR/8JT9l83P73yvtOzys9dn2f5tvTHze9HlQbP7Ywv/AAkf/CU/ZfNz+98r7Ts8rPXZ9n+bb0x83vXo/wDZGm/2p/af9n2n9obdv2ryV83GMY34zjHvR/ZGm/2p/af9n2n9obdv2ryV83GMY34zjHvQB5x5UGz+2ML/AMJH/wAJT9l83P73yvtOzys9dn2f5tvTHze9HlQbP7Ywv/CR/wDCU/ZfNz+98r7Ts8rPXZ9n+bb0x83vXo/9kab/AGp/af8AZ9p/aG3b9q8lfNxjGN+M4x70f2Rpv9qf2n/Z9p/aG3b9q8lfNxjGN+M4x70AXKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigD/9k=",
"text/plain": [
"