[petsc-users] TS Solver stops working when including ts.setDM
Eirik Jaccheri Høydalsvik
eirik.hoydalsvik at sintef.no
Tue Feb 25 02:19:10 CST 2025
Thanks again for the quick response,
I tried prining the jacobians with -ksp_view_mat as you suggested, with a system of only 3 cells (I am studying at a 1d problem). Printing the jacobian in the first timestep I got the two matrices attached at the end of this email. The jacobians are in general agreement, with some small diviations, like the final element of the matrix being 1.6e-5 in the sparse case and 3.7 In the full case.
Questions:
1. Are differences on the order of 1e-5 expected when computing the jacobians in different ways?
2. Do you think these differences can be the cause of my problems? Any suggestions for furtner debugging strategies?
Eirik
! sparse jacobian
row 0: (0, 1.1012) (1, -104.568) (2, 0.258649) (3, -0.0644364) (4, -13.1186) (5, 1.3237e-08)
row 1: (0, -0.44489) (1, 1846.04) (2, 2.12629e-07) (3, 0.445291) (4, -1846.04) (5, 7.08762e-08)
row 2: (0, 540.692) (1, -40219.1) (2, 126.734) (3, -31.5544) (4, -7023.46) (5, 6.48896e-06)
row 3: (0, -0.101197) (1, 104.568) (2, -0.258649) (3, 1.16563) (4, -91.4489) (5, 0.258649) (6, -0.0644365) (7, -13.1186) (8, -4.4809e-08)
row 4: (0, 0.) (1, 0.) (2, 0.) (3, -0.44489) (4, 1846.04) (5, -2.17357e-07) (6, 0.445291) (7, -1846.04) (8, 2.17355e-07)
row 5: (0, -49.7734) (1, 51373.8) (2, -126.734) (3, 572.246) (4, -33195.6) (5, 126.734) (6, -31.5544) (7, -7023.46) (8, -2.19026e-05)
row 6: (3, -0.101197) (4, 104.568) (5, -0.258649) (6, 1.06444) (7, 13.1186) (8, 3.32334e-08)
row 7: (3, 0.) (4, 0.) (5, 0.) (6, 0.) (7, 0.) (8, 1.)
row 8: (3, -49.7734) (4, 51373.8) (5, -126.734) (6, 522.472) (7, 18178.2) (8, 1.61503e-05)
! full jacobian
1.1011966827009450e+00 -1.0456754702270389e+02 2.5864915220241336e-01 -6.4436436239323838e-02 -1.3118626729240630e+01 7.6042484344402957e-08 2.9290438414140398e-08 2.5347494781467651e-08 7.2381179542635411e-08
-4.4488897562431995e-01 1.8460406897256150e+03 -5.0558790784552242e-07 4.4529051871980491e-01 -1.8460406898012175e+03 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
5.4069168701069862e+02 -4.0219094660763396e+04 1.2673402711669499e+02 -3.1554364136687809e+01 -7.0234605760797031e+03 3.7635960251523168e-05 1.4708306305192963e-05 1.2833718246687978e-05 3.5617173111594722e-05
-1.0119659898285956e-01 1.0456754705230254e+02 -2.5864912672499502e-01 1.1656331184446040e+00 -9.1448920317937109e+01 2.5864905777109459e-01 -6.4436443447843730e-02 -1.3118626754633008e+01 3.8772800266974086e-09
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -4.4488904650049199e-01 1.8460406899429699e+03 -2.8823314009443733e-07 4.4529051871980491e-01 -1.8460406898012175e+03 0.0000000000000000e+00
-4.9773392795905657e+01 5.1373794518018862e+04 -1.2673401444613337e+02 5.7224585844509852e+02 -3.3195615874594827e+04 1.2673393520690603e+02 -3.1554356503250116e+01 -7.0234583029005144e+03 1.9105655626471588e-06
-8.1675260962506883e-08 -2.9290438414140398e-08 -2.5347494781467651e-08 -1.0119667997558363e-01 1.0456754704720647e+02 -2.5864913361425051e-01 1.0644364161519400e+00 1.3118626729240630e+01 -7.6042484344402957e-08
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.0000000000000000e+00
-4.0087344635721997e-05 -1.4564107223769502e-05 -1.2401121002417596e-05 -4.9773414819747863e+01 5.1373776293551586e+04 -1.2673397275364130e+02 5.2247224727851687e+02 1.8178158133060850e+04 -3.7347562088676249e-05
From: Matthew Knepley <knepley at gmail.com>
Date: Monday, February 24, 2025 at 15:00
To: Eirik Jaccheri Høydalsvik <eirik.hoydalsvik at sintef.no>
Cc: petsc-users at mcs.anl.gov <petsc-users at mcs.anl.gov>
Subject: Re: [petsc-users] TS Solver stops working when including ts.setDM
On Mon, Feb 24, 2025 at 8:56 AM Eirik Jaccheri Høydalsvik <eirik.hoydalsvik at sintef.no<mailto:eirik.hoydalsvik at sintef.no>> wrote:
1. Thank you for the quick answer, I think this sounds reasonable😊 Is there any way to compare the brute-force jacobian to the one computed using the coloring information?
The easiest way we have is to print them both out:
-ksp_view_mat
on both runs. We have a way to compare the analytic and FD Jacobians (-snes_test_jacobian), but
not two different FDs.
Thanks,
Matt
1.
From: Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>>
Date: Monday, February 24, 2025 at 14:53
To: Eirik Jaccheri Høydalsvik <eirik.hoydalsvik at sintef.no<mailto:eirik.hoydalsvik at sintef.no>>
Cc: petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov> <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>>
Subject: Re: [petsc-users] TS Solver stops working when including ts.setDM
On Mon, Feb 24, 2025 at 8:41 AM Eirik Jaccheri Høydalsvik via petsc-users <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>> wrote:
Hi,
I am using the petsc4py.ts timestepper to solve a PDE. It is not easy to obtain the jacobian for my equations, so I do not provide a jacobian function. The code is given at the end of the email.
When I comment out the function call “ts.setDM(da)”, the code runs and gives reasonable results.
However, when I add this line of code, the program crashes with the error message provided at the end of the email.
Questions:
1. Do you know why adding this line of code can make the SNES solver diverge? Any suggestions for how to debug the issue?
I will not know until I run it, but here is my guess. When the DMDA is specified, PETSc uses coloring to produce the Jacobian. When it is not, it just brute-forces the entire J. My guess is that your residual does not respect the stencil in the DMDA, so the coloring is wrong, making a wrong Jacobian.
2. What is the advantage of adding the DMDA object to the ts solver? Will this speed up the calculation of the finite difference jacobian?
Yes, it speeds up the computation of the FD Jacobian.
Thanks,
Matt
Best regards,
Eirik Høydalsvik
SINTEF ER/NTNU
Error message:
[Eiriks-MacBook-Pro.local:26384] shmem: mmap: an error occurred while determining whether or not /var/folders/w1/35jw9y4n7lsbw0dhjqdhgzz80000gn/T//ompi.Eiriks-MacBook-Pro.501/jf.0/2046361600/sm_segment.Eiriks-MacBook-Pro.501.79f90000.0 could be created.
t 0 of 1 with dt = 0.2
0 TS dt 0.2 time 0.
TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 2.000e-01 retrying with dt=5.000e-02
TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 5.000e-02 retrying with dt=1.250e-02
TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.250e-02 retrying with dt=3.125e-03
TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 3.125e-03 retrying with dt=7.813e-04
TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 7.813e-04 retrying with dt=1.953e-04
TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.953e-04 retrying with dt=4.883e-05
TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 4.883e-05 retrying with dt=1.221e-05
TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.221e-05 retrying with dt=3.052e-06
TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 3.052e-06 retrying with dt=7.629e-07
TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 7.629e-07 retrying with dt=1.907e-07
TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.907e-07 retrying with dt=4.768e-08
Traceback (most recent call last):
File "/Users/iept1445/Code/tank_model/closed_tank.py", line 200, in <module>
return_dict1d = get_tank_composition_1d(tank_params)
File "/Users/iept1445/Code/tank_model/src/tank_model1d.py", line 223, in get_tank_composition_1d
ts.solve(u=x)
File "petsc4py/PETSc/TS.pyx", line 2478, in petsc4py.PETSc.TS.solve
petsc4py.PETSc.Error: error code 91
[0] TSSolve() at /Users/iept1445/.cache/uv/sdists-v6/pypi/petsc/3.22.1/svV0XnlR4s3LB4lmsaKjR/src/src/ts/interface/ts.c:4072
[0] TSStep() at /Users/iept1445/.cache/uv/sdists-v6/pypi/petsc/3.22.1/svV0XnlR4s3LB4lmsaKjR/src/src/ts/interface/ts.c:3440
[0] TSStep has failed due to DIVERGED_STEP_REJECTED
Options for solver:
COMM = PETSc.COMM_WORLD
da = PETSc.DMDA().create(
dim=(N_vertical,),
dof=3,
stencil_type=PETSc.DMDA().StencilType.STAR,
stencil_width=1,
# boundary_type=PETSc.DMDA.BoundaryType.GHOSTED,
)
x = da.createGlobalVec()
x_old = da.createGlobalVec()
f = da.createGlobalVec()
J = da.createMat()
rho_ref = rho_m[0] # kg/m3
e_ref = e_m[0] # J/mol
p_ref = p0 # Pa
x.setArray(np.array([rho_m / rho_ref, e_m / e_ref, ux_m]).T.flatten())
x_old.setArray(np.array([rho_m / rho_ref, e_m / e_ref, ux_m]).T.flatten())
optsDB = PETSc.Options()
optsDB["snes_lag_preconditioner_persists"] = False
optsDB["snes_lag_jacobian"] = 1
optsDB["snes_lag_jacobian_persists"] = False
optsDB["snes_lag_preconditioner"] = 1
optsDB["ksp_type"] = "gmres" # "gmres" # gmres"
optsDB["pc_type"] = "ilu" # "lu" # "ilu"
optsDB["snes_type"] = "newtonls"
optsDB["ksp_rtol"] = 1e-7
optsDB["ksp_atol"] = 1e-7
optsDB["ksp_max_it"] = 100
optsDB["snes_rtol"] = 1e-5
optsDB["snes_atol"] = 1e-5
optsDB["snes_stol"] = 1e-5
optsDB["snes_max_it"] = 100
optsDB["snes_mf"] = False
optsDB["ts_max_time"] = t_end
optsDB["ts_type"] = "beuler" # "bdf" #
optsDB["ts_max_snes_failures"] = -1
optsDB["ts_monitor"] = ""
optsDB["ts_adapt_monitor"] = ""
# optsDB["snes_monitor"] = ""
# optsDB["ksp_monitor"] = ""
optsDB["ts_atol"] = 1e-4
x0 = x_old
residual_wrap = residual_ts(
eos,
x0,
N_vertical,
g,
pos,
z,
mw,
dt,
dx,
p_amb,
A_nozzle,
r_tank_inner,
mph_uv_flsh_L,
rho_ref,
e_ref,
p_ref,
closed_tank,
J,
f,
da,
drift_func,
T_wall,
tank_params,
)
# optsDB["ts_adapt_type"] = "none"
ts = PETSc.TS().create(comm=COMM)
# TODO: Figure out why DM crashes the code
# ts.setDM(residual_wrap.da)
ts.setIFunction(residual_wrap.residual_ts, None)
ts.setTimeStep(dt)
ts.setMaxSteps(-1)
ts.setTime(t_start) # s
ts.setMaxTime(t_end) # s
ts.setMaxSteps(1e5)
ts.setStepLimits(1e-3, 1e5)
ts.setFromOptions()
ts.solve(u=x)
Residual function:
class residual_ts:
def __init__(
self,
eos,
x0,
N,
g,
pos,
z,
mw,
dt,
dx,
p_amb,
A_nozzle,
r_tank_inner,
mph_uv_flsh_l,
rho_ref,
e_ref,
p_ref,
closed_tank,
J,
f,
da,
drift_func,
T_wall,
tank_params,
):
self.eos = eos
self.x0 = x0
self.N = N
self.g = g
self.pos = pos
self.z = z
self.mw<https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlayXlGYI$ > = mw
self.dt = dt
self.dx = dx
self.p_amb = p_amb
self.A_nozzle = A_nozzle
self.r_tank_inner = r_tank_inner
self.mph_uv_flsh_L = mph_uv_flsh_l
self.rho_ref = rho_ref
self.e_ref = e_ref
self.p_ref = p_ref
self.closed_tank = closed_tank
self.J = J
self.f = f
self.da = da
self.drift_func = drift_func
self.T_wall = T_wall
self.tank_params = tank_params
self.Q_wall = np.zeros(N)
self.n_iter = 0
self.t_current = [0.0]
self.s_top = [0.0]
self.p_choke = [0.0]
# setting interp func # TODO figure out how to generalize this method
self._interp_func = _jalla_upwind
# allocate space for new params
self.p = np.zeros(N) # Pa
self.T = np.zeros(N) # K
self.alpha = np.zeros((2, N))
self.rho = np.zeros((2, N))
self.e = np.zeros((2, N))
# allocate space for ghost cells
self.alpha_ghost = np.zeros((2, N + 2))
self.rho_ghost = np.zeros((2, N + 2))
self.rho_m_ghost = np.zeros(N + 2)
self.u_m_ghost = np.zeros(N + 1)
self.u_ghost = np.zeros((2, N + 1))
self.e_ghost = np.zeros((2, N + 2))
self.pos_ghost = np.zeros(N + 2)
self.h_ghost = np.zeros((2, N + 2))
# allocate soace for local X and Xdot
self.X_LOCAL = da.createLocalVec()
self.XDOT_LOCAL = da.createLocalVec()
def residual_ts(self, ts, t, X, XDOT, F):
self.n_iter += 1
# TODO: Estimate time use
"""
Caculate residuals for equations
(rho, rho (e + gx))_t + (rho u, rho u (h + gx))_x = 0
P_x = - g \rho
"""
n_phase = 2
self.da.globalToLocal(X, self.X_LOCAL)
self.da.globalToLocal(XDOT, self.XDOT_LOCAL)
x = self.da.getVecArray(self.X_LOCAL)
xdot = self.da.getVecArray(self.XDOT_LOCAL)
f = self.da.getVecArray(F)
T_c, v_c, p_c = self.eos.critical(self.z) # K, m3/mol, Pa
rho_m = x[:, 0] * self.rho_ref # kg/m3
e_m = x[:, 1] * self.e_ref # J/mol
u_m = x[:-1, 2] # m/s
# derivatives
rho_m_dot = xdot[:, 0] * self.rho_ref # kg/m3
e_m_dot = xdot[:, 1] * self.e_ref # kg/m3
dt = ts.getTimeStep() # s
for i in range(self.N):
# get new parameters
self.mph_uv_flsh_L[i] = self.eos.multi_phase_uvflash(
self.z, e_m[i], self.mw<https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlayXlGYI$ > / rho_m[i], self.mph_uv_flsh_L[i]
)
betaL, betaV, betaS = _get_betas(self.mph_uv_flsh_L[i]) # mol/mol
beta = [betaL, betaV]
if betaS != 0.0:
print("there is a solid phase which is not accounted for")
self.T[i], self.p[i] = _get_tank_temperature_pressure(
self.mph_uv_flsh_L[i]
) # K, Pa)
for j, phase in enumerate([self.eos.LIQPH, self.eos.VAPPH]):
# new parameters
self.rho_ghost[:, 1:-1][j][i] = (
self.mw<https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlayXlGYI$ >
/ self.eos.specific_volume(self.T[i], self.p[i], self.z, phase)[0]
) # kg/m3
self.e_ghost[:, 1:-1][j][i] = self.eos.internal_energy_tv(
self.T[i], self.mw<https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlayXlGYI$ > / self.rho_ghost[:, 1:-1][j][i], self.z, phase
)[
0
] # J/mol
self.h_ghost[:, 1:-1][j][i] = (
self.e_ghost[:, 1:-1][j][i]
+ self.p[i] * self.mw<https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlayXlGYI$ > / self.rho_ghost[:, 1:-1][j][i]
) # J/mol
self.alpha_ghost[:, 1:-1][j][i] = (
beta[j] / self.rho_ghost[:, 1:-1][j][i] * rho_m[i]
) # m3/m3
# calculate drift velocity
for i in range(self.N - 1):
self.u_ghost[:, 1:-1][0][i], self.u_ghost[:, 1:-1][1][i] = (
calc_drift_velocity(
u_m[i],
self._interp_func(
self.rho_ghost[:, 1:-1][0][i],
self.rho_ghost[:, 1:-1][0][i + 1],
u_m[i],
),
self._interp_func(
self.rho_ghost[:, 1:-1][1][i],
self.rho_ghost[:, 1:-1][1][i + 1],
u_m[i],
),
self.g,
self._interp_func(self.T[i], self.T[i + 1], u_m[i]),
T_c,
self.r_tank_inner,
self._interp_func(
self.alpha_ghost[:, 1:-1][0][i],
self.alpha_ghost[:, 1:-1][0][i + 1],
u_m[i],
),
self._interp_func(
self.alpha_ghost[:, 1:-1][1][i],
self.alpha_ghost[:, 1:-1][1][i + 1],
u_m[i],
),
self.drift_func,
)
) # liq m / s , vapour m / s
u_bottom = 0
if self.closed_tank:
u_top = 0.0 # m/s
else:
# calc phase to skip env_isentrope_cross
if (
self.mph_uv_flsh_L[-1].liquid != None
and self.mph_uv_flsh_L[-1].vapour == None
and self.mph_uv_flsh_L[-1].solid == None
):
phase_env = self.eos.LIQPH
else:
phase_env = self.eos.TWOPH
self.h_m = e_m + self.p * self.mw<https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlayXlGYI$ > / rho_m # J / mol
self.s_top[0] = _get_s_tank(self.eos, self.mph_uv_flsh_L[-1]) # J / mol / K
mdot, self.p_choke[0] = calc_mass_outflow(
self.eos,
self.z,
self.h_m[-1],
self.s_top[0],
self.p[-1],
self.p_amb,
self.A_nozzle,
self.mw<https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlayXlGYI$ >,
phase_env,
debug_plot=False,
) # mol / s , Pa
u_top = -mdot * self.mw<https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlayXlGYI$ > / rho_m[-1] / (np.pi * self.r_tank_inner**2) # m/s
# assemble vectors with ghost cells
self.alpha_ghost[:, 0] = self.alpha_ghost[:, 1:-1][:, 0] # m3/m3
self.alpha_ghost[:, -1] = self.alpha_ghost[:, 1:-1][:, -1] # m3/m3
self.rho_ghost[:, 0] = self.rho_ghost[:, 1:-1][:, 0] # kg/m3
self.rho_ghost[:, -1] = self.rho_ghost[:, 1:-1][:, -1] # kg/m3
self.rho_m_ghost[0] = rho_m[0] # kg/m3
self.rho_m_ghost[1:-1] = rho_m # kg/m3
self.rho_m_ghost[-1] = rho_m[-1] # kg/m3
# u_ghost[:, 1:-1] = u # m/s
self.u_ghost[:, 0] = u_bottom # m/s
self.u_ghost[:, -1] = u_top # m/s
self.u_m_ghost[0] = u_bottom # m/s
self.u_m_ghost[1:-1] = u_m # m/s
self.u_m_ghost[-1] = u_top # m/s
self.e_ghost[:, 0] = self.e_ghost[:, 1:-1][:, 0] # J/mol
self.e_ghost[:, -1] = self.e_ghost[:, 1:-1][:, -1] # J/mol
self.pos_ghost[1:-1] = self.pos # m
self.pos_ghost[0] = self.pos[0] # m
self.pos_ghost[-1] = self.pos[-1] # m
self.h_ghost[:, 0] = self.h_ghost[:, 1] # J/mol
self.h_ghost[:, -1] = self.h_ghost[:, -2] # J/mol
# recalculate wall temperature and heat flux
# TODO ARE WE DOING THE STAGGERING CORRECTLY?
lz = self.tank_params["lz_tank"] / self.N # m
if ts.getTime() != self.t_current[0] and self.tank_params["heat_transfer"]:
self.t_current[0] = ts.getTime()
for i in range(self.N):
self.T_wall[i], self.Q_wall[i], h_ht = (
solve_radial_heat_conduction_implicit(
self.tank_params,
self.T[i],
self.T_wall[i],
(self.u_m_ghost[i] + self.u_m_ghost[i + 1]) / 2,
self.rho_m_ghost[i + 1],
self.mph_uv_flsh_L[i],
lz,
dt,
)
) # K, J/s, W/m2K
# Calculate residuals
f[:, :] = 0.0
f[:, 0] = dt * rho_m_dot # kg/m3
f[:-1, 1] = self.p[1:] - self.p[:-1] + self.dx * self.g * rho_m[0:-1] # Pa/m
f[:, 2] = (
dt
* (
rho_m_dot * (e_m / self.mw<https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlayXlGYI$ > + self.g * self.pos)
+ rho_m * e_m_dot / self.mw<https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlayXlGYI$ >
)
- rho_m_dot * e_m_dot / self.mw<https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlayXlGYI$ > * dt**2
- self.Q_wall / (np.pi * self.r_tank_inner**2 * lz) * dt
) # J / m3
# add contribution from space
for i in range(n_phase):
e_flux_i = np.zeros_like(self.u_ghost[i]) # J/m3 m/s
rho_flux_i = np.zeros_like(self.u_ghost[i]) # kg/m2/s
for j in range(1, self.N + 1):
if self.u_ghost[i][j] >= 0.0:
rho_flux_new = _rho_flux(
self.alpha_ghost[i][j], self.rho_ghost[i][j], self.u_ghost[i][j]
)
e_flux_new = _e_flux(
self.alpha_ghost[i][j],
self.rho_ghost[i][j],
self.h_ghost[i][j],
self.mw<https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlayXlGYI$ >,
self.g,
self.pos_ghost[j],
self.u_ghost[i][j],
)
# backward euler
rho_flux_i[j] = rho_flux_new # kg/m2/s
e_flux_i[j] = e_flux_new # J/m3 m/s
else:
rho_flux_new = _rho_flux(
self.alpha_ghost[i][j + 1],
self.rho_ghost[i][j + 1],
self.u_ghost[i][j],
)
e_flux_new = _e_flux(
self.alpha_ghost[i][j + 1],
self.rho_ghost[i][j + 1],
self.h_ghost[i][j + 1],
self.mw<https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlayXlGYI$ >,
self.g,
self.pos_ghost[j + 1],
self.u_ghost[i][j],
)
# backward euler
rho_flux_i[j] = rho_flux_new
e_flux_i[j] = e_flux_new
# mass eq
f[:, 0] += (dt / self.dx) * (rho_flux_i[1:] - rho_flux_i[:-1]) # kg/m3
# energy eq
f[:, 2] += (dt / self.dx) * (e_flux_i[1:] - e_flux_i[:-1]) # J/m3
f1_ref, f2_ref, f3_ref = self.rho_ref, self.p_ref, self.e_ref
f[:, 0] /= f1_ref
f[:-1, 1] /= f2_ref
f[:, 2] /= f3_ref
# dummy eq
f[-1, 1] = x[-1, 2]
--
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener
https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlT6sU48M$ <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlics9Kn4$ >
--
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener
https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlT6sU48M$ <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!YHklSzZKWeUyyqPeA2cK5USsTLG_tjjMs0innpBaQeO9qFz_YqnRf-SHAGLf2GZ-Ku35MtOoNkWCBhU44E22_6AQznSlics9Kn4$ >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250225/4c485e02/attachment-0001.html>
More information about the petsc-users
mailing list