Simulation of a Jet with Euler Equations

import numpy as np

from josie.bc import Dirichlet, Neumann, NeumannDirichlet
from josie.boundary import Line
from josie.euler.eos import PerfectGas
from josie.euler.solver import EulerSolver
from josie.euler.state import EulerState

from josie.mesh import Mesh
from josie.mesh.cell import SimpleCell

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

# EOS
eos = PerfectGas(gamma=1.4)

# Params
JET_CENTER = 0.5
JET_RADIUS = 0.05

# Inlet
U_JET = 1
V_JET = 0
RHO_JET = 1000
P_JET = 1e3

RHOe_JET = eos.rhoe(RHO_JET, P_JET)
e_JET = RHOe_JET / RHO_JET
E_JET = e_JET + 0.5 * (U_JET**2 + V_JET**2)
C_JET = eos.sound_velocity(RHO_JET, P_JET)
print(f"Mach: {U_JET/C_JET}")
Q_JET = EulerState(RHO_JET, RHO_JET*U_JET, RHO_JET*V_JET, RHO_JET*E_JET, RHOe_JET, 
                   U_JET, V_JET, P_JET, C_JET, e_JET)

# Field conditions at init
U_INIT = 0
V_INIT = 0
RHO_INIT = RHO_JET/1000
P_INIT = P_JET
RHOe_INIT = eos.rhoe(RHO_INIT, P_INIT)
e_INIT = RHOe_INIT / RHO_INIT
E_INIT = e_INIT + 0.5 * (U_INIT**2 + V_INIT**2)
C_INIT = eos.sound_velocity(RHO_INIT, P_INIT)

print(f"Mach: {U_INIT/C_INIT}")

Q_INIT = EulerState(RHO_INIT, RHO_INIT*U_INIT, RHO_INIT*V_INIT, RHO_INIT*E_INIT, RHOe_INIT, 
           U_INIT, V_INIT, P_INIT, C_INIT, e_INIT)

# Neumann
dQ = np.zeros(len(EulerState.fields)).view(EulerState)

def partition_fun(centroids: np.ndarray):
    yc = centroids[..., 1].flatten()
    
    # Partition cells of the inlet
    idx = np.where((yc - JET_CENTER)**2 < JET_RADIUS**2)
    
    return idx


# Assign BC to boundaries
left.bc = NeumannDirichlet(dirichlet_value=Q_JET, neumann_value=dQ, partition_fun=partition_fun)
# left.bc = Dirichlet(Q_JET)
top.bc = Dirichlet(Q_INIT)
right.bc = Neumann(dQ)
bottom.bc = Dirichlet(Q_INIT)

mesh = Mesh(left, bottom, right, top, SimpleCell)
mesh.interpolate(25, 100)
mesh.generate()
mesh.write('mesh.xdmf')
# mesh.plot()

Mach: 0.8451542547285166
Mach: 0.0

# Scheme
from josie.euler.schemes import Rusanov
from josie.general.schemes.time import ExplicitEuler

class MyScheme(Rusanov, ExplicitEuler): 
    pass

scheme = MyScheme(eos)



solver = EulerSolver(mesh, scheme)


def init_fun(solver): 
    solver.values[:] = Q_INIT
    

solver.init(init_fun)

from josie.io.write.strategy import TimeStrategy
from josie.io.write.writer import XDMFWriter
writer = XDMFWriter("euler.xdmf", TimeStrategy(dt_save=0.01, animate=True), solver, final_time=1, CFL=0.2)
writer.solve()

solver.show("rho")

/builds/rubendibattista/josiepy/josie/plot/matplotlib.py:224: MatplotlibDeprecationWarning: Starting from Matplotlib 3.6, colorbar() will steal space from the mappable's axes, rather than from the current axes, to place the colorbar.  To silence this warning, explicitly pass the 'ax' argument to colorbar().
  fig.colorbar(patch_coll)