[petsc-users] TS Solver stops working when including ts.setDM
Matthew Knepley
knepley at gmail.com
Tue Feb 25 13:48:02 CST 2025
On Tue, Feb 25, 2025 at 9:51 AM Eirik Jaccheri Høydalsvik <
eirik.hoydalsvik at sintef.no> wrote:
> Hi,
>
> After sending you the email, I rescaled the residual function and got the
> two jacobians to agree down to e-7.
>
> I have tried with “lu” and ”ilu” as preconditioners, and this does not
> work. However, I just tried to use “sor” as a preconditioner, and using sor
> using the da jacobian works just fine!
>
> Why should it work with sor and not with ilu or lu?
>
ILU fails all the time, so that is not surprising. However, I do not
understand why SOR would succeed and LU would fail,
except that SOR is functioning as a kind of globalization by solving very
inexactly. Can you run with
-snes_monitor -snes_converged_reason -ksp_monitor_true_solution
-ksp_converged_reason -snes_linesearch_monitor
and send the output?
Thanks,
Matt
> Eirik
>
> Jacobians:
> row 0: (0, 1.1012) (1, -51356.3) (2, 0.258649) (3, -0.0644364) (4,
> -6402.63) (5, 6.19796e-08)
>
> row 1: (0, -0.445291) (1, 901708.) (2, 0.) (3, 0.44529) (4, -901708.)
> (5, 3.63946e-07)
>
> row 2: (0, 1.10139) (1, -40239.6) (2, 0.258157) (3, -0.0642761) (4,
> -6985.51) (5, 6.19796e-08)
>
> row 3: (0, -0.101197) (1, 51356.3) (2, -0.258649) (3, 1.16563) (4,
> -44953.7) (5, 0.258649) (6, -0.0644364) (7, -6402.63) (8, 8.23293e-08)
>
> row 4: (0, 0.) (1, 0.) (2, 0.) (3, -0.44529) (4, 901708.) (5,
> -3.63946e-07) (6, 0.44529) (7, -901708.) (8, -2.8832e-07)
>
> row 5: (0, -0.101388) (1, 51394.3) (2, -0.258157) (3, 1.16566) (4,
> -33254.1) (5, 0.258157) (6, -0.0642762) (7, -6985.51) (8, 8.27566e-08)
>
> row 6: (3, -0.101197) (4, 51356.3) (5, -0.258649) (6, 1.06444) (7,
> 6402.63) (8, -5.80354e-08)
>
> row 7: (3, 0.) (4, 0.) (5, 0.) (6, 0.) (7, 0.) (8, 1.)
>
> row 8: (3, -0.101388) (4, 51394.3) (5, -0.258157) (6, 1.06428) (7,
> 18140.2) (8, -5.88806e-08)
>
> 1.1011966721737030e+00 -5.1356338141763350e+04 2.5864916418200712e-01
> -6.4436434390614972e-02 -6.4026259743175688e+03 7.0149583230222760e-08
> -1.2114185055821600e-08 7.4938912205780135e-08 -1.2114185055821600e-08
>
> -4.4529045442108389e-01 9.0170829570587131e+05 -3.6394551278146074e-07
> 4.4529038352260736e-01 -9.0170829570607911e+05 -2.8832047116453381e-07
> 2.8832047116453381e-07 -2.8832047116453381e-07 2.8832047116453381e-07
>
> 1.1013882301195626e+00 -4.0239613965471210e+04 2.5815705499526392e-01
> -6.4276195275552644e-02 -6.9855091616075770e+03 7.0994758931791700e-08
> -1.2114185055821600e-08 7.5502362673492770e-08 -1.2114185055821600e-08
>
> -1.0119660581208668e-01 5.1356338141756765e+04 -2.5864909514315410e-01
> 1.1656331102401210e+00 -4.4953712167476042e+04 2.5864905556513557e-01
> -6.4436357255260146e-02 -6.4026259743998889e+03 8.6207799859128616e-08
>
> 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
> -4.4529067184307852e-01 9.0170829570636747e+05 0.0000000000000000e+00
> 4.4529045442108389e-01 -9.0170829570587131e+05 3.6394551278146074e-07
>
> -1.0138834089854361e-01 5.1394313808568040e+04 -2.5815698520133573e-01
> 1.1656640380374783e+00 -3.3254086573823952e+04 2.5815685372183378e-01
> -6.4276095304361430e-02 -6.9855068889361701e+03 8.6288440103988163e-08
>
> -6.9304407528653807e-08 0.0000000000000000e+00 -6.9304407528653807e-08
> -1.0119667848877271e-01 5.1356338141793494e+04 -2.5864912586737532e-01
> 1.0644363666594443e+00 6.4026259743251758e+03 -7.4093736504211182e-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
>
> -6.9586132762510124e-08 0.0000000000000000e+00 -6.9586132762510124e-08
> -1.0138837786069978e-01 5.1394295578534489e+04 -2.5815692427475545e-01
> 1.0642752170084262e+00 1.8140206731963925e+04 -7.4375461738067500e-08
>
>
>
>
> *From: *Matthew Knepley <knepley at gmail.com>
> *Date: *Tuesday, February 25, 2025 at 15:27
> *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 Tue, Feb 25, 2025 at 3:19 AM Eirik Jaccheri Høydalsvik <
> eirik.hoydalsvik at sintef.no> wrote:
>
> 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.
>
>
>
> We usually expect to see single precision accuracy (1e-7), so this
> indicates that your condition number is high.
>
>
>
> If you use LU (-pc_type lu) to solve the linear system, do you get similar
> results?
>
>
>
> Thanks,
>
>
>
> Matt
>
>
>
> 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> 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>
> *Date: *Monday, February 24, 2025 at 14:53
> *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:41 AM Eirik Jaccheri Høydalsvik via petsc-users
> <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 = 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!Z3QqpIz5I2DwEJ5HG6GhoIsIxTOL2irfY5RMgPjrfCc99V9ZTfdxb5k_Tx49NclrnyAR3XQI-OmF1Y9QeBdH$
> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Z3QqpIz5I2DwEJ5HG6GhoIsIxTOL2irfY5RMgPjrfCc99V9ZTfdxb5k_Tx49NclrnyAR3XQI-OmF1ehytAxM$ >
>
>
>
>
> --
>
> 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!Z3QqpIz5I2DwEJ5HG6GhoIsIxTOL2irfY5RMgPjrfCc99V9ZTfdxb5k_Tx49NclrnyAR3XQI-OmF1Y9QeBdH$
> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Z3QqpIz5I2DwEJ5HG6GhoIsIxTOL2irfY5RMgPjrfCc99V9ZTfdxb5k_Tx49NclrnyAR3XQI-OmF1ehytAxM$ >
>
>
>
>
> --
>
> 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!Z3QqpIz5I2DwEJ5HG6GhoIsIxTOL2irfY5RMgPjrfCc99V9ZTfdxb5k_Tx49NclrnyAR3XQI-OmF1Y9QeBdH$
> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Z3QqpIz5I2DwEJ5HG6GhoIsIxTOL2irfY5RMgPjrfCc99V9ZTfdxb5k_Tx49NclrnyAR3XQI-OmF1ehytAxM$ >
>
--
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!Z3QqpIz5I2DwEJ5HG6GhoIsIxTOL2irfY5RMgPjrfCc99V9ZTfdxb5k_Tx49NclrnyAR3XQI-OmF1Y9QeBdH$ <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!Z3QqpIz5I2DwEJ5HG6GhoIsIxTOL2irfY5RMgPjrfCc99V9ZTfdxb5k_Tx49NclrnyAR3XQI-OmF1ehytAxM$ >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250225/46093c72/attachment-0001.html>
More information about the petsc-users
mailing list