#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 7 11:17:17 2024
@author: nchevaug
"""
import time
import numpy as np
import matplotlib.pyplot as plt
import h5py
USE_NUMBA = True
#import os
#import pyxfem.postpro.saveloaddata
import pyxfem.core.graph as graph
from pyxfem.core.exec_tools import Profile
from pyxfem.fluidmechanics.exact_riemann_solver import PerfectGaz, ExactRiemannSolver, init_shock_tube
[docs]
def Cjr_1d(edgescoords):
# In 1d the Cjr are constant ...
L = edgescoords[...,1] - edgescoords[...,0]
nedges = len(edgescoords)
C = np.ones((nedges, 2))
C[:,0] = -1.
#C = np.stack([-np.ones(nedges), np.ones(nedges)], axis =1)
return L, C
[docs]
class GlassExplicitEuler_Lagrangian1D:
''' 1d gas dynamic Euler equation using a lagrangian mesh'''
def __init__(self, x, edges, boundary_conditions, gamma = 1.4):
self.gaz = PerfectGaz(gamma)
self.x = x
self.e2v = edges
self.e2v_flat = edges.flatten()
self.v2e = graph.upward_adjacencies(edges)
self.ncells = edges.shape[0]
self.nverts = x.shape[0]
self.bcs = boundary_conditions
self.L, self.C = Cjr_1d(x[edges]) #volumes(lenght), Cjr = dLj/dxr
self.m = None
self.mu = None
self.mE = None
self.t = 0.
self.rho = np.zeros(self.ncells)
self.u = np.zeros(self.ncells)
self.p = np.zeros(self.ncells)
self.c = np.zeros(self.ncells)
[docs]
def set_initial_state(self, density, velocity, pressure):
''' input : density, velocity and pressure are function of x
the state of each cell : in array of Ncells line
( m : masse,
mu quantity of movment (impulsion)
mE : total energy
)
'''
xedges = self.x[self.e2v]
xcog = np.sum(xedges, axis =1)/2.
densitycog = density(xcog)
velocitycog = velocity(xcog)
pressurecog = pressure(xcog)
#mass of each cell :
self.m = densitycog*self.L
#impulsion of each cell :
self.mu = velocitycog*self.m
#kinetic energy of each cell :
ec = 0.5*self.m*(velocitycog)**2
#internal energy of each cell :
ei = self.m*self.gaz.internal_specific_energy(pressurecog,densitycog) # internal energy of each cell
#Total energy of each cell:
self.mE = ei+ec
# set natural variables
self.set_rho_u_p_c()
[docs]
def set_rho_u_p_c(self):
''' return density velocity and pressure for each cell '''
m = self.m
np.divide(self.m, self.L, out = self.rho)
np.divide(self.mu, self.m, out = self.u)
self.gaz.pressure(self.mE/m-0.5*self.u**2, self.rho, out= self.p)
self.gaz.soundspeed(self.p, self.rho, out= self.c)
return
#return self.rho, self.u, self.p, self.c
[docs]
def critical_timestep(self):
return np.min(self.L/self.c)
[docs]
def riemannSolver(self):
# return nodal velocities u_r and flux F_{jr}
nv = self.nverts
edges = self.e2v
# every variable whose name start with c is an array indexed by cell number
cvC = self.C #cell * 2
crho, cu, cp, cc = self.rho, self.u, self.p, self.c
#cvA = (cc*crho)[..., np.newaxis]*np.ones(cvC.shape)
cvA = (cc*crho)[..., np.newaxis]*np.array([1.,1.])[np.newaxis]
vA = np.zeros((nv))
vb = np.zeros((nv))
np.add.at(vA, self.e2v_flat, cvA.flatten())
np.add.at(vb, self.e2v_flat, (cp[:, np.newaxis]*cvC + (crho*cc*cu)[:, np.newaxis]).flatten() )
#correction for external pressure :
def setbc():
for bc in self.bcs :
if bc["type"] !="external_pressure":
raise NotImplementedError("Only external_pressure is currently implemented")
p = bc["value"]
vids = bc["nodes"]
for vid in vids :
cells = self.v2e[vid]
cell = cells[0]
if len(cells) != 1 :
raise ValueError(f"node {vid} is not on a boundary")
cv = 0 if self.e2v[cell,0] == vid else 1
vb[vid]-= cvC[cell, cv]*p
setbc()
vu = vb/vA
cvn = cvC
#cvp = np.stack([cp, cp], axis =1) - (crho*cc)[:, np.newaxis]*(vu[edges]-cu[:,np.newaxis])* cvn
cvp = cp[:, np.newaxis] - (crho*cc)[:, np.newaxis]*(vu[edges]-cu[:,np.newaxis])* cvn
return vu, cvp
[docs]
def step(self, dt):
cvC = self.C
vu, cvp = self.riemannSolver()
cvu = vu[self.e2v]
qflux = np.einsum('jr, jr -> j', cvC, cvp )
eflux = np.einsum('jr, jr, jr -> j', cvC, cvu, cvp)
# tmp = cvC*cvp
# qflux = np.sum(tmp, axis =1)
# eflux = np.sum(tmp*cvu, axis =1)
#cellfluxsum= np.stack([ qflux, eflux], axis = 1)
self.mu -= dt*qflux
self.mE -= dt*eflux
self.x += dt*vu
#self.L += dt* np.einsum('jr, jr -> j', cvC, vu[self.edges])
#self.L , self.C = Cjr_1d(self.x[self.e2v])
#in 1d, cjr i s constant
ex = self.x[self.e2v]
self.L = ex[...,1] - ex[..., 0]
#print('error l', np.linalg.norm(self.L-self.L2))
self.t+=dt
self.set_rho_u_p_c()
#self.gaz.soundspeed(self.p, self.rho, out = self.c)
[docs]
def save_model(self, filepath):
with h5py.File(filepath, 'w') as f:
group = f.create_group("simulation_data")
group.create_dataset("x0", data=self.x)
group.create_dataset("edges", data=self.e2v)
group.create_dataset("gamma", data=self.gaz.gamma)
#group.create_dataset("bcs", data=self.bcs)
[docs]
def save_step(self, filepath, timestep):
timestepname = f't_{timestep:06d}'
data = {'time':self.t,
'x':self.x,
'L':self.L,
'm':self.m,
'mu':self.mu,
'mE':self.mE}
postpro.saveloaddata.save_timestep_data(filepath, timestepname, data)
[docs]
def load_step(self, filepath, timestep):
time_step_name = f't_{timestep:06d}'
data = postpro.saveloaddata.load_timestep_data(filepath, time_step_name)
self.t = data['time']
self.x = data['x']
self.L = data['L']
self.m = data['m']
self.mu = data['mu']
self.mE = data['mE']
[docs]
def load1d_fromfile(filepath):
with h5py.File(filepath, 'r') as f:
if "simulation_data" not in f:
print(f"simulation data do not exist in {filepath}")
return None
group = f["simulation_data"]
x0 = group["x0"][()]
edges = group["edges"][()]
gamma = group["gamma"][()]
#boundary_nodes = group["boundarynodes"][()]
return GlassExplicitEuler_Lagrangian1D(x0, edges, boundary_nodes)
return None
[docs]
def shock_tube1d(pb ="sod", N=50):
gamma = 1.4
L = 1.
x = np.linspace(-L/2., L/2, N+1)
edges = np.stack([np.arange(N), np.arange(1,N+1)], axis =1)
(rhol, ul, pl), (rhor, ur, pr), tend = init_shock_tube(pb)
exactriem = ExactRiemannSolver(gamma = gamma,
rho0 = rhol, rho1 = rhor,
u0 = ul, u1 = ur,
p0 = pl, p1 = pr)
bc_left = {"type" : "external_pressure",
"value" :pl,
"nodes" : [0]
}
bc_right = {"type" : "external_pressure",
"value" :pr,
"nodes" : [N]
}
bcs= [bc_left, bc_right]
pb = GlassExplicitEuler_Lagrangian1D(x, edges, bcs)
pb.set_initial_state(density = lambda x : np.where (x < 0., rhol, rhor),
velocity = lambda x : np.where (x < 0.,ul, ur),
pressure = lambda x : np.where (x < 0.,pl, pr),
)
fig, axes = plt.subplots(2,2)
def plot():
fields = ["rho", "u", "p", "e"]
xmin = np.min(pb.x)
xmax = np.max(pb.x)
LL = xmax - xmin
xe, rhoe, ue, pe = exactriem.solution_sample(x0 = -xmin, t = pb.t, l = LL, n = 3*N)
xe = xe+xmin
xedges = pb.x[pb.e2v]
x0, x1= xedges[:,0], xedges[:,1]
for field, ax in zip(fields, axes.flatten()):
ax.cla()
ax.set_ylabel(field)
if field == "rho" :
ax.hlines(pb.rho, x0, x1, label = "num")
ax.plot(xe, rhoe, label= "exact")
elif field == "u" :
ax.hlines(pb.u, x0, x1, label = "num")
ax.plot(xe, ue, label= "exact")
elif field == "p" :
ax.hlines(pb.p, x0, x1, label = "num")
ax.plot(xe, pe, label= "exact")
elif field =="e":
enum = pb.gaz.internal_specific_energy(pb.p, pb.rho)
ee = pb.gaz.internal_specific_energy(pe, rhoe)
ax.hlines(enum, x0, x1, label = "num")
ax.plot(xe, ee, label= "exact")
else :
print(field, 'is unknown in {__name__}')
ax.legend()
return pb, exactriem, tend, plot
[docs]
def test1d(pb = "sod", N = 100, animate = True, plot_periode = 10 ):
pb, exacsol, tend, plot = shock_tube1d(pb = pb, N = N)
x = pb.x
i = 0
filepath = "res_choctube1d.hdf5"
pb.save_model(filepath)
if (animate):
plt.ion()
plot()
while pb.t < tend :
Profile.enable()
try:
dt = 0.5*pb.critical_timestep()
except :
print("I can' t comput dt!!!")
raise ValueError()
pb.step(dt)
i+=1
Profile.disable()
if i%plot_periode == 0 :
print(i, dt, pb.t)
if animate :
plot()
plt.pause(0.001)
#PB.save_step(filepath, i)
if animate :
plt.ioff()
plot()
plt.show()
Profile.print_stats()
if __name__ == '__main__' :
start = time.process_time()
test1d(pb= "sod", animate = True, N = 200)
# test1d(pb= "toro1", animate = True, N = 1000, plot_periode=10)
# test1d(pb= "toro2", animate = True)
# test1d(pb= "toro3", animate = True) #failed !
# test1d(pb= "toro4", animate = True)
print('time :', time.process_time()-start )