<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;color:black;mso-ligatures:none">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,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none">Eirik H�ydalsvik</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none">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.</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none">t 0 of 1 with dt = 0.2</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none">0 TS dt 0.2 time 0.</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 2.000e-01 retrying with dt=5.000e-02 </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 5.000e-02 retrying with dt=1.250e-02 </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.250e-02 retrying with dt=3.125e-03 </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 3.125e-03 retrying with dt=7.813e-04 </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 7.813e-04 retrying with dt=1.953e-04 </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.953e-04 retrying with dt=4.883e-05 </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 4.883e-05 retrying with dt=1.221e-05 </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.221e-05 retrying with dt=3.052e-06 </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 3.052e-06 retrying with dt=7.629e-07 </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 7.629e-07 retrying with dt=1.907e-07 </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.907e-07 retrying with dt=4.768e-08 </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none">Traceback (most recent call last):</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> File "/Users/iept1445/Code/tank_model/closed_tank.py", line 200, in <module></span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> return_dict1d = get_tank_composition_1d(tank_params)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> File "/Users/iept1445/Code/tank_model/src/tank_model1d.py", line 223, in get_tank_composition_1d</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ts.solve(u=x)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> File "petsc4py/PETSc/TS.pyx", line 2478, in petsc4py.PETSc.TS.solve</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none">petsc4py.PETSc.Error: error code 91</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none">[0] TSSolve() at /Users/iept1445/.cache/uv/sdists-v6/pypi/petsc/3.22.1/svV0XnlR4s3LB4lmsaKjR/src/src/ts/interface/ts.c:4072</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none">[0] TSStep() at /Users/iept1445/.cache/uv/sdists-v6/pypi/petsc/3.22.1/svV0XnlR4s3LB4lmsaKjR/src/src/ts/interface/ts.c:3440</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none">[0] TSStep has failed due to DIVERGED_STEP_REJECTED<br>
<br>
Options for solver:<br>
<br>
COMM = PETSc.COMM_WORLD</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> da = PETSc.DMDA().create(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> dim=(N_vertical,),</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> dof=3,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> stencil_type=PETSc.DMDA().StencilType.STAR,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> stencil_width=1,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # boundary_type=PETSc.DMDA.BoundaryType.GHOSTED,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> )</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> x = da.createGlobalVec()</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> x_old = da.createGlobalVec()</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f = da.createGlobalVec()</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> J = da.createMat()</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> rho_ref = rho_m[0] # kg/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> e_ref = e_m[0] # J/mol</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> p_ref = p0 # Pa</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> x.setArray(np.array([rho_m / rho_ref, e_m / e_ref, ux_m]).T.flatten())</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> x_old.setArray(np.array([rho_m / rho_ref, e_m / e_ref, ux_m]).T.flatten())</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB = PETSc.Options()</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["snes_lag_preconditioner_persists"] = False</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["snes_lag_jacobian"] = 1</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["snes_lag_jacobian_persists"] = False</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["snes_lag_preconditioner"] = 1</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none">optsDB["ksp_type"] = "gmres" # "gmres" # gmres"</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="font-size:11.0pt;color:black;mso-ligatures:none">optsDB["pc_type"] = "ilu" # "lu" # "ilu"</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["snes_type"] = "newtonls"</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["ksp_rtol"] = 1e-7</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["ksp_atol"] = 1e-7</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["ksp_max_it"] = 100</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["snes_rtol"] = 1e-5</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["snes_atol"] = 1e-5</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["snes_stol"] = 1e-5</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["snes_max_it"] = 100</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["snes_mf"] = False</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["ts_max_time"] = t_end</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none">optsDB["ts_type"] = "beuler" # "bdf" #</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="font-size:11.0pt;color:black;mso-ligatures:none">optsDB["ts_max_snes_failures"] = -1</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["ts_monitor"] = ""</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["ts_adapt_monitor"] = ""</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # optsDB["snes_monitor"] = ""</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # optsDB["ksp_monitor"] = ""</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> optsDB["ts_atol"] = 1e-4</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> x0 = x_old</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> residual_wrap = residual_ts(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> eos,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> x0,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> N_vertical,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> g,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> pos,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> z,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> mw,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> dt,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> dx,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> p_amb,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> A_nozzle,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none">r_tank_inner,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none"> mph_uv_flsh_L,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="font-size:11.0pt;color:black;mso-ligatures:none">rho_ref,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> e_ref,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> p_ref,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> closed_tank,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> J,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> da,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> drift_func,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> T_wall,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> tank_params,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> )</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # optsDB["ts_adapt_type"] = "none"</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ts = PETSc.TS().create(comm=COMM)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # TODO: Figure out why DM crashes the code</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # ts.setDM(residual_wrap.da)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ts.setIFunction(residual_wrap.residual_ts, None)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ts.setTimeStep(dt)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ts.setMaxSteps(-1)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ts.setTime(t_start) # s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ts.setMaxTime(t_end) # s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ts.setMaxSteps(1e5)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ts.setStepLimits(1e-3, 1e5)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ts.setFromOptions()</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ts.solve(u=x)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"><br>
<br>
Residual function:<br>
<br>
class residual_ts:</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> def __init__(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> eos,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> x0,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> N,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> g,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> pos,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> z,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> mw,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> dt,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> dx,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> p_amb,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> A_nozzle,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none">r_tank_inner,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none"> mph_uv_flsh_l,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="font-size:11.0pt;color:black;mso-ligatures:none">rho_ref,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> e_ref,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> p_ref,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> closed_tank,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> J,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> da,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> drift_func,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> T_wall,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> tank_params,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ):</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.eos = eos</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.x0 = x0</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.N = N</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.g = g</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.pos = pos</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.z = z</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.mw = mw</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none">self.dt = dt</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none"> self.dx = dx</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="font-size:11.0pt;color:black;mso-ligatures:none">self.p_amb = p_amb</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.A_nozzle = A_nozzle</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none">self.r_tank_inner = r_tank_inner</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none"> self.mph_uv_flsh_L = mph_uv_flsh_l</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="font-size:11.0pt;color:black;mso-ligatures:none">self.rho_ref = rho_ref</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.e_ref = e_ref</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.p_ref = p_ref</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.closed_tank = closed_tank</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none">self.J = J</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none"> self.f = f</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none"> self.da = da</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="font-size:11.0pt;color:black;mso-ligatures:none">self.drift_func = drift_func</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.T_wall = T_wall</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.tank_params = tank_params</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.Q_wall = np.zeros(N)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.n_iter = 0</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.t_current = [0.0]</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.s_top = [0.0]</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.p_choke = [0.0]</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # setting interp func # TODO figure out how to generalize this method</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self._interp_func = _jalla_upwind</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # allocate space for new params</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.p = np.zeros(N) # Pa</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.T = np.zeros(N) # K</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.alpha = np.zeros((2, N))</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho = np.zeros((2, N))</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.e = np.zeros((2, N))</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # allocate space for ghost cells</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.alpha_ghost = np.zeros((2, N + 2))</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_ghost = np.zeros((2, N + 2))</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_m_ghost = np.zeros(N + 2)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.u_m_ghost = np.zeros(N + 1)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.u_ghost = np.zeros((2, N + 1))</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.e_ghost = np.zeros((2, N + 2))</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.pos_ghost = np.zeros(N + 2)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.h_ghost = np.zeros((2, N + 2))</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # allocate soace for local X and Xdot</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.X_LOCAL = da.createLocalVec()</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.XDOT_LOCAL = da.createLocalVec()</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> def residual_ts(self, ts, t, X, XDOT, F):</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.n_iter += 1</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # TODO: Estimate time use</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> """</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> Caculate residuals for equations</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> (rho, rho (e + gx))_t + (rho u, rho u (h + gx))_x = 0</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> P_x = - g \rho</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> """</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> n_phase = 2</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.da.globalToLocal(X, self.X_LOCAL)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.da.globalToLocal(XDOT, self.XDOT_LOCAL)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> x = self.da.getVecArray(self.X_LOCAL)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> xdot = self.da.getVecArray(self.XDOT_LOCAL)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f = self.da.getVecArray(F)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> T_c, v_c, p_c = self.eos.critical(self.z) # K, m3/mol, Pa</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> rho_m = x[:, 0] * self.rho_ref # kg/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> e_m = x[:, 1] * self.e_ref # J/mol</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> u_m = x[:-1, 2] # m/s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # derivatives</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> rho_m_dot = xdot[:, 0] * self.rho_ref # kg/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> e_m_dot = xdot[:, 1] * self.e_ref # kg/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> dt = ts.getTimeStep() # s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> for i in range(self.N):</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # get new parameters</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.mph_uv_flsh_L[i] = self.eos.multi_phase_uvflash(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.z, e_m[i], self.mw / rho_m[i], self.mph_uv_flsh_L[i]</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> )</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> betaL, betaV, betaS = _get_betas(self.mph_uv_flsh_L[i]) # mol/mol</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> beta = [betaL, betaV]</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> if betaS != 0.0:</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> print("there is a solid phase which is not accounted for")</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.T[i], self.p[i] = _get_tank_temperature_pressure(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.mph_uv_flsh_L[i]</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ) # K, Pa)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> for j, phase in enumerate([self.eos.LIQPH, self.eos.VAPPH]):</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # new parameters</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_ghost[:, 1:-1][j][i] = (</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.mw</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> / self.eos.specific_volume(self.T[i], self.p[i], self.z, phase)[0]</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ) # kg/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.e_ghost[:, 1:-1][j][i] = self.eos.internal_energy_tv(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.T[i], self.mw / self.rho_ghost[:, 1:-1][j][i], self.z, phase</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> )[</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> 0</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ] # J/mol</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.h_ghost[:, 1:-1][j][i] = (</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.e_ghost[:, 1:-1][j][i]</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> + self.p[i] * self.mw / self.rho_ghost[:, 1:-1][j][i]</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ) # J/mol</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.alpha_ghost[:, 1:-1][j][i] = (</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> beta[j] / self.rho_ghost[:, 1:-1][j][i] * rho_m[i]</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ) # m3/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # calculate drift velocity</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> for i in range(self.N - 1):</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.u_ghost[:, 1:-1][0][i], self.u_ghost[:, 1:-1][1][i] = (</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> calc_drift_velocity(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> u_m[i],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self._interp_func(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_ghost[:, 1:-1][0][i],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_ghost[:, 1:-1][0][i + 1],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> u_m[i],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ),</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self._interp_func(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_ghost[:, 1:-1][1][i],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_ghost[:, 1:-1][1][i + 1],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> u_m[i],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ),</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.g,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self._interp_func(self.T[i], self.T[i + 1], u_m[i]),</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> T_c,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.r_tank_inner,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self._interp_func(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.alpha_ghost[:, 1:-1][0][i],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.alpha_ghost[:, 1:-1][0][i + 1],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> u_m[i],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ),</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self._interp_func(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.alpha_ghost[:, 1:-1][1][i],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.alpha_ghost[:, 1:-1][1][i + 1],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> u_m[i],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ),</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.drift_func,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> )</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ) # liq m / s , vapour m / s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> u_bottom = 0</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> if self.closed_tank:</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> u_top = 0.0 # m/s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> else:</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # calc phase to skip env_isentrope_cross</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> if (</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.mph_uv_flsh_L[-1].liquid != None</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> and self.mph_uv_flsh_L[-1].vapour == None</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> and self.mph_uv_flsh_L[-1].solid == None</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ):</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> phase_env = self.eos.LIQPH</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> else:</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> phase_env = self.eos.TWOPH</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.h_m = e_m + self.p * self.mw / rho_m # J / mol</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.s_top[0] = _get_s_tank(self.eos, self.mph_uv_flsh_L[-1]) # J / mol / K</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> mdot, self.p_choke[0] = calc_mass_outflow(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.eos,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.z,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.h_m[-1],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.s_top[0],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.p[-1],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.p_amb,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.A_nozzle,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.mw,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> phase_env,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> debug_plot=False,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ) # mol / s , Pa</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> u_top = -mdot * self.mw / rho_m[-1] / (np.pi * self.r_tank_inner**2) # m/s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # assemble vectors with ghost cells</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.alpha_ghost[:, 0] = self.alpha_ghost[:, 1:-1][:, 0] # m3/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.alpha_ghost[:, -1] = self.alpha_ghost[:, 1:-1][:, -1] # m3/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_ghost[:, 0] = self.rho_ghost[:, 1:-1][:, 0] # kg/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_ghost[:, -1] = self.rho_ghost[:, 1:-1][:, -1] # kg/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_m_ghost[0] = rho_m[0] # kg/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_m_ghost[1:-1] = rho_m # kg/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_m_ghost[-1] = rho_m[-1] # kg/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # u_ghost[:, 1:-1] = u # m/s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.u_ghost[:, 0] = u_bottom # m/s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.u_ghost[:, -1] = u_top # m/s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.u_m_ghost[0] = u_bottom # m/s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.u_m_ghost[1:-1] = u_m # m/s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.u_m_ghost[-1] = u_top # m/s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.e_ghost[:, 0] = self.e_ghost[:, 1:-1][:, 0] # J/mol</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.e_ghost[:, -1] = self.e_ghost[:, 1:-1][:, -1] # J/mol</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.pos_ghost[1:-1] = self.pos # m</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.pos_ghost[0] = self.pos[0] # m</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.pos_ghost[-1] = self.pos[-1] # m</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.h_ghost[:, 0] = self.h_ghost[:, 1] # J/mol</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.h_ghost[:, -1] = self.h_ghost[:, -2] # J/mol</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # recalculate wall temperature and heat flux</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # TODO ARE WE DOING THE STAGGERING CORRECTLY?</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> lz = self.tank_params["lz_tank"] / self.N # m</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> if ts.getTime() != self.t_current[0] and self.tank_params["heat_transfer"]:</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.t_current[0] = ts.getTime()</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> for i in range(self.N):</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.T_wall[i], self.Q_wall[i], h_ht = (</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> solve_radial_heat_conduction_implicit(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.tank_params,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.T[i],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.T_wall[i],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> (self.u_m_ghost[i] + self.u_m_ghost[i + 1]) / 2,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_m_ghost[i + 1],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.mph_uv_flsh_L[i],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> lz,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> dt,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> )</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ) # K, J/s, W/m2K</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # Calculate residuals</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f[:, :] = 0.0</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f[:, 0] = dt * rho_m_dot # kg/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f[:-1, 1] = self.p[1:] - self.p[:-1] + self.dx * self.g * rho_m[0:-1] # Pa/m</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f[:, 2] = (</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> dt</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> * (</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> rho_m_dot * (e_m / self.mw + self.g * self.pos)</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> + rho_m * e_m_dot / self.mw</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> )</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> - rho_m_dot * e_m_dot / self.mw * dt**2</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> - self.Q_wall / (np.pi * self.r_tank_inner**2 * lz) * dt</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> ) # J / m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # add contribution from space</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> for i in range(n_phase):</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> e_flux_i = np.zeros_like(self.u_ghost[i]) # J/m3 m/s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> rho_flux_i = np.zeros_like(self.u_ghost[i]) # kg/m2/s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> for j in range(1, self.N + 1):</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> if self.u_ghost[i][j] >= 0.0:</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> rho_flux_new = _rho_flux(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.alpha_ghost[i][j], self.rho_ghost[i][j], self.u_ghost[i][j]</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> )</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> e_flux_new = _e_flux(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.alpha_ghost[i][j],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_ghost[i][j],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.h_ghost[i][j],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.mw,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.g,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.pos_ghost[j],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.u_ghost[i][j],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> )</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # backward euler</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> rho_flux_i[j] = rho_flux_new # kg/m2/s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> e_flux_i[j] = e_flux_new # J/m3 m/s</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> else:</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> rho_flux_new = _rho_flux(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.alpha_ghost[i][j + 1],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_ghost[i][j + 1],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.u_ghost[i][j],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> )</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> e_flux_new = _e_flux(</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.alpha_ghost[i][j + 1],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.rho_ghost[i][j + 1],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.h_ghost[i][j + 1],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.mw,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.g,</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.pos_ghost[j + 1],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> self.u_ghost[i][j],</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> )</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # backward euler</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> rho_flux_i[j] = rho_flux_new</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> e_flux_i[j] = e_flux_new</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # mass eq</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f[:, 0] += (dt / self.dx) * (rho_flux_i[1:] - rho_flux_i[:-1]) # kg/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # energy eq</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f[:, 2] += (dt / self.dx) * (e_flux_i[1:] - e_flux_i[:-1]) # J/m3</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f1_ref, f2_ref, f3_ref = self.rho_ref, self.p_ref, self.e_ref</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f[:, 0] /= f1_ref</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f[:-1, 1] /= f2_ref</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f[:, 2] /= f3_ref</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> # dummy eq</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> f[-1, 1] = x[-1, 2]</span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt;color:black;mso-ligatures:none"> </span><span style="color:black;mso-ligatures:none"><o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
</div>
</body>
</html>