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?
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?
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 = 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 / 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
/ 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 / 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 / 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 / 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,
phase_env,
debug_plot=False,
) # mol / s , Pa
u_top = -mdot * self.mw / 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 + self.g * self.pos)
+ rho_m * e_m_dot / self.mw
)
- rho_m_dot * e_m_dot / self.mw * 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,
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,
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]
