examples/004_Two-phase-sigma.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1b18d0fa",
   "metadata": {},
   "source": [
    "## Define fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e5256c92",
   "metadata": {},
   "outputs": [],
   "source": [
    "from josie.fluid.fields import FluidFields\n",
    "from josie.fields import Fields\n",
    "\n",
    "class SigmaFields(FluidFields):\n",
    "    # Conservative\n",
    "    rho = 0\n",
    "    rhoU = 1\n",
    "    rhoV = 2\n",
    "    rhoY = 3\n",
    "    rhoYw = 4\n",
    "    rhoa = 5\n",
    "    rhoz = 6\n",
    "    \n",
    "    # Auxiliary\n",
    "    U = 7\n",
    "    V = 8\n",
    "    P = 9\n",
    "    c = 10\n",
    "    Y = 11\n",
    "    w = 12\n",
    "    alpha = 13\n",
    "    z = 14\n",
    "    p1 = 15\n",
    "    p2 = 16\n",
    "\n",
    "    \n",
    "class SigmaConsFields(Fields):\n",
    "    rho = 0\n",
    "    rhoU = 1\n",
    "    rhoV = 2\n",
    "    rhoY = 3\n",
    "    rhoYw = 4\n",
    "    rhoa = 5\n",
    "    rhoz = 6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a0c55d0",
   "metadata": {},
   "source": [
    "## Define State"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5572ea8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from josie.fluid.state import SingleFluidState\n",
    "from josie.state import SubsetState\n",
    "\n",
    "class SigmaConsState(SubsetState):\n",
    "    full_state_fields = SigmaFields\n",
    "    fields = SigmaConsFields\n",
    "    \n",
    "class SigmaState(SingleFluidState):\n",
    "    fields = SigmaFields\n",
    "    cons_state = SigmaConsState\n",
    "    prim_state = SigmaConsState\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4582417",
   "metadata": {},
   "source": [
    "## Define Barotropic EOS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1623c4a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Barotropic:\n",
    "    \"\"\"\n",
    "    .. math::\n",
    "    \n",
    "        \\pressure = a\\rho^b\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, a, b):\n",
    "        self._constant = a\n",
    "        self._exp = b\n",
    "    \n",
    "    def p(self, rho):\n",
    "        return self._constant * rho ** self._exp\n",
    "    \n",
    "    def rho(self, p):\n",
    "        return (p / self._constant) ** (1/self._exp)\n",
    "    \n",
    "    def sound_velocity(self, rho, p):\n",
    "        return self._constant * self._exp * rho ** (self._exp - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dff0be1",
   "metadata": {},
   "source": [
    "## Define Mixture EOS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ff95629b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from josie.twofluid.state import PhasePair\n",
    "\n",
    "class MixtureEOS(PhasePair):\n",
    "    def __init__(self, eos1, eos2, surface_tension):\n",
    "        super().__init__(eos1, eos2)\n",
    "        self.surface_tension = surface_tension"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93527dba",
   "metadata": {},
   "source": [
    "## Define Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3c9267c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from josie.problem import Problem\n",
    "from josie.dimension import MAX_DIMENSIONALITY\n",
    "from josie.math import Direction\n",
    "\n",
    "class SigmaProblem(Problem):\n",
    "    \n",
    "    # Small scale inertia\n",
    "    sqrtm = np.sqrt(1)\n",
    "    \n",
    "    def __init__(self, eos: MixtureEOS):\n",
    "        self.eos = eos\n",
    "        \n",
    "    def s(self, cells, t):\n",
    "        values = cells.values.view(SigmaState)\n",
    "        fields = values.fields\n",
    "        nx, ny, num_dofs, _ = values.shape\n",
    "        \n",
    "        rhoYw = values[..., [0], fields.rhoYw]\n",
    "        p1 = values[..., [0], fields.p1]\n",
    "        p2 = values[..., [0], fields.p2]\n",
    "        \n",
    "        \n",
    "        # Store source\n",
    "        ss = np.zeros_like(values[..., [0], :])\n",
    "        \n",
    "        ss[..., fields.rhoa] = - rhoYw / self.sqrtm\n",
    "        ss[..., fields.rhoYw] = (p2 - p1)/self.sqrtm\n",
    "        \n",
    "        return ss\n",
    "        \n",
    "        \n",
    "    def F(self, values: SigmaState):\n",
    "        nx, ny, num_dofs, _ = values.shape\n",
    "        \n",
    "        # Flux tensor\n",
    "        F = np.zeros(\n",
    "            (\n",
    "                nx,\n",
    "                ny,\n",
    "                num_dofs,\n",
    "                len(SigmaConsFields),\n",
    "                MAX_DIMENSIONALITY,\n",
    "            )\n",
    "        )\n",
    "        \n",
    "        fields = values.view(SigmaState).fields\n",
    "                  \n",
    "        rho = values[..., fields.rho]\n",
    "        rhoU = values[..., fields.rhoU]\n",
    "        rhoV = values[..., fields.rhoV]\n",
    "        rhoY = values[..., fields.rhoY]\n",
    "        rhoYw = values[..., fields.rhoYw]\n",
    "        rhoa = values[..., fields.rhoa]\n",
    "        rhoz = values[..., fields.rhoz]\n",
    "        \n",
    "        U = values[..., fields.U]\n",
    "        V = values[..., fields.V]\n",
    "        P = values[..., fields.P]\n",
    "        \n",
    "       \n",
    "        F[..., fields.rho, Direction.X] = rhoU\n",
    "        F[..., fields.rho, Direction.Y] = rhoV\n",
    "        F[..., fields.rhoU, Direction.X] = rhoU * U + P\n",
    "        F[..., fields.rhoU, Direction.Y] = rhoU * V\n",
    "        F[..., fields.rhoV, Direction.X] = rhoV * U\n",
    "        F[..., fields.rhoV, Direction.Y] = rhoV * V + P\n",
    "        F[..., fields.rhoY, Direction.X] = rhoY * U\n",
    "        F[..., fields.rhoY, Direction.Y] = rhoY * V\n",
    "        F[..., fields.rhoYw, Direction.X] = rhoYw * U\n",
    "        F[..., fields.rhoYw, Direction.Y] = rhoYw * V\n",
    "        F[..., fields.rhoa, Direction.X] = rhoa * U\n",
    "        F[..., fields.rhoa, Direction.Y] = rhoa * V\n",
    "        F[..., fields.rhoz, Direction.X] = rhoz * U\n",
    "        F[..., fields.rhoz, Direction.Y] = rhoz * V\n",
    "        \n",
    "        return F\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af378d54",
   "metadata": {},
   "source": [
    "## Define Scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b3dba309",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as cp\n",
    "\n",
    "from josie.scheme.convective import ConvectiveScheme\n",
    "from josie.twofluid.fields import Phases\n",
    "\n",
    "class SigmaScheme(ConvectiveScheme):\n",
    "    def __init__(self, eos: MixtureEOS):\n",
    "        super().__init__(SigmaProblem(eos))\n",
    "        \n",
    "    def post_step(self, values: SigmaState):\n",
    "        fields = values.view(SigmaState).fields\n",
    "        \n",
    "        rho = values[..., fields.rho]\n",
    "        rhoU = values[..., fields.rhoU]\n",
    "        rhoV = values[..., fields.rhoV]\n",
    "        rhoY = values[..., fields.rhoY]\n",
    "        rhoYw = values[..., fields.rhoYw]\n",
    "        rhoa = values[..., fields.rhoa]\n",
    "        rhoz = values[..., fields.rhoz]\n",
    "        \n",
    "        U = rhoU / rho\n",
    "        V = rhoV / rho\n",
    "        Y = rhoY / rho\n",
    "        Y2 = 1 - Y\n",
    "        w = rhoYw / rho / Y\n",
    "        alpha = rhoa / rho\n",
    "        alpha2 = 1 - alpha\n",
    "        z = rhoz / rho\n",
    "        Sigma = (rho ** (1/3)) * (z ** (2/3))\n",
    "        \n",
    "        surface_tension = self.problem.eos.surface_tension\n",
    "        \n",
    "        eos1 = self.problem.eos[Phases.PHASE1]\n",
    "        rho1 = rho * Y / alpha\n",
    "        p1 = eos1.p(rho1)\n",
    "        c1 = eos1.sound_velocity(rho1, p1)\n",
    "        \n",
    "        \n",
    "        rho2 = rho * Y2 / alpha2 \n",
    "        eos2 = self.problem.eos[Phases.PHASE2]\n",
    "        p2 = eos2.p(rho2)\n",
    "        c2 = eos2.sound_velocity(rho1, p1)\n",
    "        \n",
    "        P = alpha * p1 + alpha2 * p2 + 0.5 * rhoYw**2 \\\n",
    "            - 2/3 * surface_tension * Sigma\n",
    "        \n",
    "        c = cp.sqrt(Y * c1**2 + Y2 * c2**2  + rho * (Y * w)**2 \\\n",
    "            + 2/9 * surface_tension * Sigma / rho)\n",
    "               \n",
    "        values[..., fields.U] = U\n",
    "        values[..., fields.V] = V\n",
    "        values[..., fields.P] = P\n",
    "        values[..., fields.c] = c\n",
    "        values[..., fields.Y] = Y\n",
    "        values[..., fields.w] = w\n",
    "        values[..., fields.alpha] = alpha\n",
    "        values[..., fields.z] = z\n",
    "        values[..., fields.p1] = p1\n",
    "        values[..., fields.p2] = p2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea2e84c9",
   "metadata": {},
   "source": [
    "## Rusanov Scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0231a559",
   "metadata": {},
   "outputs": [],
   "source": [
    "from josie.euler.schemes import Rusanov as EulerRusanov\n",
    "from josie.general.schemes.time.rk import RK2\n",
    "from josie.general.schemes.source import ConstantSource\n",
    "from josie.general.schemes.space.muscl import MUSCL\n",
    "from josie.general.schemes.space.limiters import MinMod\n",
    "\n",
    "\n",
    "\n",
    "class Rusanov(ConstantSource, RK2, SigmaScheme):\n",
    "    @staticmethod\n",
    "    def compute_U_norm(values: SigmaState, normals):\n",
    "        fields = values.view(SigmaState).fields\n",
    "\n",
    "        # Get the velocity components\n",
    "        UV_slice = slice(fields.U, fields.V + 1)\n",
    "        UV = values[..., cp.newaxis, UV_slice]\n",
    "\n",
    "        # Find the normal velocity\n",
    "        U = cp.einsum(\"...mkl,...l->...mk\", UV, normals)\n",
    "\n",
    "        return U\n",
    "\n",
    "    def CFL(self, cells, CFL_value):\n",
    "\n",
    "        #dt = super().CFL(cells, CFL_value)\n",
    "\n",
    "        values = cells.values.view(SigmaState)\n",
    "        fields = values.fields\n",
    "\n",
    "        # Get the velocity components\n",
    "        UV_slice = slice(fields.U, fields.V + 1)\n",
    "        UV = cells.values[..., UV_slice]\n",
    "\n",
    "        U = cp.linalg.norm(UV, axis=-1, keepdims=True)\n",
    "        c = cells.values[..., fields.c]\n",
    "\n",
    "        sigma = cp.max(cp.abs(U) + c[..., cp.newaxis])\n",
    "\n",
    "        # Min mesh dx\n",
    "        dx = cells.min_length\n",
    "\n",
    "        new_dt = CFL_value * dx / sigma\n",
    "\n",
    "        #return cp.min((dt, new_dt))\n",
    "        return new_dt\n",
    "    \n",
    "    @staticmethod\n",
    "    def compute_sigma(U_L, U_R, c_L, c_R):\n",
    "\n",
    "        sigma_L = cp.abs(U_L) + c_L[..., cp.newaxis]\n",
    "\n",
    "        sigma_R = cp.abs(U_R) + c_R[..., cp.newaxis]\n",
    "\n",
    "        # Concatenate everything in a single array\n",
    "        sigma_array = cp.concatenate((sigma_L, sigma_R), axis=-1)\n",
    "\n",
    "        # And the we found the max on the last axis (i.e. the maximum value\n",
    "        # of sigma for each cell)\n",
    "        sigma = cp.max(sigma_array, axis=-1, keepdims=True)\n",
    "\n",
    "        return sigma\n",
    "\n",
    "    def intercellFlux(\n",
    "        self, Q_L: SigmaState, Q_R: SigmaState, normals: np.ndarray, surfaces: np.ndarray\n",
    "    ) -> SigmaState:\n",
    "\n",
    "        fields = SigmaState.fields\n",
    "\n",
    "        FS = cp.zeros_like(Q_L).view(SigmaState)\n",
    "        \n",
    "        # Get normal velocities\n",
    "        U_L = self.compute_U_norm(Q_L, normals)\n",
    "        U_R = self.compute_U_norm(Q_R, normals)\n",
    "\n",
    "        # Speed of sound\n",
    "        c_L = Q_L[..., fields.c]\n",
    "        c_R = Q_R[..., fields.c]\n",
    "\n",
    "        sigma = self.compute_sigma(U_L, U_R, c_L, c_R)\n",
    "        \n",
    "        DeltaF = 0.5 * (self.problem.F(Q_L) + self.problem.F(Q_R))\n",
    "\n",
    "        # This is the flux tensor dot the normal\n",
    "        DeltaF = cp.einsum(\"...mkl,...l->...mk\", DeltaF, normals)\n",
    "\n",
    "        # First four variables of the total state are the conservative\n",
    "        # variables\n",
    "        Q_L_cons = Q_L.view(SigmaState).get_conservative()\n",
    "        Q_R_cons = Q_R.view(SigmaState).get_conservative()\n",
    "\n",
    "        DeltaQ = 0.5 * sigma * (Q_R_cons - Q_L_cons)\n",
    "\n",
    "        FS.set_conservative(\n",
    "            surfaces[..., cp.newaxis, cp.newaxis] * (DeltaF - DeltaQ)\n",
    "        )\n",
    "\n",
    "        return FS\n",
    "    \n",
    "class FinalScheme(\n",
    "        Rusanov, MUSCL, MinMod,\n",
    "    ):\n",
    "        pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cb36013",
   "metadata": {},
   "source": [
    "## Define Solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5d3bafcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "from josie.solver import Solver\n",
    "\n",
    "class SigmaSolver(Solver):\n",
    "    def __init__(self, mesh, scheme):\n",
    "\n",
    "        super().__init__(mesh, SigmaState, scheme)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09141c12",
   "metadata": {},
   "source": [
    "## Define BCs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ea5a3cd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from josie.bc import Dirichlet, Neumann, make_periodic\n",
    "\n",
    "    \n",
    "# Shortcut to define all Neumann conditions for everything except \n",
    "# vertical velocity\n",
    "Q0 = np.zeros(len(SigmaFields)).view(SigmaState)\n",
    "zero_normal_velocity = Neumann(Q0)\n",
    "\n",
    "# We override the bc for rhoV and V\n",
    "zero_normal_velocity.bc[..., [SigmaFields.rhoV, SigmaFields.V]] = \\\n",
    "    Dirichlet(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb1fd8ca",
   "metadata": {},
   "source": [
    "## Prepare the Domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "b5b88916",
   "metadata": {},
   "outputs": [],
   "source": [
    "from josie.boundary import Line\n",
    "from josie.math import Direction\n",
    "\n",
    "left = Line([0, 0], [0, 1])\n",
    "bottom = Line([0, 0], [1, 0])\n",
    "right = Line([1, 0], [1, 1])\n",
    "top = Line([0, 1], [1, 1])\n",
    "\n",
    "# Left-right is periodic\n",
    "left, right = make_periodic(left, right, Direction.X)\n",
    "\n",
    "# Top-bottom\n",
    "top.bc = zero_normal_velocity\n",
    "bottom.bc = zero_normal_velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c1e649a",
   "metadata": {},
   "source": [
    "## Generate the Mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "78c67ea5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from josie.mesh import Mesh\n",
    "from josie.mesh.cell import MUSCLCell\n",
    "\n",
    "mesh = Mesh(left, bottom, right, top, MUSCLCell)\n",
    "mesh.interpolate(150, 150)\n",
    "mesh.generate()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "366252c0",
   "metadata": {},
   "source": [
    "## EOS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "7ade9cfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "eos1 = Barotropic(a=1, b=1.4)\n",
    "eos2 = Barotropic(a=1, b=1.67)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "462fa019",
   "metadata": {},
   "source": [
    "## Initial Conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8f60d060",
   "metadata": {},
   "outputs": [],
   "source": [
    "# General\n",
    "P = 0.71\n",
    "speed = 0.25\n",
    "\n",
    "l = 0.05 # Reference length for the interface initial thickness\n",
    "We = 100 # Weber number\n",
    "\n",
    "\n",
    "\n",
    "# Top value\n",
    "alpha = 0.8\n",
    "alpha2 = 1 - alpha\n",
    "w = 0\n",
    "U = -speed\n",
    "V = 0\n",
    "z = 0\n",
    "p1 = p2 = P\n",
    "rho1 = eos1.rho(p1)\n",
    "c1 = eos1.sound_velocity(rho1, p1)\n",
    "\n",
    "rho2 = eos2.rho(p2)\n",
    "c2 = eos2.sound_velocity(rho2, p2)\n",
    "\n",
    "rho = alpha * rho1 + (1 - alpha) * rho2\n",
    "Sigma = (rho ** (1/3)) * (z ** (2/3))\n",
    "Y = rho1 * alpha / rho\n",
    "Y2 = 1 - Y\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "rhoU = rho * U\n",
    "rhoV = rho * V\n",
    "rhoY = rho * Y\n",
    "rhoYw = rho * Y * w\n",
    "rhoa = rho * alpha\n",
    "rhoz = rho * z\n",
    "\n",
    "# EOS\n",
    "surface_tension = rho * speed**2 / We\n",
    "\n",
    "eos = MixtureEOS(eos1, eos2, surface_tension=0)\n",
    "\n",
    "c = np.sqrt(Y * c1**2 + Y2 * c2**2  + rho * (Y * w)**2 \\\n",
    "            - 2/9 * eos.surface_tension * Sigma / rho)\n",
    "\n",
    "Q_top = SigmaState(rho, rhoU, rhoV, rhoY, rhoYw, rhoa, rhoz,\n",
    "                   U, V, P, c, Y, w, alpha, z, p1, p2)\n",
    "\n",
    "# Btm value\n",
    "alpha = 0.2\n",
    "alpha2 = 1 - alpha\n",
    "w = 0\n",
    "U = speed\n",
    "V = 0\n",
    "z = 0\n",
    "p1 = p2 = P\n",
    "rho1 = eos1.rho(p1)\n",
    "c1 = eos1.sound_velocity(rho1, p1)\n",
    "\n",
    "rho2 = eos2.rho(p2)\n",
    "c2 = eos2.sound_velocity(rho2, p2)\n",
    "\n",
    "rho = alpha * rho1 + (1 - alpha) * rho2\n",
    "Sigma = (rho ** (1/3)) * (z ** (2/3))\n",
    "Y = rho1 * alpha / rho\n",
    "Y2 = 1 - Y\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "rhoU = rho * U\n",
    "rhoV = rho * V\n",
    "rhoY = rho * Y\n",
    "rhoYw = rho * Y * w\n",
    "rhoa = rho * alpha\n",
    "rhoz = rho * z\n",
    "\n",
    "c = np.sqrt(Y * c1**2 + Y2 * c2**2  + rho * (Y * w)**2 \\\n",
    "            - 2/9 * eos.surface_tension * Sigma / rho)\n",
    "\n",
    "\n",
    "Q_btm = SigmaState(rho, rhoU, rhoV, rhoY, rhoYw, rhoa, rhoz,\n",
    "                   U, V, P, c, Y, w, alpha, z, p1, p2)\n",
    "\n",
    "# Interface value\n",
    "alpha = 0.5\n",
    "alpha2 = 1 - alpha\n",
    "w = 0\n",
    "U = 0\n",
    "V = 0\n",
    "z = 1e-2\n",
    "\n",
    "p1 = p2 = P\n",
    "rho1 = eos1.rho(p1)\n",
    "c1 = eos1.sound_velocity(rho1, p1)\n",
    "\n",
    "rho2 = eos2.rho(p2)\n",
    "c2 = eos2.sound_velocity(rho2, p2)\n",
    "\n",
    "rho = alpha * rho1 + (1 - alpha) * rho2\n",
    "Sigma = (rho ** (1/3)) * (z ** (2/3))\n",
    "\n",
    "Y = rho1 * alpha / rho\n",
    "Y2 = 1 - Y\n",
    "\n",
    "\n",
    "rhoU = rho * U\n",
    "rhoV = rho * V\n",
    "rhoY = rho * Y\n",
    "rhoYw = rho * Y * w\n",
    "rhoa = rho * alpha\n",
    "rhoz = rho * z\n",
    "\n",
    "c = np.sqrt(Y * c1**2 + Y2 * c2**2  + rho * (Y * w)**2 \\\n",
    "            - 2/9 * eos.surface_tension * Sigma / rho)\n",
    "\n",
    "\n",
    "Q_int = SigmaState(rho, rhoU, rhoV, rhoY, rhoYw, rhoa, rhoz,\n",
    "                   U, V, P, c, Y, w, alpha, z, p1, p2)\n",
    "\n",
    "def interface(x):\n",
    "    \"\"\"The curve defining the interface\"\"\"\n",
    "    \n",
    "    y = 0.5 * np.ones_like(x)\n",
    "    \n",
    "    # xa, xb are the interval where to put the perturbation\n",
    "    xa = 0.65\n",
    "    xb = 0.85\n",
    "    K = 0.03\n",
    "    \n",
    "    # Perturbation\n",
    "  \n",
    "    idx_pert = np.where((x >= xa) & (x <= xb))\n",
    "    y[idx_pert] += K * np.sin(np.pi * ((x[idx_pert] - xa)/(xb - xa)))\n",
    "    \n",
    "    return y\n",
    "\n",
    "\n",
    "def init_fun(cells):\n",
    "    # First set everything to the interface value\n",
    "    cells.values = Q_int\n",
    "    \n",
    "    dx = cells.min_length\n",
    "    xc = cells.centroids[..., 0]\n",
    "    yc = cells.centroids[..., 1]\n",
    "    \n",
    "    # Interface midline\n",
    "    y_int = interface(xc)\n",
    "    \n",
    "    idx_top = np.where(yc > (y_int + dx))\n",
    "    idx_btm = np.where(yc < (y_int - dx))\n",
    "\n",
    "    \n",
    "    # Change the value in the various part of the domain\n",
    "    cells.values[idx_top[0], idx_top[1], ...] = Q_top\n",
    "    cells.values[idx_btm[0], idx_btm[1], ...] = Q_btm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25f0dc44",
   "metadata": {},
   "source": [
    "## Initialize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0896682d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAATt0lEQVR4nO2dX6hlZ3mHn3fmdCLSGKWpEOKfjHRSHO1F4tSam1ZxKJMUJhcOkgFjI9FBS0qpImSw0JJetSU2iGl1tGIV1KQtyECVgdqEgDghU9LE/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==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from josie.io.write.writer import XDMFWriter\n",
    "from josie.io.write.strategy import TimeStrategy\n",
    "\n",
    "# Instantiate the solver\n",
    "solver = SigmaSolver(mesh, FinalScheme(eos))\n",
    "solver.init(init_fun)\n",
    "\n",
    "solver.plot()\n",
    "solver.show(\"rho\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "clear-aluminum",
   "metadata": {},
   "source": [
    "## Logging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "provincial-masters",
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "from datetime import datetime\n",
    "\n",
    "now = datetime.now().strftime('%Y%m%d%H%M%S')\n",
    "\n",
    "logger = logging.getLogger(\"josie\")\n",
    "logger.setLevel(logging.DEBUG)\n",
    "\n",
    "fh = logging.FileHandler(f\"sigma-{now}.log\")\n",
    "fh.setLevel(logging.DEBUG)\n",
    "formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')\n",
    "fh.setFormatter(formatter)\n",
    "\n",
    "logger.addHandler(fh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "regular-digest",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/jq/xw5g5cn92_x68gwfg8krm8m80000gn/T/ipykernel_21380/1453061573.py:21: RuntimeWarning: invalid value encountered in true_divide\n",
      "  U = rhoU / rho\n",
      "/var/folders/jq/xw5g5cn92_x68gwfg8krm8m80000gn/T/ipykernel_21380/1453061573.py:22: RuntimeWarning: invalid value encountered in true_divide\n",
      "  V = rhoV / rho\n",
      "/var/folders/jq/xw5g5cn92_x68gwfg8krm8m80000gn/T/ipykernel_21380/1453061573.py:23: RuntimeWarning: invalid value encountered in true_divide\n",
      "  Y = rhoY / rho\n",
      "/var/folders/jq/xw5g5cn92_x68gwfg8krm8m80000gn/T/ipykernel_21380/1453061573.py:25: RuntimeWarning: invalid value encountered in true_divide\n",
      "  w = rhoYw / rho / Y\n",
      "/var/folders/jq/xw5g5cn92_x68gwfg8krm8m80000gn/T/ipykernel_21380/1453061573.py:26: RuntimeWarning: invalid value encountered in true_divide\n",
      "  alpha = rhoa / rho\n",
      "/var/folders/jq/xw5g5cn92_x68gwfg8krm8m80000gn/T/ipykernel_21380/1453061573.py:28: RuntimeWarning: invalid value encountered in true_divide\n",
      "  z = rhoz / rho\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# Write strategy\n",
    "strategy = TimeStrategy(dt_save=0.1, animate=False)\n",
    "writer = XDMFWriter(f\"sigma-{now}.xdmf\", strategy, solver, \n",
    "                    final_time=5, CFL=0.5)\n",
    "\n",
    "writer.solve()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}