On Mon, Feb 24, 2025 at 8:56 AM Eirik Jaccheri Høydalsvik <
eirik.hoydals...@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 <knep...@gmail.com>
> *Date: *Monday, February 24, 2025 at 14:53
> *To: *Eirik Jaccheri Høydalsvik <eirik.hoydals...@sintef.no>
> *Cc: *petsc-users@mcs.anl.gov <petsc-users@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@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 = 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]
>
>
>
>
>
>
>
>
> --
>
> 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!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzok2yPpno$
>
> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzokOWx5lx$
> >
>
--
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!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzok2yPpno$
<https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzokOWx5lx$
>
