<html xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=Windows-1252">
<meta name="Generator" content="Microsoft Word 15 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:Aptos;
panose-1:2 11 0 4 2 2 2 2 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
font-size:12.0pt;
font-family:"Aptos",sans-serif;
mso-ligatures:standardcontextual;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;
mso-ligatures:none;}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
--></style>
</head>
<body lang="EN-US" link="#467886" vlink="#96607D" style="word-wrap:break-word">
<div class="WordSection1">
<p class="MsoNormal"><span style="font-size:11.0pt">Hi,<br>
<br>
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.<br>
<br>
When I comment out the function call “ts.setDM(da)”, the code runs and gives reasonable results.<br>
<br>
However, when I add this line of code, the program crashes with the error message provided at the end of the email.<br>
<br>
Questions:<br>
<br>
1. Do you know why adding this line of code can make the SNES solver diverge? Any suggestions for how to debug the issue?<br>
<br>
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?<br>
<br>
Best regards,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">Eirik Høydalsvik<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">SINTEF ER/NTNU <br>
<br>
Error message:<br>
<br>
[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.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">t 0 of 1 with dt = 0.2<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">0 TS dt 0.2 time 0.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 2.000e-01 retrying with dt=5.000e-02
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 5.000e-02 retrying with dt=1.250e-02
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.250e-02 retrying with dt=3.125e-03
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 3.125e-03 retrying with dt=7.813e-04
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 7.813e-04 retrying with dt=1.953e-04
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.953e-04 retrying with dt=4.883e-05
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 4.883e-05 retrying with dt=1.221e-05
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.221e-05 retrying with dt=3.052e-06
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 3.052e-06 retrying with dt=7.629e-07
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 7.629e-07 retrying with dt=1.907e-07
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.907e-07 retrying with dt=4.768e-08
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">Traceback (most recent call last):<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> File "/Users/iept1445/Code/tank_model/closed_tank.py", line 200, in <module><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> return_dict1d = get_tank_composition_1d(tank_params)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> File "/Users/iept1445/Code/tank_model/src/tank_model1d.py", line 223, in get_tank_composition_1d<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ts.solve(u=x)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> File "petsc4py/PETSc/TS.pyx", line 2478, in petsc4py.PETSc.TS.solve<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">petsc4py.PETSc.Error: error code 91<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">[0] TSSolve() at /Users/iept1445/.cache/uv/sdists-v6/pypi/petsc/3.22.1/svV0XnlR4s3LB4lmsaKjR/src/src/ts/interface/ts.c:4072<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">[0] TSStep() at /Users/iept1445/.cache/uv/sdists-v6/pypi/petsc/3.22.1/svV0XnlR4s3LB4lmsaKjR/src/src/ts/interface/ts.c:3440<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">[0] TSStep has failed due to DIVERGED_STEP_REJECTED<br>
<br>
Options for solver:<br>
<br>
COMM = PETSc.COMM_WORLD<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> da = PETSc.DMDA().create(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> dim=(N_vertical,),<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> dof=3,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> stencil_type=PETSc.DMDA().StencilType.STAR,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> stencil_width=1,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # boundary_type=PETSc.DMDA.BoundaryType.GHOSTED,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> )<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> x = da.createGlobalVec()<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> x_old = da.createGlobalVec()<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f = da.createGlobalVec()<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> J = da.createMat()<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> rho_ref = rho_m[0] # kg/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> e_ref = e_m[0] # J/mol<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> p_ref = p0 # Pa<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> x.setArray(np.array([rho_m / rho_ref, e_m / e_ref, ux_m]).T.flatten())<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> x_old.setArray(np.array([rho_m / rho_ref, e_m / e_ref, ux_m]).T.flatten())<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB = PETSc.Options()<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["snes_lag_preconditioner_persists"] = False<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["snes_lag_jacobian"] = 1<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["snes_lag_jacobian_persists"] = False<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["snes_lag_preconditioner"] = 1<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> </span><span lang="NO-BOK" style="font-size:11.0pt">optsDB["ksp_type"] = "gmres" # "gmres" # gmres"<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt"> </span><span style="font-size:11.0pt">optsDB["pc_type"] = "ilu" # "lu" # "ilu"<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["snes_type"] = "newtonls"<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["ksp_rtol"] = 1e-7<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["ksp_atol"] = 1e-7<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["ksp_max_it"] = 100<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["snes_rtol"] = 1e-5<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["snes_atol"] = 1e-5<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["snes_stol"] = 1e-5<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["snes_max_it"] = 100<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["snes_mf"] = False<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["ts_max_time"] = t_end<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> </span><span lang="NO-BOK" style="font-size:11.0pt">optsDB["ts_type"] = "beuler" # "bdf" #<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt"> </span><span style="font-size:11.0pt">optsDB["ts_max_snes_failures"] = -1<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["ts_monitor"] = ""<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["ts_adapt_monitor"] = ""<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # optsDB["snes_monitor"] = ""<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # optsDB["ksp_monitor"] = ""<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> optsDB["ts_atol"] = 1e-4<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> x0 = x_old<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> residual_wrap = residual_ts(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> eos,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> x0,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> N_vertical,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> g,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> pos,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> z,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> mw,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> dt,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> dx,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> p_amb,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> A_nozzle,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> </span><span lang="NO-BOK" style="font-size:11.0pt">r_tank_inner,<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt"> mph_uv_flsh_L,<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt"> </span>
<span style="font-size:11.0pt">rho_ref,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> e_ref,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> p_ref,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> closed_tank,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> J,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> da,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> drift_func,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> T_wall,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> tank_params,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> )<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # optsDB["ts_adapt_type"] = "none"<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ts = PETSc.TS().create(comm=COMM)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # TODO: Figure out why DM crashes the code<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # ts.setDM(residual_wrap.da)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ts.setIFunction(residual_wrap.residual_ts, None)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ts.setTimeStep(dt)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ts.setMaxSteps(-1)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ts.setTime(t_start) # s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ts.setMaxTime(t_end) # s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ts.setMaxSteps(1e5)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ts.setStepLimits(1e-3, 1e5)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ts.setFromOptions()<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ts.solve(u=x)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><br>
<br>
Residual function:<br>
<br>
class residual_ts:<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> def __init__(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> eos,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> x0,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> N,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> g,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> pos,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> z,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> mw,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> dt,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> dx,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> p_amb,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> A_nozzle,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> </span><span lang="NO-BOK" style="font-size:11.0pt">r_tank_inner,<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt"> mph_uv_flsh_l,<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt"> </span>
<span style="font-size:11.0pt">rho_ref,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> e_ref,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> p_ref,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> closed_tank,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> J,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> da,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> drift_func,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> T_wall,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> tank_params,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ):<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.eos = eos<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.x0 = x0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.N = N<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.g = g<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.pos = pos<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.z = z<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.mw = mw<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> </span><span lang="NO-BOK" style="font-size:11.0pt">self.dt = dt<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt"> self.dx = dx<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt"> </span>
<span style="font-size:11.0pt">self.p_amb = p_amb<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.A_nozzle = A_nozzle<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> </span><span lang="NO-BOK" style="font-size:11.0pt">self.r_tank_inner = r_tank_inner<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt"> self.mph_uv_flsh_L = mph_uv_flsh_l<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt"> </span>
<span style="font-size:11.0pt">self.rho_ref = rho_ref<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.e_ref = e_ref<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.p_ref = p_ref<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.closed_tank = closed_tank<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> </span><span lang="NO-BOK" style="font-size:11.0pt">self.J = J<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt"> self.f = f<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt"> self.da = da<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt"> </span>
<span style="font-size:11.0pt">self.drift_func = drift_func<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.T_wall = T_wall<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.tank_params = tank_params<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.Q_wall = np.zeros(N)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.n_iter = 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.t_current = [0.0]<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.s_top = [0.0]<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.p_choke = [0.0]<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # setting interp func # TODO figure out how to generalize this method<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self._interp_func = _jalla_upwind<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # allocate space for new params<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.p = np.zeros(N) # Pa<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.T = np.zeros(N) # K<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.alpha = np.zeros((2, N))<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho = np.zeros((2, N))<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.e = np.zeros((2, N))<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # allocate space for ghost cells<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.alpha_ghost = np.zeros((2, N + 2))<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_ghost = np.zeros((2, N + 2))<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_m_ghost = np.zeros(N + 2)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.u_m_ghost = np.zeros(N + 1)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.u_ghost = np.zeros((2, N + 1))<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.e_ghost = np.zeros((2, N + 2))<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.pos_ghost = np.zeros(N + 2)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.h_ghost = np.zeros((2, N + 2))<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # allocate soace for local X and Xdot<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.X_LOCAL = da.createLocalVec()<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.XDOT_LOCAL = da.createLocalVec()<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> def residual_ts(self, ts, t, X, XDOT, F):<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.n_iter += 1<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # TODO: Estimate time use<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> """<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Caculate residuals for equations<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> (rho, rho (e + gx))_t + (rho u, rho u (h + gx))_x = 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> P_x = - g \rho<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> """<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> n_phase = 2<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.da.globalToLocal(X, self.X_LOCAL)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.da.globalToLocal(XDOT, self.XDOT_LOCAL)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> x = self.da.getVecArray(self.X_LOCAL)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> xdot = self.da.getVecArray(self.XDOT_LOCAL)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f = self.da.getVecArray(F)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> T_c, v_c, p_c = self.eos.critical(self.z) # K, m3/mol, Pa<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> rho_m = x[:, 0] * self.rho_ref # kg/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> e_m = x[:, 1] * self.e_ref # J/mol<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> u_m = x[:-1, 2] # m/s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # derivatives<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> rho_m_dot = xdot[:, 0] * self.rho_ref # kg/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> e_m_dot = xdot[:, 1] * self.e_ref # kg/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> dt = ts.getTimeStep() # s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> for i in range(self.N):<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # get new parameters<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.mph_uv_flsh_L[i] = self.eos.multi_phase_uvflash(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.z, e_m[i], self.mw / rho_m[i], self.mph_uv_flsh_L[i]<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> )<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> betaL, betaV, betaS = _get_betas(self.mph_uv_flsh_L[i]) # mol/mol<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> beta = [betaL, betaV]<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> if betaS != 0.0:<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> print("there is a solid phase which is not accounted for")<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.T[i], self.p[i] = _get_tank_temperature_pressure(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.mph_uv_flsh_L[i]<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ) # K, Pa)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> for j, phase in enumerate([self.eos.LIQPH, self.eos.VAPPH]):<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # new parameters<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_ghost[:, 1:-1][j][i] = (<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.mw<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> / self.eos.specific_volume(self.T[i], self.p[i], self.z, phase)[0]<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ) # kg/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.e_ghost[:, 1:-1][j][i] = self.eos.internal_energy_tv(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.T[i], self.mw / self.rho_ghost[:, 1:-1][j][i], self.z, phase<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> )[<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ] # J/mol<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.h_ghost[:, 1:-1][j][i] = (<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.e_ghost[:, 1:-1][j][i]<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> + self.p[i] * self.mw / self.rho_ghost[:, 1:-1][j][i]<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ) # J/mol<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.alpha_ghost[:, 1:-1][j][i] = (<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> beta[j] / self.rho_ghost[:, 1:-1][j][i] * rho_m[i]<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ) # m3/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # calculate drift velocity<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> for i in range(self.N - 1):<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.u_ghost[:, 1:-1][0][i], self.u_ghost[:, 1:-1][1][i] = (<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> calc_drift_velocity(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> u_m[i],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self._interp_func(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_ghost[:, 1:-1][0][i],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_ghost[:, 1:-1][0][i + 1],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> u_m[i],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ),<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self._interp_func(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_ghost[:, 1:-1][1][i],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_ghost[:, 1:-1][1][i + 1],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> u_m[i],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ),<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.g,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self._interp_func(self.T[i], self.T[i + 1], u_m[i]),<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> T_c,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.r_tank_inner,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self._interp_func(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.alpha_ghost[:, 1:-1][0][i],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.alpha_ghost[:, 1:-1][0][i + 1],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> u_m[i],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ),<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self._interp_func(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.alpha_ghost[:, 1:-1][1][i],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.alpha_ghost[:, 1:-1][1][i + 1],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> u_m[i],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ),<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.drift_func,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> )<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ) # liq m / s , vapour m / s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> u_bottom = 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> if self.closed_tank:<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> u_top = 0.0 # m/s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> else:<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # calc phase to skip env_isentrope_cross<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> if (<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.mph_uv_flsh_L[-1].liquid != None<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> and self.mph_uv_flsh_L[-1].vapour == None<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> and self.mph_uv_flsh_L[-1].solid == None<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ):<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> phase_env = self.eos.LIQPH<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> else:<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> phase_env = self.eos.TWOPH<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.h_m = e_m + self.p * self.mw / rho_m # J / mol<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.s_top[0] = _get_s_tank(self.eos, self.mph_uv_flsh_L[-1]) # J / mol / K<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> mdot, self.p_choke[0] = calc_mass_outflow(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.eos,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.z,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.h_m[-1],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.s_top[0],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.p[-1],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.p_amb,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.A_nozzle,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.mw,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> phase_env,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> debug_plot=False,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ) # mol / s , Pa<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> u_top = -mdot * self.mw / rho_m[-1] / (np.pi * self.r_tank_inner**2) # m/s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # assemble vectors with ghost cells<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.alpha_ghost[:, 0] = self.alpha_ghost[:, 1:-1][:, 0] # m3/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.alpha_ghost[:, -1] = self.alpha_ghost[:, 1:-1][:, -1] # m3/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_ghost[:, 0] = self.rho_ghost[:, 1:-1][:, 0] # kg/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_ghost[:, -1] = self.rho_ghost[:, 1:-1][:, -1] # kg/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_m_ghost[0] = rho_m[0] # kg/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_m_ghost[1:-1] = rho_m # kg/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_m_ghost[-1] = rho_m[-1] # kg/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # u_ghost[:, 1:-1] = u # m/s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.u_ghost[:, 0] = u_bottom # m/s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.u_ghost[:, -1] = u_top # m/s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.u_m_ghost[0] = u_bottom # m/s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.u_m_ghost[1:-1] = u_m # m/s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.u_m_ghost[-1] = u_top # m/s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.e_ghost[:, 0] = self.e_ghost[:, 1:-1][:, 0] # J/mol<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.e_ghost[:, -1] = self.e_ghost[:, 1:-1][:, -1] # J/mol<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.pos_ghost[1:-1] = self.pos # m<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.pos_ghost[0] = self.pos[0] # m<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.pos_ghost[-1] = self.pos[-1] # m<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.h_ghost[:, 0] = self.h_ghost[:, 1] # J/mol<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.h_ghost[:, -1] = self.h_ghost[:, -2] # J/mol<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # recalculate wall temperature and heat flux<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # TODO ARE WE DOING THE STAGGERING CORRECTLY?<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> lz = self.tank_params["lz_tank"] / self.N # m<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> if ts.getTime() != self.t_current[0] and self.tank_params["heat_transfer"]:<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.t_current[0] = ts.getTime()<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> for i in range(self.N):<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.T_wall[i], self.Q_wall[i], h_ht = (<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> solve_radial_heat_conduction_implicit(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.tank_params,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.T[i],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.T_wall[i],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> (self.u_m_ghost[i] + self.u_m_ghost[i + 1]) / 2,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_m_ghost[i + 1],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.mph_uv_flsh_L[i],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> lz,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> dt,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> )<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ) # K, J/s, W/m2K<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # Calculate residuals<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f[:, :] = 0.0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f[:, 0] = dt * rho_m_dot # kg/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f[:-1, 1] = self.p[1:] - self.p[:-1] + self.dx * self.g * rho_m[0:-1] # Pa/m<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f[:, 2] = (<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> dt<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> * (<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> rho_m_dot * (e_m / self.mw + self.g * self.pos)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> + rho_m * e_m_dot / self.mw<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> )<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> - rho_m_dot * e_m_dot / self.mw * dt**2<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> - self.Q_wall / (np.pi * self.r_tank_inner**2 * lz) * dt<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ) # J / m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # add contribution from space<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> for i in range(n_phase):<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> e_flux_i = np.zeros_like(self.u_ghost[i]) # J/m3 m/s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> rho_flux_i = np.zeros_like(self.u_ghost[i]) # kg/m2/s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> for j in range(1, self.N + 1):<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> if self.u_ghost[i][j] >= 0.0:<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> rho_flux_new = _rho_flux(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.alpha_ghost[i][j], self.rho_ghost[i][j], self.u_ghost[i][j]<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> )<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> e_flux_new = _e_flux(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.alpha_ghost[i][j],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_ghost[i][j],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.h_ghost[i][j],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.mw,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.g,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.pos_ghost[j],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.u_ghost[i][j],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> )<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # backward euler<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> rho_flux_i[j] = rho_flux_new # kg/m2/s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> e_flux_i[j] = e_flux_new # J/m3 m/s<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> else:<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> rho_flux_new = _rho_flux(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.alpha_ghost[i][j + 1],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_ghost[i][j + 1],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.u_ghost[i][j],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> )<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> e_flux_new = _e_flux(<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.alpha_ghost[i][j + 1],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.rho_ghost[i][j + 1],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.h_ghost[i][j + 1],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.mw,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.g,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.pos_ghost[j + 1],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> self.u_ghost[i][j],<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> )<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # backward euler<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> rho_flux_i[j] = rho_flux_new<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> e_flux_i[j] = e_flux_new<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # mass eq<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f[:, 0] += (dt / self.dx) * (rho_flux_i[1:] - rho_flux_i[:-1]) # kg/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # energy eq<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f[:, 2] += (dt / self.dx) * (e_flux_i[1:] - e_flux_i[:-1]) # J/m3<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f1_ref, f2_ref, f3_ref = self.rho_ref, self.p_ref, self.e_ref<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f[:, 0] /= f1_ref<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f[:-1, 1] /= f2_ref<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f[:, 2] /= f3_ref<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> # dummy eq<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> f[-1, 1] = x[-1, 2]<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
</div>
</body>
</html>