Coverage for pyxfem/fluidmechanics/euler1d_finite_volume_1d.py: 0%
200 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""
4Created on Thu Mar 7 11:17:17 2024
6@author: nchevaug
7"""
8import time
9import numpy as np
10import matplotlib.pyplot as plt
11import h5py
12USE_NUMBA = True
15#import os
16#import pyxfem.postpro.saveloaddata
17import pyxfem.core.graph as graph
18from pyxfem.core.exec_tools import Profile
19from pyxfem.fluidmechanics.exact_riemann_solver import PerfectGaz, ExactRiemannSolver, init_shock_tube
21def Cjr_1d(edgescoords):
22 # In 1d the Cjr are constant ...
23 L = edgescoords[...,1] - edgescoords[...,0]
24 nedges = len(edgescoords)
25 C = np.ones((nedges, 2))
26 C[:,0] = -1.
27 #C = np.stack([-np.ones(nedges), np.ones(nedges)], axis =1)
28 return L, C
30class GlassExplicitEuler_Lagrangian1D:
31 ''' 1d gas dynamic Euler equation using a lagrangian mesh'''
32 def __init__(self, x, edges, boundary_conditions, gamma = 1.4):
33 self.gaz = PerfectGaz(gamma)
34 self.x = x
35 self.e2v = edges
36 self.e2v_flat = edges.flatten()
38 self.v2e = graph.upward_adjacencies(edges)
39 self.ncells = edges.shape[0]
40 self.nverts = x.shape[0]
41 self.bcs = boundary_conditions
42 self.L, self.C = Cjr_1d(x[edges]) #volumes(lenght), Cjr = dLj/dxr
43 self.m = None
44 self.mu = None
45 self.mE = None
46 self.t = 0.
47 self.rho = np.zeros(self.ncells)
48 self.u = np.zeros(self.ncells)
49 self.p = np.zeros(self.ncells)
50 self.c = np.zeros(self.ncells)
53 def set_initial_state(self, density, velocity, pressure):
54 ''' input : density, velocity and pressure are function of x
55 the state of each cell : in array of Ncells line
56 ( m : masse,
57 mu quantity of movment (impulsion)
58 mE : total energy
59 )
60 '''
61 xedges = self.x[self.e2v]
62 xcog = np.sum(xedges, axis =1)/2.
63 densitycog = density(xcog)
64 velocitycog = velocity(xcog)
65 pressurecog = pressure(xcog)
66 #mass of each cell :
67 self.m = densitycog*self.L
68 #impulsion of each cell :
69 self.mu = velocitycog*self.m
70 #kinetic energy of each cell :
71 ec = 0.5*self.m*(velocitycog)**2
72 #internal energy of each cell :
73 ei = self.m*self.gaz.internal_specific_energy(pressurecog,densitycog) # internal energy of each cell
74 #Total energy of each cell:
75 self.mE = ei+ec
76 # set natural variables
77 self.set_rho_u_p_c()
79 def set_rho_u_p_c(self):
80 ''' return density velocity and pressure for each cell '''
81 m = self.m
82 np.divide(self.m, self.L, out = self.rho)
83 np.divide(self.mu, self.m, out = self.u)
84 self.gaz.pressure(self.mE/m-0.5*self.u**2, self.rho, out= self.p)
85 self.gaz.soundspeed(self.p, self.rho, out= self.c)
86 return
87 #return self.rho, self.u, self.p, self.c
89 def critical_timestep(self):
90 return np.min(self.L/self.c)
92 def riemannSolver(self):
93 # return nodal velocities u_r and flux F_{jr}
94 nv = self.nverts
95 edges = self.e2v
96 # every variable whose name start with c is an array indexed by cell number
97 cvC = self.C #cell * 2
98 crho, cu, cp, cc = self.rho, self.u, self.p, self.c
99 #cvA = (cc*crho)[..., np.newaxis]*np.ones(cvC.shape)
100 cvA = (cc*crho)[..., np.newaxis]*np.array([1.,1.])[np.newaxis]
101 vA = np.zeros((nv))
102 vb = np.zeros((nv))
103 np.add.at(vA, self.e2v_flat, cvA.flatten())
104 np.add.at(vb, self.e2v_flat, (cp[:, np.newaxis]*cvC + (crho*cc*cu)[:, np.newaxis]).flatten() )
105 #correction for external pressure :
106 def setbc():
107 for bc in self.bcs :
108 if bc["type"] !="external_pressure":
109 raise NotImplementedError("Only external_pressure is currently implemented")
110 p = bc["value"]
111 vids = bc["nodes"]
112 for vid in vids :
113 cells = self.v2e[vid]
114 cell = cells[0]
115 if len(cells) != 1 :
116 raise ValueError(f"node {vid} is not on a boundary")
117 cv = 0 if self.e2v[cell,0] == vid else 1
118 vb[vid]-= cvC[cell, cv]*p
119 setbc()
120 vu = vb/vA
121 cvn = cvC
122 #cvp = np.stack([cp, cp], axis =1) - (crho*cc)[:, np.newaxis]*(vu[edges]-cu[:,np.newaxis])* cvn
123 cvp = cp[:, np.newaxis] - (crho*cc)[:, np.newaxis]*(vu[edges]-cu[:,np.newaxis])* cvn
125 return vu, cvp
127 def step(self, dt):
128 cvC = self.C
129 vu, cvp = self.riemannSolver()
130 cvu = vu[self.e2v]
131 qflux = np.einsum('jr, jr -> j', cvC, cvp )
132 eflux = np.einsum('jr, jr, jr -> j', cvC, cvu, cvp)
133 # tmp = cvC*cvp
134 # qflux = np.sum(tmp, axis =1)
135 # eflux = np.sum(tmp*cvu, axis =1)
138 #cellfluxsum= np.stack([ qflux, eflux], axis = 1)
139 self.mu -= dt*qflux
140 self.mE -= dt*eflux
141 self.x += dt*vu
142 #self.L += dt* np.einsum('jr, jr -> j', cvC, vu[self.edges])
143 #self.L , self.C = Cjr_1d(self.x[self.e2v])
144 #in 1d, cjr i s constant
145 ex = self.x[self.e2v]
146 self.L = ex[...,1] - ex[..., 0]
148 #print('error l', np.linalg.norm(self.L-self.L2))
149 self.t+=dt
150 self.set_rho_u_p_c()
151 #self.gaz.soundspeed(self.p, self.rho, out = self.c)
153 def save_model(self, filepath):
154 with h5py.File(filepath, 'w') as f:
155 group = f.create_group("simulation_data")
156 group.create_dataset("x0", data=self.x)
157 group.create_dataset("edges", data=self.e2v)
158 group.create_dataset("gamma", data=self.gaz.gamma)
159 #group.create_dataset("bcs", data=self.bcs)
161 def save_step(self, filepath, timestep):
162 timestepname = f't_{timestep:06d}'
163 data = {'time':self.t,
164 'x':self.x,
165 'L':self.L,
166 'm':self.m,
167 'mu':self.mu,
168 'mE':self.mE}
169 postpro.saveloaddata.save_timestep_data(filepath, timestepname, data)
171 def load_step(self, filepath, timestep):
172 time_step_name = f't_{timestep:06d}'
173 data = postpro.saveloaddata.load_timestep_data(filepath, time_step_name)
174 self.t = data['time']
175 self.x = data['x']
176 self.L = data['L']
177 self.m = data['m']
178 self.mu = data['mu']
179 self.mE = data['mE']
181def load1d_fromfile(filepath):
182 with h5py.File(filepath, 'r') as f:
183 if "simulation_data" not in f:
184 print(f"simulation data do not exist in {filepath}")
185 return None
186 group = f["simulation_data"]
187 x0 = group["x0"][()]
188 edges = group["edges"][()]
189 gamma = group["gamma"][()]
190 #boundary_nodes = group["boundarynodes"][()]
191 return GlassExplicitEuler_Lagrangian1D(x0, edges, boundary_nodes)
192 return None
194def shock_tube1d(pb ="sod", N=50):
195 gamma = 1.4
196 L = 1.
197 x = np.linspace(-L/2., L/2, N+1)
198 edges = np.stack([np.arange(N), np.arange(1,N+1)], axis =1)
199 (rhol, ul, pl), (rhor, ur, pr), tend = init_shock_tube(pb)
201 exactriem = ExactRiemannSolver(gamma = gamma,
202 rho0 = rhol, rho1 = rhor,
203 u0 = ul, u1 = ur,
204 p0 = pl, p1 = pr)
205 bc_left = {"type" : "external_pressure",
206 "value" :pl,
207 "nodes" : [0]
208 }
209 bc_right = {"type" : "external_pressure",
210 "value" :pr,
211 "nodes" : [N]
212 }
213 bcs= [bc_left, bc_right]
214 pb = GlassExplicitEuler_Lagrangian1D(x, edges, bcs)
215 pb.set_initial_state(density = lambda x : np.where (x < 0., rhol, rhor),
216 velocity = lambda x : np.where (x < 0.,ul, ur),
217 pressure = lambda x : np.where (x < 0.,pl, pr),
218 )
219 fig, axes = plt.subplots(2,2)
221 def plot():
222 fields = ["rho", "u", "p", "e"]
223 xmin = np.min(pb.x)
224 xmax = np.max(pb.x)
225 LL = xmax - xmin
226 xe, rhoe, ue, pe = exactriem.solution_sample(x0 = -xmin, t = pb.t, l = LL, n = 3*N)
227 xe = xe+xmin
228 xedges = pb.x[pb.e2v]
229 x0, x1= xedges[:,0], xedges[:,1]
230 for field, ax in zip(fields, axes.flatten()):
231 ax.cla()
232 ax.set_ylabel(field)
233 if field == "rho" :
234 ax.hlines(pb.rho, x0, x1, label = "num")
235 ax.plot(xe, rhoe, label= "exact")
236 elif field == "u" :
237 ax.hlines(pb.u, x0, x1, label = "num")
238 ax.plot(xe, ue, label= "exact")
239 elif field == "p" :
240 ax.hlines(pb.p, x0, x1, label = "num")
241 ax.plot(xe, pe, label= "exact")
242 elif field =="e":
243 enum = pb.gaz.internal_specific_energy(pb.p, pb.rho)
244 ee = pb.gaz.internal_specific_energy(pe, rhoe)
245 ax.hlines(enum, x0, x1, label = "num")
246 ax.plot(xe, ee, label= "exact")
247 else :
248 print(field, 'is unknown in {__name__}')
249 ax.legend()
251 return pb, exactriem, tend, plot
253def test1d(pb = "sod", N = 100, animate = True, plot_periode = 10 ):
254 pb, exacsol, tend, plot = shock_tube1d(pb = pb, N = N)
255 x = pb.x
256 i = 0
257 filepath = "res_choctube1d.hdf5"
258 pb.save_model(filepath)
260 if (animate):
261 plt.ion()
262 plot()
263 while pb.t < tend :
264 Profile.enable()
265 try:
266 dt = 0.5*pb.critical_timestep()
267 except :
268 print("I can' t comput dt!!!")
269 raise ValueError()
270 pb.step(dt)
271 i+=1
272 Profile.disable()
273 if i%plot_periode == 0 :
274 print(i, dt, pb.t)
275 if animate :
276 plot()
277 plt.pause(0.001)
278 #PB.save_step(filepath, i)
280 if animate :
281 plt.ioff()
282 plot()
283 plt.show()
284 Profile.print_stats()
286if __name__ == '__main__' :
287 start = time.process_time()
288 test1d(pb= "sod", animate = True, N = 200)
289 # test1d(pb= "toro1", animate = True, N = 1000, plot_periode=10)
290 # test1d(pb= "toro2", animate = True)
291 # test1d(pb= "toro3", animate = True) #failed !
292 # test1d(pb= "toro4", animate = True)
293 print('time :', time.process_time()-start )