Implementing a solver for the advection equation¶

We’re going to write all the necessary objects to solve the following equation:

$\partial_t \varphi + \mathbf{u} \cdot \nabla \varphi = 0$

That, compared to our generalized model:

$\partial_t \mathbf{q} + \nabla \cdot \left(\underline{\underline{\mathbf{F}}}(\mathbf{q}) + \underline{\underline{\underline{\underline{D(\mathbf{q})}}}} \cdot \nabla \mathbf{q} \right) + \underline{\underline{\underline{\mathbf{B}(\mathbf{q})}}} \cdot \nabla\mathbf{q} = \mathbf{s}(\mathbf{q})$

Gives $\mathbf{q} = \left\{\varphi\right\}$ $\mathbf{F} = \underline{\underline{\mathbf{0}}}$ , $\mathbf{D} = \underline{\underline{\underline{\underline{\mathbf{0}}}}}$

So our problem needs just to provide $\mathbf{B}$ .

$B_{pqr} = u_r$

`State`¶

from josie.state import State
from josie.fields import Fields

class AdvectionFields(Fields):
    phi = 0
    
class AdvectionState(State):
    fields = AdvectionFields

`Problem`¶

We just need to implement the $\mathbf{B}(\mathbf{q})$ operator

from josie.problem import Problem

Problem??

from josie.mesh.cellset import CellSet, MeshCellSet

class AdvectionProblem(Problem):
    V = np.array([1, 0])
    
    DIMENSIONALITY = 2 ## 2D
    
    def B(self, cells: CellSet):
        values: AdvectionState = cells.values
        nx, ny, num_fields = values.shape
        fields = values.fields
        
        B = np.zeros((nx, ny, num_fields, num_fields, self.DIMENSIONALITY))
        
        B[..., fields.phi, fields.phi, : ] = self.V
        
        return B

`Scheme`¶

cell_scheme

We just need to implement the

$\mathbf{G}_{\frac{1}{2}}(\mathbf{q}) = \left| \mathbf{q}\hat{\mathbf{n}} \right|_f S_f \Rightarrow \left| \varphi \hat{\mathbf{n}}\right|_f S_f$

from josie.scheme.nonconservative import NonConservativeScheme

NonConservativeScheme??

class Upwind(NonConservativeScheme):
    def G(self, cells: MeshCellSet, neighs: CellSet):
        V = self.problem.V
        
        # Get the normal velocities (that are equal everywhere) per each cell
        Vn = np.einsum("k,...kl->...l", V, neighs.normals[..., np.newaxis])
        
        values_face = np.zeros_like(cells.values)
        
        # We use the cell value where Vn > 0 
        np.copyto(values_face, cells.values, where=(Vn>0))
        
        # We use the neighbour value otherwise
        np.copyto(values_face, neighs.values, where=(Vn<=0))
        
        # Multiply by the normal 
        valuesn_face = np.einsum("...i,...j->...ij", values_face, neighs.normals)
        
        # Multiply by the surface
        G = valuesn_face * neighs.surfaces[..., np.newaxis, np.newaxis]
        
        return G

`Mesh`¶

We generated a 1D mesh

from josie.boundary import Line

left = Line([0, 0], [0, 1])
bottom = Line([0, 0], [1, 0])
right = Line([1, 0], [1, 1])
top = Line([0, 1], [1, 1])

We apply periodic boundary condition along the x-axis (no BC on y-axis since it’s a 1D simulation)

from josie.bc import make_periodic, Direction

left, right = make_periodic(left, right, Direction.X)
top.bc = None
bottom.bc = None

print(left, right, top, bottom)

<josie.boundary.boundary.Line object at 0x7f30ec5cf3a0> <josie.boundary.boundary.Line object at 0x7f30ec5cf340> <josie.boundary.boundary.Line object at 0x7f30ec5cf580> <josie.boundary.boundary.Line object at 0x7f30ec5cf460>

from josie.mesh import Mesh
from josie.mesh.cell import SimpleCell
mesh = Mesh(left, bottom, right, top, SimpleCell)
mesh.interpolate(300, 1)
mesh.generate()

`Solver`¶

Let’s assemble our scheme, that still needs the time scheme, and a CFL method

from josie.general.schemes.time import ExplicitEuler

class AdvectionScheme(Upwind, ExplicitEuler):
    def CFL(
        self,
        cells: MeshCellSet,
        CFL_value: float,
    ) -> float:

        U_abs = np.linalg.norm(self.problem.V)
        dx = np.min(cells.surfaces)

        return CFL_value * dx / U_abs


scheme = AdvectionScheme(AdvectionProblem())

from josie.solver import Solver

solver = Solver(mesh, AdvectionState, scheme)

We need now to define an initialization function to initialize the domain

def init_fun(cells: MeshCellSet):
    xc = cells.centroids[..., 0]

    xc_r = np.where(xc >= 0.5)
    xc_l = np.where(xc < 0.5)

    cells.values[xc_r[0], xc_r[1], ...] = 1
    cells.values[xc_l[0], xc_l[1], ...] = 0

Run the simulation¶

Let’s initialize the solver state

solver.init(init_fun)

# Plotting stuff
import matplotlib.pyplot as plt
import copy

from matplotlib.animation import ArtistAnimation
solution = []

t = 0
final_time = 1
CFL = 0.5

Let’s iterate in time choosing a time-based writing strategy (save every $0.01s$ ) to memory

from josie.io.write.strategy import TimeStrategy
from josie.io.write.writer import MemoryWriter 

strategy = TimeStrategy(dt_save=0.01)
writer = MemoryWriter(strategy, solver, final_time, CFL)
writer.solve()

fig, ax = plt.subplots()
ims = []

for solution in writer.data:
    cells = solution.mesh.cells
    x = cells.centroids[..., 0]
    x = x.reshape(x.size)
    phi = cells.values[..., AdvectionFields.phi]
    (im,) = ax.plot(x, phi, 'ko-')
    ims.append([im])
    
    
ani = ArtistAnimation(fig, ims, interval=50)

png

Implementing a solver for the advection equation¶

State¶

Problem¶

Scheme¶

Mesh¶

Solver¶