[![binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gl/rubendibattista%2Fjosiepy/master/?filepath=examples%2F004_Two-phase-sigma.ipynb) # Multiphase flow including capillarity effects from Cordesse This notebook (partially) implements the numerical strategies and the model presented in the thesis by Pierre Cordesse, [*Contribution to the study of combustion instabilities in cryotechnic rocket engines : coupling diffuse interface models with kinetic-based moment methods for primary atomization simulations*](http://theses.fr/2020UPASC016). In that work, the author proposes a series of submodles that are solved sequentially with a splitting approach. This notebook for the moment implements the homogeneous Euler part with the instantaneous relaxation of the phase pressures. ## Define fields ```python from josie.fluid.fields import FluidFields from josie.fields import Fields class SigmaFields(FluidFields): # Conservative rho = 0 rhoU = 1 rhoV = 2 rhoY = 3 rhoYw = 4 rhoa = 5 rhoz = 6 # Auxiliary U = 7 V = 8 P = 9 c = 10 Y = 11 w = 12 alpha = 13 z = 14 p1 = 15 p2 = 16 class SigmaConsFields(Fields): rho = 0 rhoU = 1 rhoV = 2 rhoY = 3 rhoYw = 4 rhoa = 5 rhoz = 6 ``` ## Define State ```python from josie.fluid.state import SingleFluidState from josie.state import SubsetState class SigmaConsState(SubsetState): full_state_fields = SigmaFields fields = SigmaConsFields class SigmaState(SingleFluidState): fields = SigmaFields cons_state = SigmaConsState ``` ## Define Barotropic EOS ```python class Barotropic: """ .. math:: \pressure = a\rho^b """ def __init__(self, a, b): self._constant = a self._exp = b def p(self, rho): return self._constant * rho ** self._exp def rho(self, p): return (p / self._constant) ** (1/self._exp) def sound_velocity(self, rho, p): return self._constant * self._exp * rho ** (self._exp - 1) ``` ## Define Mixture EOS ```python from josie.twofluid.state import PhasePair class MixtureEOS(PhasePair): def __init__(self, eos1, eos2, surface_tension): super().__init__(eos1, eos2) self.surface_tension = surface_tension ``` ## Define Problem ```python import numpy as np from josie.problem import Problem from josie.dimension import MAX_DIMENSIONALITY from josie.math import Direction class SigmaProblem(Problem): # Small scale inertia sqrtm = np.sqrt(1) def __init__(self, eos: MixtureEOS): self.eos = eos def s(self, cells, t): values = cells.values.view(SigmaState) fields = values.fields nx, ny, _ = values.shape rhoYw = values[..., fields.rhoYw] p1 = values[..., fields.p1] p2 = values[..., fields.p2] # Store source ss = np.zeros_like(values) ss[..., fields.rhoa] = - rhoYw / self.sqrtm ss[..., fields.rhoYw] = (p2 - p1)/self.sqrtm return ss def F(self, cells): values = cells.values.view(SigmaState) nx, ny, _ = values.shape # Flux tensor F = np.zeros( ( nx, ny, len(SigmaConsFields), MAX_DIMENSIONALITY, ) ) fields = values.fields rho = values[..., fields.rho] rhoU = values[..., fields.rhoU] rhoV = values[..., fields.rhoV] rhoY = values[..., fields.rhoY] rhoYw = values[..., fields.rhoYw] rhoa = values[..., fields.rhoa] rhoz = values[..., fields.rhoz] U = values[..., fields.U] V = values[..., fields.V] P = values[..., fields.P] F[..., fields.rho, Direction.X] = rhoU F[..., fields.rho, Direction.Y] = rhoV F[..., fields.rhoU, Direction.X] = rhoU * U + P F[..., fields.rhoU, Direction.Y] = rhoU * V F[..., fields.rhoV, Direction.X] = rhoV * U F[..., fields.rhoV, Direction.Y] = rhoV * V + P F[..., fields.rhoY, Direction.X] = rhoY * U F[..., fields.rhoY, Direction.Y] = rhoY * V F[..., fields.rhoYw, Direction.X] = rhoYw * U F[..., fields.rhoYw, Direction.Y] = rhoYw * V F[..., fields.rhoa, Direction.X] = rhoa * U F[..., fields.rhoa, Direction.Y] = rhoa * V F[..., fields.rhoz, Direction.X] = rhoz * U F[..., fields.rhoz, Direction.Y] = rhoz * V return F ``` ## Define Scheme ```python import numpy as cp from josie.scheme.convective import ConvectiveScheme from josie.twofluid.fields import Phases class SigmaScheme(ConvectiveScheme): def __init__(self, eos: MixtureEOS): super().__init__(SigmaProblem(eos)) def post_step(self, cells): values = cells.values.view(SigmaState) fields = values.fields rho = values[..., fields.rho] rhoU = values[..., fields.rhoU] rhoV = values[..., fields.rhoV] rhoY = values[..., fields.rhoY] rhoYw = values[..., fields.rhoYw] rhoa = values[..., fields.rhoa] rhoz = values[..., fields.rhoz] U = rhoU / rho V = rhoV / rho Y = rhoY / rho Y2 = 1 - Y w = rhoYw / rho / Y alpha = rhoa / rho alpha2 = 1 - alpha z = rhoz / rho Sigma = (rho ** (1/3)) * (z ** (2/3)) surface_tension = self.problem.eos.surface_tension eos1 = self.problem.eos[Phases.PHASE1] rho1 = rho * Y / alpha p1 = eos1.p(rho1) c1 = eos1.sound_velocity(rho1, p1) rho2 = rho * Y2 / alpha2 eos2 = self.problem.eos[Phases.PHASE2] p2 = eos2.p(rho2) c2 = eos2.sound_velocity(rho1, p1) P = alpha * p1 + alpha2 * p2 + 0.5 * rhoYw**2 \ - 2/3 * surface_tension * Sigma c = cp.sqrt(Y * c1**2 + Y2 * c2**2 + rho * (Y * w)**2 \ + 2/9 * surface_tension * Sigma / rho) values[..., fields.U] = U values[..., fields.V] = V values[..., fields.P] = P values[..., fields.c] = c values[..., fields.Y] = Y values[..., fields.w] = w values[..., fields.alpha] = alpha values[..., fields.z] = z values[..., fields.p1] = p1 values[..., fields.p2] = p2 ``` ## Rusanov Scheme ```python from josie.euler.schemes import Rusanov as EulerRusanov from josie.general.schemes.time import ExplicitEuler from josie.general.schemes.source import ConstantSource class Rusanov(SigmaScheme, ExplicitEuler, ConstantSource): @staticmethod def compute_U_norm(values, normals): fields = values.fields # Get the velocity components UV_slice = slice(fields.U, fields.V + 1) UV = values[..., cp.newaxis, UV_slice] # Find the normal velocity U = cp.einsum("...kl,...l->...k", UV, normals) return U def CFL(self, cells, CFL_value): dt = super().CFL(cells, CFL_value) values = cells.values.view(SigmaState) fields = values.fields # Get the velocity components UV_slice = slice(fields.U, fields.V + 1) UV = cells.values[..., UV_slice] U = cp.linalg.norm(UV, axis=-1, keepdims=True) c = cells.values[..., fields.c] sigma = cp.max(cp.abs(U) + c[..., cp.newaxis]) # Min mesh dx dx = cells.min_length new_dt = CFL_value * dx / sigma return cp.min((dt, new_dt)) @staticmethod def compute_sigma(U_L, U_R, c_L, c_R): sigma_L = cp.abs(U_L) + c_L[..., cp.newaxis] sigma_R = cp.abs(U_R) + c_R[..., cp.newaxis] # Concatenate everything in a single array sigma_array = cp.concatenate((sigma_L, sigma_R), axis=-1) # And the we found the max on the last axis (i.e. the maximum value # of sigma for each cell) sigma = cp.max(sigma_array, axis=-1, keepdims=True) return sigma def F(self, cells, neighs): Q_L = cells.values.view(SigmaState) Q_R = neighs.values.view(SigmaState) fields = SigmaState.fields FS = cp.zeros_like(Q_L).view(SigmaState) # Get normal velocities U_L = self.compute_U_norm(Q_L, neighs.normals) U_R = self.compute_U_norm(Q_R, neighs.normals) # Speed of sound c_L = Q_L[..., fields.c] c_R = Q_R[..., fields.c] sigma = self.compute_sigma(U_L, U_R, c_L, c_R) DeltaF = 0.5 * (self.problem.F(cells) + self.problem.F(neighs)) # This is the flux tensor dot the normal DeltaF = cp.einsum("...kl,...l->...k", DeltaF, neighs.normals) # First four variables of the total state are the conservative # variables Q_L_cons = Q_L.get_conservative() Q_R_cons = Q_R.get_conservative() DeltaQ = 0.5 * sigma * (Q_R_cons - Q_L_cons) FS.set_conservative( neighs.surfaces[..., cp.newaxis] * (DeltaF - DeltaQ) ) return FS ``` ## Define Solver ```python from josie.solver import Solver class SigmaSolver(Solver): def __init__(self, mesh, scheme): super().__init__(mesh, SigmaState, scheme) ``` ## Define BCs ```python import numpy as np from josie.bc import Dirichlet, Neumann, make_periodic # Shortcut to define all Neumann conditions for everything except # vertical velocity Q0 = np.zeros(len(SigmaFields)).view(SigmaState) zero_normal_velocity = Neumann(Q0) # We override the bc for rhoV and V zero_normal_velocity.bc[..., [SigmaFields.rhoV, SigmaFields.V]] = \ Dirichlet(0) ``` ## Prepare the Domain ```python from josie.boundary import Line from josie.math import Direction left = Line([0, 0], [0, 1]) bottom = Line([0, 0], [1, 0]) right = Line([1, 0], [1, 1]) top = Line([0, 1], [1, 1]) # Left-right is periodic left, right = make_periodic(left, right, Direction.X) # Top-bottom top.bc = zero_normal_velocity bottom.bc = zero_normal_velocity ``` ## Generate the Mesh ```python from josie.mesh import Mesh from josie.mesh.cell import SimpleCell mesh = Mesh(left, bottom, right, top, SimpleCell) mesh.interpolate(150, 150) mesh.generate() ``` ## EOS ```python eos1 = Barotropic(a=1, b=1.4) eos2 = Barotropic(a=1, b=1.67) ``` ## Initial Conditions ```python # General P = 0.71 speed = 0.25 l = 0.05 # Reference length for the interface initial thickness We = 100 # Weber number # Top value alpha = 0.8 alpha2 = 1 - alpha w = 0 U = -speed V = 0 z = 0 p1 = p2 = P rho1 = eos1.rho(p1) c1 = eos1.sound_velocity(rho1, p1) rho2 = eos2.rho(p2) c2 = eos2.sound_velocity(rho2, p2) rho = alpha * rho1 + (1 - alpha) * rho2 Sigma = (rho ** (1/3)) * (z ** (2/3)) Y = rho1 * alpha / rho Y2 = 1 - Y rhoU = rho * U rhoV = rho * V rhoY = rho * Y rhoYw = rho * Y * w rhoa = rho * alpha rhoz = rho * z # EOS surface_tension = rho * speed**2 / We eos = MixtureEOS(eos1, eos2, surface_tension=0) c = np.sqrt(Y * c1**2 + Y2 * c2**2 + rho * (Y * w)**2 \ - 2/9 * eos.surface_tension * Sigma / rho) Q_top = SigmaState(rho, rhoU, rhoV, rhoY, rhoYw, rhoa, rhoz, U, V, P, c, Y, w, alpha, z, p1, p2) # Btm value alpha = 0.2 alpha2 = 1 - alpha w = 0 U = speed V = 0 z = 0 p1 = p2 = P rho1 = eos1.rho(p1) c1 = eos1.sound_velocity(rho1, p1) rho2 = eos2.rho(p2) c2 = eos2.sound_velocity(rho2, p2) rho = alpha * rho1 + (1 - alpha) * rho2 Sigma = (rho ** (1/3)) * (z ** (2/3)) Y = rho1 * alpha / rho Y2 = 1 - Y rhoU = rho * U rhoV = rho * V rhoY = rho * Y rhoYw = rho * Y * w rhoa = rho * alpha rhoz = rho * z c = np.sqrt(Y * c1**2 + Y2 * c2**2 + rho * (Y * w)**2 \ - 2/9 * eos.surface_tension * Sigma / rho) Q_btm = SigmaState(rho, rhoU, rhoV, rhoY, rhoYw, rhoa, rhoz, U, V, P, c, Y, w, alpha, z, p1, p2) # Interface value alpha = 0.5 alpha2 = 1 - alpha w = 0 U = 0 V = 0 z = 1e-2 p1 = p2 = P rho1 = eos1.rho(p1) c1 = eos1.sound_velocity(rho1, p1) rho2 = eos2.rho(p2) c2 = eos2.sound_velocity(rho2, p2) rho = alpha * rho1 + (1 - alpha) * rho2 Sigma = (rho ** (1/3)) * (z ** (2/3)) Y = rho1 * alpha / rho Y2 = 1 - Y rhoU = rho * U rhoV = rho * V rhoY = rho * Y rhoYw = rho * Y * w rhoa = rho * alpha rhoz = rho * z c = np.sqrt(Y * c1**2 + Y2 * c2**2 + rho * (Y * w)**2 \ - 2/9 * eos.surface_tension * Sigma / rho) Q_int = SigmaState(rho, rhoU, rhoV, rhoY, rhoYw, rhoa, rhoz, U, V, P, c, Y, w, alpha, z, p1, p2) def interface(x): """The curve defining the interface""" y = 0.5 * np.ones_like(x) # xa, xb are the interval where to put the perturbation xa = 0.65 xb = 0.85 K = 0.03 # Perturbation idx_pert = np.where((x >= xa) & (x <= xb)) y[idx_pert] += K * np.sin(np.pi * ((x[idx_pert] - xa)/(xb - xa))) return y def init_fun(cells): # First set everything to the interface value cells.values = Q_int dx = cells.min_length xc = cells.centroids[..., 0] yc = cells.centroids[..., 1] # Interface midline y_int = interface(xc) idx_top = np.where(yc > (y_int + dx)) idx_btm = np.where(yc < (y_int - dx)) # Change the value in the various part of the domain cells.values[idx_top[0], idx_top[1], ...] = Q_top cells.values[idx_btm[0], idx_btm[1], ...] = Q_btm ``` ## Initialize ```python from josie.io.write.writer import XDMFWriter from josie.io.write.strategy import TimeStrategy # Instantiate the solver solver = SigmaSolver(mesh, Rusanov(eos)) solver.init(init_fun) solver.plot() solver.show("rho") ``` ![png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAATt0lEQVR4nO2dX6hlZ3mHn3fmdCLSGKWpEOKfjHRSHO1F4tSam1ZxKJMUJhcOkgFjI9FBS0qpImSw0JJetSU2iGl1tGIV1KQtyECVgdqEgDghU9LE/GFkklozGppGbW7E6iRvL/YeZ/fknPWtnLX/fOt8zwObtfbZ76znl5Uz39nrt/fZE5mJiIhsf3asOoCIiCwHF3wRkUZwwRcRaQQXfBGRRnDBFxFpBBd8EZFGKC74EfG5iHgmIh7Z5PGIiE9ExJmIeDgirp5/TBERGUqfZ/ifBw50PH4tsGd6OwL87fBYIiIyb4oLfmbeB/yoY+R64As54STwyoi4bF4BRURkPqzN4RiXA0/N3D87/drT6wcj4giTqwCAt8zBLSLSGs9m5q9u5Q8u9UXbzDyWmfsycx/A/jjE+e3sfmm7zJlWnLXm0tleLp3Fmf9kq2Rm8QZcATyyyWOfBg7P3D8NXNbjmAlkELl+v7Rd5kwrzlpz6Wwvl87uGeBUn3V7o9s8nuEfB947fbfO24DnMvNFdY6IiKyYHs/Ev8ykj/85k37+ZuCDwAenjwdwJ/AE8G1gX8+rhtwfh/L8dna/tF3mTCvOWnPpbC+XzuLMlp/hx6o+HjkiEiAIkkmG8/ulbZ/Zec204qw1l872cunsngH+7fzroC8Vf9NWRKQRVrrg749DJPmLV6HP75e2y5xpxVlrLp3t5dLZPTMEKx2dVefS2V4unVY6IiIyECsdnVXn0tleLp3dM0Ow0tFZdS6d7eXSaaUjIiIDccEXEWkEO3ydVefS2V4und0zQ7DD11l1Lp3t5dJphy8iIgOx0tFZdS6d7eXS2T0zBCsdnVXn0tleLp1WOiIiMhArHZ1V59LZXi6d3TNDsNLRWXUune3l0mmlIyIiA7HS0Vl1Lp3t5dLZPTMEKx2dVefS2V4unVY6IiIyEBd8EZFGsMPXWXUune3l0tk9MwQ7fJ1V59LZXi6ddvgiIjIQKx2dVefS2V4und0zQ7DS0Vl1Lp3t5dJppSMiIgOx0tFZdS6d7eXS2T0zBCsdnVXn0tleLp1WOiIiMhArHZ1V59LZXi6d3TNDsNLRWXUune3l0mmlIyIiA3HBFxFpBDt8nVXn0tleLp3dM4PIzOINOACcBs4At27w+OuAe4AHgYeB63ocM4EMItfvl7bLnGnFWWsune3l0tk9A5zqs25vdCs+w4+IncCdwLXAXuBwROxdN/YnwN2ZeRVwA/A3peOKiMhy6VPpvBU4k5lPZubPgK8A16+bSeAV0/1LgB/0kW/ny66xOWvNpbO9XDq7ZwbRo3o5BHx25v6NwCfXzVwGfBs4C/wYeMsmxzoCnJretvVl19ictebS2V4und0zLLLS6clh4POZ+RrgOuCLEfGiY2fmsczct9X3kIqIyNbps+B/H3jtzP3XTL82y83A3QCZ+S3gZcClpQNv58uusTlrzaWzvVw6u2cG0aPSWQOeBHYDu4CHgDetm/k6cNN0/41MOvzwXTrjcdaaS2d7uXR2z7DISiczzwG3ACeAx5m8G+fRiLgtIg5Oxz4CfCAiHgK+zGTxz9KxRURkiWz1J8XQG5D741Ce387ul7bLnGnFWWsune3l0lmc2fIzfD88TWfVuXS2l0unH54mIiID8bN0dFadS2d7uXR2zwzBSkdn1bl0tpdLp5WOiIgMxAVfRKQR7PB1Vp1LZ3u5dHbPDMEOX2fVuXS2l0unHb6IiAzESkdn1bl0tpdLZ/fMEKx0dFadS2d7uXRa6YiIyECsdHRWnUtne7l0ds8MwUpHZ9W5dLaXS6eVjoiIDMRKR2fVuXS2l0tn98wQrHR0Vp1LZ3u5dFrpiIjIQFzwRUQawQ5fZ9W5dLaXS2f3zBDs8HVWnUtne7l02uGLiMhArHR0Vp1LZ3u5dHbPDMFKR2fVuXS2l0unlY6IiAzESkdn1bl0tpdLZ/fMEKx0dFadS2d7uXRa6YiIyECsdHRWnUtne7l0ds8MwUpHZ9W5dLaXS6eVjoiIDMRKR2fVuXS2l0tn98wQrHR0Vp1LZ3u5dFrpiIjIQHot+BFxICJOR8SZiLh1k5l3R8RjEfFoRHxpvjFFRGQwmdl5A3YCTwBvAHYBDwF7183sAR4EXjW9/+oex839cSjPb2f3S9tlzrTirDWXzvZy6SzOnCqtr5uuuz0W5muAEzP3jwJH1838JfD+lySGBDKIXL9f2i5zphVnrbl0tpdLZ/cMAxb8PpXO5cBTM/fPTr82y5XAlRHxzYg4GREHNjpQRByJiFMRcaqHV0RE5si8XrRdY1LrvB04DHwmIl65figzj2XmvvOvMO/fxm+dGpuz1lw628uls3tmEHOqdD4FvG/m/jeA37TSGY+z1lw628uls3uGBVc6DwB7ImJ3ROwCbgCOr5v5KpNn90TEpUwqnid7HFtERJZEccHPzHPALcAJ4HHg7sx8NCJui4iD07ETwA8j4jHgHuCjmfnD0rG382XX2Jy15tLZXi6d3TND8DdtdVadS2d7uXT6m7YiIjIQPzxNZ9W5dLaXS2f3zBCsdHRWnUtne7l0WumIiMhAXPBFRBrBDl9n1bl0tpdLZ/fMEOzwdVadS2d7uXTa4YuIyECsdHRWnUtne7l0ds8MwUpHZ9W5dLaXS6eVjoismNwRndt5zczOynyx0tFZdS6dy5/ZfcftL9ruvuN2eOGFzu1mjz3/9J7izGbH2yzPdjj/Wz3eEKx0dFadS+dyZ9ixA1544cVb2Pyxjpm1NTh3jpd2nB6z2+H8b/V4WOmIyFZ4KTXK2o4XOrebfa0003W8Uu6+2WWClY7OqnPpXFxdU6pgZrfPP72Hc+fo3M5rZnZ2szwbVUTr//tqPP/zON4gtvpPZQ298eJ/tmtb/TNkY3PWmkvn/GbYseP/bzf42toaG267Hpv3TJ/Zrv+G89vazv+8jseC/4nDhbLVV+9X8Y6B7eysNZfO+c1sxEupUWpio4poPbWd/3kdbxCrfIa/+47b8/x2dr+0XeZMK85ac+mc78zzT+/5xXZ2v7Rd5sw8j1fb+Z/T8bb8DH/llU6fS7M+l6ILm2nFWWsunXObWVS9sopKZ161zyj/n4+x0jl/cbLVV++XNdOKs9ZcOuc30xq1nf95HW8IK+/wRURkSayy0mmtT6zZWWsune3l0lmcGV+HH7TVJ9burDWXzvZy6eyeYYwdvoiILBkrHZ0159LZXi6dxRkrndouu8bmrDWXzvZy6eyewUpHRESKWOnorDmXzvZy6SzOWOnUdtk1NmetuXS2l0tn9wxWOiIiUsRKR2fNuXS2l0tnccZKp7bLrrE5a82ls71cOrtnsNIREZESLvgiIq1gh6+z5lw628ulsziz2A4fOACcBs4At3bMvQtIYJ8d/rictebS2V4und0zLLLDj4idwJ3AtcBe4HBE7N1g7mLgj4D7S8cUEZEV0OPZ/TXAiZn7R4GjG8zdAfwecC89nuHDtr/sGpWz1lw628ulszizuEoHOAR8dub+jcAn181cDfzTdP9eNlnwgSPAqektYftedo3NWWsune3l0tk9w4AFf42BRMQO4OPATaXZzDwGHAPYETH5ESAiIsthaKUDXAI8C3x3evsp8AMKtQ5s+8uuUTlrzaWzvVw6izMLrXTWgCeB3cAu4CHgTR3z9+K7dEbnrDWXzvZy6eyeYZHv0snMc8AtwAngceDuzHw0Im6LiIOlPy8iIpWw1Z8UQ2+w7S+7RuWsNZfO9nLpLM744Wm1XXaNzVlrLp3t5dLZPYMfniYiIkWsdHTWnEtne7l0FmesdGq77Bqbs9ZcOtvLpbN7BisdEREp4YIvItIKdvg6a86ls71cOoszdvi19Wxjc9aaS2d7uXR2z2CHLyIiRax0dNacS2d7uXQWZ6x0arvsGpuz1lw628uls3sGKx0RESlipaOz5lw628ulszhjpVPbZdfYnLXm0tleLp3dM1jpiIhIESsdnTXn0tleLp3FGSud2i67xuasNZfO9nLp7J7BSkdEREq44IuItIIdvs6ac+lsL5fO4owdfm0929ictebS2V4und0z2OGLiEgRKx2dNefS2V4uncUZK53aLrvG5qw1l872cunsnsFKR0REiljp6Kw5l872cukszljp1HbZNTZnrbl0tpdLZ/cMVjoiIlLESkdnzbl0tpdLZ3HGSqe2y66xOWvNpbO9XDq7Z7DSERGRIlY6OmvOpbO9XDqLM1Y6tV12jc1Zay6d7eXS2T2DlY6IiJToteBHxIGIOB0RZyLi1g0e/3BEPBYRD0fENyLi9fOPKiIig+jRte8EngDeAOwCHgL2rpt5B/Dy6f6HgLvs8MflrDWXzvZy6SzOLK7DB64BTszcPwoc7Zi/CvimHf64nLXm0tleLp3dMyy4w78ceGrm/tnp1zbjZuDrGz0QEUci4lREnMoeYhERmSM9nuEfAj47c/9G4JObzL4HOAlcZKUzLmetuXS2l0tncWb1lQ6wH3gceHUfsZVOXc5ac+lsL5fO7hkWXOk8AOyJiN0RsQu4ATg+OxARVwGfBg5m5jM9jikiIsum12UAXAd8h8m7dT42/dptTBZ4gH8B/gv49+ntuJXOuJy15tLZXi6dxRl/07a2y66xOWvNpbO9XDq7Z/A3bUVEpMiqnuHDtr/sGpWz1lw628ulszhjpVPbZdfYnLXm0tleLp3dM1jpiIhICRd8EZFWsMPXWXMune3l0lmcscOvrWcbm7PWXDrby6WzewY7fBERKWKlo7PmXDrby6WzOGOlU9tl19ictebS2V4und0zWOmIiEgRKx2dNefS2V4uncUZK53aLrvG5qw1l872cunsnsFKR0REiljp6Kw5l872cukszljp1HbZNTZnrbl0tpdLZ/cMVjoiIlLESkdnzbl0tpdLZ3HGSqe2y66xOWvNpbO9XDq7Z7DSERGREi74IiKtYIevs+ZcOtvLpbM4Y4dfW882NmetuXS2l0tn9wx2+CIiUsRKR2fNuXS2l0tnccZKp7bLrrE5a82ls71cOrtnsNIREZEiVjo6a86ls71cOoszVjq1XXaNzVlrLp3t5dLZPYOVjoiIFLHS0VlzLp3t5dJZnLHSqe2ya2zOWnPpbC+Xzu4ZrHRERKSEC76ISCvY4eusOZfO9nLpLM4stsMHDgCngTPArRs8fhFw1/Tx+4Er7PDH5aw1l872cunsnmGRHX5E7ATuBK4F9gKHI2LvurGbgR9n5q8Bfw38Rem4IiKyZHo8u78GODFz/yhwdN3MCeCa6f4a8CwQVjrjcdaaS2d7uXQWZ7b8DD+mi++mRMQh4EBmvn96/0bgtzLzlpmZR6YzZ6f3n5jOPLvuWEeAI9O7bwYe6ZS3w6VMfkiK52IWz8UFPBcX+PXMvHgrf3Bt3km6yMxjwDGAiDiVmfuW6a8Vz8UFPBcX8FxcwHNxgYg4tdU/2+dtmd8HXjtz/zXTr204ExFrwCXAD7caSkRE5k+fBf8BYE9E7I6IXcANwPF1M8eB35/uHwL+NUtdkYiILJVipZOZ5yLiFiYvzO4EPpeZj0bEbUxePDgO/B3wxYg4A/yIyQ+FEscG5N5ueC4u4Lm4gOfiAp6LC2z5XBRftBURke2BH60gItIILvgiIo2w8AU/Ig5ExOmIOBMRt27w+EURcdf08fsj4opFZ1oVPc7FhyPisYh4OCK+ERGvX0XOZVA6FzNz74qIjIht+5a8PuciIt49/d54NCK+tOyMy6LH35HXRcQ9EfHg9O/JdavIuWgi4nMR8cz0d5w2ejwi4hPT8/RwRFzd68AL/oC0ncATwBuAXcBDwN51M38AfGq6fwNw1yIzrerW81y8A3j5dP9DLZ+L6dzFwH3ASWDfqnOv8PtiD/Ag8Krp/VevOvcKz8Ux4EPT/b3Ad1ede0Hn4reBq4FHNnn8OuDrQABvA+7vc9xFP8N/K3AmM5/MzJ8BXwGuXzdzPfD30/1/BN4ZEbHgXKugeC4y857M/Mn07kkmv/OwHenzfQHw50w+l+mnywy3ZPqciw8Ad2bmjwEy85klZ1wWfc5FAq+Y7l8C/GCJ+ZZGZt7H5B2Pm3E98IWccBJ4ZURcVjruohf8y4GnZu6fnX5tw5nMPAc8B/zKgnOtgj7nYpabmfwE344Uz8X0EvW1mfnPywy2Avp8X1wJXBkR34yIkxFxYGnplkufc/FnwHsi4izwNeAPlxOtOl7qegIs+aMVpB8R8R5gH/A7q86yCiJiB/Bx4KYVR6mFNSa1ztuZXPXdFxG/kZn/s8pQK+Iw8PnMvD0irmHy+z9vzswXVh1sDCz6Gb4fy3CBPueCiNgPfAw4mJn/u6Rsy6Z0Li5m8uF690bEd5l0lMe36Qu3fb4vzgLHM/PnmfkfwHeY/ADYbvQ5FzcDdwNk5reAlzH5YLXW6LWerGfRC74fy3CB4rmIiKuATzNZ7LdrTwuFc5GZz2XmpZl5RWZeweT1jIOZueUPjaqYPn9Hvsrk2T0RcSmTiufJJWZcFn3OxfeAdwJExBuZLPj/vdSUdXAceO/03TpvA57LzKdLf2ihlU4u7mMZRkfPc/FXwC8D/zB93fp7mXlwZaEXRM9z0QQ9z8UJ4Hcj4jHgeeCjmbntroJ7nouPAJ+JiD9m8gLuTdvxCWJEfJnJD/lLp69X/CnwSwCZ+Skmr19cx+RfGfwJ8L5ex92G50pERDbA37QVEWkEF3wRkUZwwRcRaQQXfBGRRnDBFxFpBBd8EZFGcMEXEWmE/wMzsRPeojTk3wAAAABJRU5ErkJggg== ) ## Logging ```python import logging from datetime import datetime now = datetime.now().strftime('%Y%m%d%H%M%S') logger = logging.getLogger("josie") logger.setLevel(logging.DEBUG) fh = logging.FileHandler(f"sigma-{now}.log") fh.setLevel(logging.DEBUG) formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') fh.setFormatter(formatter) logger.addHandler(fh) ``` ```python # Write strategy strategy = TimeStrategy(dt_save=0.1, animate=False) writer = XDMFWriter(f"sigma-{now}.xdmf", strategy, solver, final_time=5, CFL=0.5) writer.solve() ```