<div dir="ltr"><div dir="ltr">On Mon, Feb 24, 2025 at 8:56 AM Eirik Jaccheri Høydalsvik <<a href="mailto:eirik.hoydalsvik@sintef.no">eirik.hoydalsvik@sintef.no</a>> wrote:</div><div class="gmail_quote gmail_quote_container"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="msg-4417416066651576480">
<div lang="EN-US" style="overflow-wrap: break-word;">
<div class="m_3355120538688373402WordSection1">
<ol style="margin-top:0in" start="1" type="1">
<li class="m_3355120538688373402MsoListParagraph" style="margin-left:0in"><span style="font-size:11pt">Thank you for the quick answer, I think this sounds reasonable</span><span style="font-size:11pt;font-family:"Apple Color Emoji"">😊</span><span style="font-size:11pt">
Is there any way to compare the brute-force jacobian to the one computed using the coloring information?</span></li></ol></div></div></div></blockquote><div><br></div><div>The easiest way we have is to print them both out:</div><div><br></div><div> -ksp_view_mat</div><div><br></div><div>on both runs. We have a way to compare the analytic and FD Jacobians (-snes_test_jacobian), but</div><div>not two different FDs.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="msg-4417416066651576480"><div lang="EN-US" style="overflow-wrap: break-word;"><div class="m_3355120538688373402WordSection1"><ol style="margin-top:0in" start="1" type="1"><li class="m_3355120538688373402MsoListParagraph" style="margin-left:0in"> </li></ol>
<div id="m_3355120538688373402mail-editor-reference-message-container">
<div>
<div>
<div style="border-right:none;border-bottom:none;border-left:none;border-top:1pt solid rgb(181,196,223);padding:3pt 0in 0in">
<p class="MsoNormal" style="margin-bottom:12pt"><b><span style="font-size:12pt;color:black">From:
</span></b><span style="font-size:12pt;color:black">Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>><br>
<b>Date: </b>Monday, February 24, 2025 at 14:53<br>
<b>To: </b>Eirik Jaccheri Høydalsvik <<a href="mailto:eirik.hoydalsvik@sintef.no" target="_blank">eirik.hoydalsvik@sintef.no</a>><br>
<b>Cc: </b><a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a> <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>><br>
<b>Subject: </b>Re: [petsc-users] TS Solver stops working when including ts.setDM<u></u><u></u></span></p>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12pt">On Mon, Feb 24, 2025 at 8:41</span><span style="font-size:12pt;font-family:Arial,sans-serif"> </span><span style="font-size:12pt">AM Eirik Jaccheri Høydalsvik via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>>
wrote:<u></u><u></u></span></p>
</div>
<div>
<blockquote style="border-top:none;border-right:none;border-bottom:none;border-left:1pt solid rgb(204,204,204);padding:0in 0in 0in 6pt;margin-left:4.8pt;margin-right:0in">
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:11pt;color:black">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?</span><span style="font-size:12pt"><u></u><u></u></span></p>
</div>
</div>
</div>
</blockquote>
<div>
<p class="MsoNormal"><span style="font-size:12pt"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12pt">I will not know until I run it, but here is my guess. When the DMDA is specified, PETSc uses coloring to produce the Jacobian. When it is not, it just brute-forces the entire J. My guess is that your residual
does not respect the stencil in the DMDA, so the coloring is wrong, making a wrong Jacobian.<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12pt"> <u></u><u></u></span></p>
</div>
<blockquote style="border-top:none;border-right:none;border-bottom:none;border-left:1pt solid rgb(204,204,204);padding:0in 0in 0in 6pt;margin-left:4.8pt;margin-right:0in">
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:11pt;color:black">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?</span><span style="font-size:12pt"><u></u><u></u></span></p>
</div>
</div>
</div>
</blockquote>
<div>
<p class="MsoNormal"><span style="font-size:12pt"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12pt">Yes, it speeds up the computation of the FD Jacobian.<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12pt"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12pt"> Thanks,<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12pt"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12pt"> Matt<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12pt"> <u></u><u></u></span></p>
</div>
<blockquote style="border-top:none;border-right:none;border-bottom:none;border-left:1pt solid rgb(204,204,204);padding:0in 0in 0in 6pt;margin-left:4.8pt;margin-right:0in">
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:11pt;color:black">Best regards,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black">Eirik Høydalsvik</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black">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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black">t 0 of 1 with dt = 0.2</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black">0 TS dt 0.2 time 0.</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> 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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> 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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> 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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> 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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> 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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> 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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> 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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> 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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> 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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> 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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> 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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black">Traceback (most recent call last):</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> File "/Users/iept1445/Code/tank_model/closed_tank.py", line 200, in <module></span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> return_dict1d = get_tank_composition_1d(tank_params)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> File "/Users/iept1445/Code/tank_model/src/tank_model1d.py", line 223, in get_tank_composition_1d</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ts.solve(u=x)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> File "petsc4py/PETSc/TS.pyx", line 2478, in petsc4py.PETSc.TS.solve</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black">petsc4py.PETSc.Error: error code 91</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black">[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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black">[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="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black">[0] TSStep has failed due to DIVERGED_STEP_REJECTED<br>
<br>
Options for solver:<br>
<br>
COMM = PETSc.COMM_WORLD</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> da = PETSc.DMDA().create(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> dim=(N_vertical,),</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> dof=3,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> stencil_type=PETSc.DMDA().StencilType.STAR,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> stencil_width=1,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # boundary_type=PETSc.DMDA.BoundaryType.GHOSTED,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> )</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> x = da.createGlobalVec()</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> x_old = da.createGlobalVec()</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f = da.createGlobalVec()</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> J = da.createMat()</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> rho_ref = rho_m[0] # kg/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> e_ref = e_m[0] # J/mol</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> p_ref = p0 # Pa</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> x.setArray(np.array([rho_m / rho_ref, e_m / e_ref, ux_m]).T.flatten())</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> x_old.setArray(np.array([rho_m / rho_ref, e_m / e_ref, ux_m]).T.flatten())</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB = PETSc.Options()</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["snes_lag_preconditioner_persists"] = False</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["snes_lag_jacobian"] = 1</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["snes_lag_jacobian_persists"] = False</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["snes_lag_preconditioner"] = 1</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span lang="NO-BOK" style="font-size:11pt;color:black">optsDB["ksp_type"] = "gmres" # "gmres" # gmres"</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11pt;color:black"> </span><span style="font-size:11pt;color:black">optsDB["pc_type"] = "ilu" # "lu" # "ilu"</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["snes_type"] = "newtonls"</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["ksp_rtol"] = 1e-7</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["ksp_atol"] = 1e-7</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["ksp_max_it"] = 100</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["snes_rtol"] = 1e-5</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["snes_atol"] = 1e-5</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["snes_stol"] = 1e-5</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["snes_max_it"] = 100</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["snes_mf"] = False</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["ts_max_time"] = t_end</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span lang="NO-BOK" style="font-size:11pt;color:black">optsDB["ts_type"] = "beuler" # "bdf" #</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11pt;color:black"> </span><span style="font-size:11pt;color:black">optsDB["ts_max_snes_failures"] = -1</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["ts_monitor"] = ""</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["ts_adapt_monitor"] = ""</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # optsDB["snes_monitor"] = ""</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # optsDB["ksp_monitor"] = ""</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> optsDB["ts_atol"] = 1e-4</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> x0 = x_old</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> residual_wrap = residual_ts(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> eos,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> x0,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> N_vertical,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> g,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> pos,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> z,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> mw,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> dt,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> dx,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> p_amb,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> A_nozzle,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span lang="NO-BOK" style="font-size:11pt;color:black">r_tank_inner,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11pt;color:black"> mph_uv_flsh_L,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11pt;color:black"> </span><span style="font-size:11pt;color:black">rho_ref,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> e_ref,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> p_ref,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> closed_tank,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> J,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> da,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> drift_func,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> T_wall,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> tank_params,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> )</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # optsDB["ts_adapt_type"] = "none"</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ts = PETSc.TS().create(comm=COMM)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # TODO: Figure out why DM crashes the code</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # ts.setDM(residual_wrap.da)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ts.setIFunction(residual_wrap.residual_ts, None)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ts.setTimeStep(dt)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ts.setMaxSteps(-1)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ts.setTime(t_start) # s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ts.setMaxTime(t_end) # s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ts.setMaxSteps(1e5)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ts.setStepLimits(1e-3, 1e5)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ts.setFromOptions()</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ts.solve(u=x)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"><br>
<br>
Residual function:<br>
<br>
class residual_ts:</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> def __init__(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> eos,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> x0,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> N,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> g,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> pos,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> z,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> mw,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> dt,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> dx,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> p_amb,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> A_nozzle,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span lang="NO-BOK" style="font-size:11pt;color:black">r_tank_inner,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11pt;color:black"> mph_uv_flsh_l,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11pt;color:black"> </span><span style="font-size:11pt;color:black">rho_ref,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> e_ref,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> p_ref,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> closed_tank,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> J,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> da,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> drift_func,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> T_wall,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> tank_params,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ):</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.eos = eos</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.x0 = x0</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.N = N</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.g = g</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.pos = pos</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.z = z</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzogOW71Eq$" target="_blank">self.mw</a> = mw</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span lang="NO-BOK" style="font-size:11pt;color:black">self.dt = dt</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11pt;color:black"> self.dx = dx</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11pt;color:black"> </span><span style="font-size:11pt;color:black">self.p_amb = p_amb</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.A_nozzle = A_nozzle</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span lang="NO-BOK" style="font-size:11pt;color:black">self.r_tank_inner = r_tank_inner</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11pt;color:black"> self.mph_uv_flsh_L = mph_uv_flsh_l</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11pt;color:black"> </span><span style="font-size:11pt;color:black">self.rho_ref = rho_ref</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.e_ref = e_ref</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.p_ref = p_ref</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.closed_tank = closed_tank</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span lang="NO-BOK" style="font-size:11pt;color:black">self.J = J</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11pt;color:black"> self.f = f</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11pt;color:black"> self.da = da</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span lang="NO-BOK" style="font-size:11pt;color:black"> </span><span style="font-size:11pt;color:black">self.drift_func = drift_func</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.T_wall = T_wall</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.tank_params = tank_params</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.Q_wall = np.zeros(N)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.n_iter = 0</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.t_current = [0.0]</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.s_top = [0.0]</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.p_choke = [0.0]</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # setting interp func # TODO figure out how to generalize this method</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self._interp_func = _jalla_upwind</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # allocate space for new params</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.p = np.zeros(N) # Pa</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.T = np.zeros(N) # K</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.alpha = np.zeros((2, N))</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho = np.zeros((2, N))</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.e = np.zeros((2, N))</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # allocate space for ghost cells</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.alpha_ghost = np.zeros((2, N + 2))</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_ghost = np.zeros((2, N + 2))</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_m_ghost = np.zeros(N + 2)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.u_m_ghost = np.zeros(N + 1)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.u_ghost = np.zeros((2, N + 1))</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.e_ghost = np.zeros((2, N + 2))</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.pos_ghost = np.zeros(N + 2)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.h_ghost = np.zeros((2, N + 2))</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # allocate soace for local X and Xdot</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.X_LOCAL = da.createLocalVec()</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.XDOT_LOCAL = da.createLocalVec()</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> def residual_ts(self, ts, t, X, XDOT, F):</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.n_iter += 1</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # TODO: Estimate time use</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> """</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> Caculate residuals for equations</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> (rho, rho (e + gx))_t + (rho u, rho u (h + gx))_x = 0</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> P_x = - g \rho</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> """</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> n_phase = 2</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.da.globalToLocal(X, self.X_LOCAL)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.da.globalToLocal(XDOT, self.XDOT_LOCAL)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> x = self.da.getVecArray(self.X_LOCAL)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> xdot = self.da.getVecArray(self.XDOT_LOCAL)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f = self.da.getVecArray(F)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> T_c, v_c, p_c = self.eos.critical(self.z) # K, m3/mol, Pa</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> rho_m = x[:, 0] * self.rho_ref # kg/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> e_m = x[:, 1] * self.e_ref # J/mol</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> u_m = x[:-1, 2] # m/s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # derivatives</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> rho_m_dot = xdot[:, 0] * self.rho_ref # kg/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> e_m_dot = xdot[:, 1] * self.e_ref # kg/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> dt = ts.getTimeStep() # s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> for i in range(self.N):</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # get new parameters</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.mph_uv_flsh_L[i] = self.eos.multi_phase_uvflash(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.z, e_m[i], <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzogOW71Eq$" target="_blank">self.mw</a> / rho_m[i], self.mph_uv_flsh_L[i]</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> )</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> betaL, betaV, betaS = _get_betas(self.mph_uv_flsh_L[i]) # mol/mol</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> beta = [betaL, betaV]</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> if betaS != 0.0:</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> print("there is a solid phase which is not accounted for")</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.T[i], self.p[i] = _get_tank_temperature_pressure(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.mph_uv_flsh_L[i]</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ) # K, Pa)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> for j, phase in enumerate([self.eos.LIQPH, self.eos.VAPPH]):</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # new parameters</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_ghost[:, 1:-1][j][i] = (</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzogOW71Eq$" target="_blank">self.mw</a></span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> / self.eos.specific_volume(self.T[i], self.p[i], self.z, phase)[0]</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ) # kg/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.e_ghost[:, 1:-1][j][i] = self.eos.internal_energy_tv(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.T[i], <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzogOW71Eq$" target="_blank">self.mw</a> / self.rho_ghost[:, 1:-1][j][i], self.z, phase</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> )[</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> 0</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ] # J/mol</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.h_ghost[:, 1:-1][j][i] = (</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.e_ghost[:, 1:-1][j][i]</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> + self.p[i] * <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzogOW71Eq$" target="_blank">self.mw</a> / self.rho_ghost[:, 1:-1][j][i]</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ) # J/mol</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.alpha_ghost[:, 1:-1][j][i] = (</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> beta[j] / self.rho_ghost[:, 1:-1][j][i] * rho_m[i]</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ) # m3/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # calculate drift velocity</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> for i in range(self.N - 1):</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.u_ghost[:, 1:-1][0][i], self.u_ghost[:, 1:-1][1][i] = (</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> calc_drift_velocity(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> u_m[i],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self._interp_func(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_ghost[:, 1:-1][0][i],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_ghost[:, 1:-1][0][i + 1],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> u_m[i],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ),</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self._interp_func(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_ghost[:, 1:-1][1][i],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_ghost[:, 1:-1][1][i + 1],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> u_m[i],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ),</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.g,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self._interp_func(self.T[i], self.T[i + 1], u_m[i]),</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> T_c,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.r_tank_inner,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self._interp_func(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.alpha_ghost[:, 1:-1][0][i],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.alpha_ghost[:, 1:-1][0][i + 1],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> u_m[i],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ),</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self._interp_func(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.alpha_ghost[:, 1:-1][1][i],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.alpha_ghost[:, 1:-1][1][i + 1],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> u_m[i],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ),</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.drift_func,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> )</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ) # liq m / s , vapour m / s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> u_bottom = 0</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> if self.closed_tank:</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> u_top = 0.0 # m/s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> else:</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # calc phase to skip env_isentrope_cross</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> if (</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.mph_uv_flsh_L[-1].liquid != None</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> and self.mph_uv_flsh_L[-1].vapour == None</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> and self.mph_uv_flsh_L[-1].solid == None</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ):</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> phase_env = self.eos.LIQPH</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> else:</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> phase_env = self.eos.TWOPH</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.h_m = e_m + self.p * <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzogOW71Eq$" target="_blank">self.mw</a> / rho_m # J / mol</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.s_top[0] = _get_s_tank(self.eos, self.mph_uv_flsh_L[-1]) # J / mol / K</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> mdot, self.p_choke[0] = calc_mass_outflow(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.eos,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.z,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.h_m[-1],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.s_top[0],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.p[-1],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.p_amb,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.A_nozzle,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzogOW71Eq$" target="_blank">self.mw</a>,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> phase_env,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> debug_plot=False,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ) # mol / s , Pa</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> u_top = -mdot * <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzogOW71Eq$" target="_blank">self.mw</a> / rho_m[-1] / (np.pi * self.r_tank_inner**2) # m/s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # assemble vectors with ghost cells</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.alpha_ghost[:, 0] = self.alpha_ghost[:, 1:-1][:, 0] # m3/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.alpha_ghost[:, -1] = self.alpha_ghost[:, 1:-1][:, -1] # m3/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_ghost[:, 0] = self.rho_ghost[:, 1:-1][:, 0] # kg/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_ghost[:, -1] = self.rho_ghost[:, 1:-1][:, -1] # kg/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_m_ghost[0] = rho_m[0] # kg/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_m_ghost[1:-1] = rho_m # kg/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_m_ghost[-1] = rho_m[-1] # kg/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # u_ghost[:, 1:-1] = u # m/s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.u_ghost[:, 0] = u_bottom # m/s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.u_ghost[:, -1] = u_top # m/s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.u_m_ghost[0] = u_bottom # m/s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.u_m_ghost[1:-1] = u_m # m/s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.u_m_ghost[-1] = u_top # m/s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.e_ghost[:, 0] = self.e_ghost[:, 1:-1][:, 0] # J/mol</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.e_ghost[:, -1] = self.e_ghost[:, 1:-1][:, -1] # J/mol</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.pos_ghost[1:-1] = self.pos # m</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.pos_ghost[0] = self.pos[0] # m</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.pos_ghost[-1] = self.pos[-1] # m</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.h_ghost[:, 0] = self.h_ghost[:, 1] # J/mol</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.h_ghost[:, -1] = self.h_ghost[:, -2] # J/mol</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # recalculate wall temperature and heat flux</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # TODO ARE WE DOING THE STAGGERING CORRECTLY?</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> lz = self.tank_params["lz_tank"] / self.N # m</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> if ts.getTime() != self.t_current[0] and self.tank_params["heat_transfer"]:</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.t_current[0] = ts.getTime()</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> for i in range(self.N):</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.T_wall[i], self.Q_wall[i], h_ht = (</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> solve_radial_heat_conduction_implicit(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.tank_params,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.T[i],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.T_wall[i],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> (self.u_m_ghost[i] + self.u_m_ghost[i + 1]) / 2,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_m_ghost[i + 1],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.mph_uv_flsh_L[i],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> lz,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> dt,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> )</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ) # K, J/s, W/m2K</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # Calculate residuals</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f[:, :] = 0.0</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f[:, 0] = dt * rho_m_dot # kg/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f[:-1, 1] = self.p[1:] - self.p[:-1] + self.dx * self.g * rho_m[0:-1] # Pa/m</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f[:, 2] = (</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> dt</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> * (</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> rho_m_dot * (e_m / <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzogOW71Eq$" target="_blank">self.mw</a> + self.g * self.pos)</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> + rho_m * e_m_dot / <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzogOW71Eq$" target="_blank">self.mw</a></span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> )</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> - rho_m_dot * e_m_dot / <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzogOW71Eq$" target="_blank">self.mw</a> * dt**2</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> - self.Q_wall / (np.pi * self.r_tank_inner**2 * lz) * dt</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> ) # J / m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # add contribution from space</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> for i in range(n_phase):</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> e_flux_i = np.zeros_like(self.u_ghost[i]) # J/m3 m/s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> rho_flux_i = np.zeros_like(self.u_ghost[i]) # kg/m2/s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> for j in range(1, self.N + 1):</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> if self.u_ghost[i][j] >= 0.0:</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> rho_flux_new = _rho_flux(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.alpha_ghost[i][j], self.rho_ghost[i][j], self.u_ghost[i][j]</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> )</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> e_flux_new = _e_flux(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.alpha_ghost[i][j],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_ghost[i][j],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.h_ghost[i][j],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzogOW71Eq$" target="_blank">self.mw</a>,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.g,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.pos_ghost[j],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.u_ghost[i][j],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> )</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # backward euler</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> rho_flux_i[j] = rho_flux_new # kg/m2/s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> e_flux_i[j] = e_flux_new # J/m3 m/s</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> else:</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> rho_flux_new = _rho_flux(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.alpha_ghost[i][j + 1],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_ghost[i][j + 1],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.u_ghost[i][j],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> )</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> e_flux_new = _e_flux(</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.alpha_ghost[i][j + 1],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.rho_ghost[i][j + 1],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.h_ghost[i][j + 1],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzogOW71Eq$" target="_blank">self.mw</a>,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.g,</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.pos_ghost[j + 1],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> self.u_ghost[i][j],</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> )</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # backward euler</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> rho_flux_i[j] = rho_flux_new</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> e_flux_i[j] = e_flux_new</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # mass eq</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f[:, 0] += (dt / self.dx) * (rho_flux_i[1:] - rho_flux_i[:-1]) # kg/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # energy eq</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f[:, 2] += (dt / self.dx) * (e_flux_i[1:] - e_flux_i[:-1]) # J/m3</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f1_ref, f2_ref, f3_ref = self.rho_ref, self.p_ref, self.e_ref</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f[:, 0] /= f1_ref</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f[:-1, 1] /= f2_ref</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f[:, 2] /= f3_ref</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> # dummy eq</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> f[-1, 1] = x[-1, 2]</span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt;color:black"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:11pt"> </span><span style="font-size:12pt"><u></u><u></u></span></p>
</div>
</div>
</div>
</blockquote>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12pt"><br clear="all">
<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12pt"><u></u> <u></u></span></p>
</div>
<p class="MsoNormal"><span class="m_3355120538688373402gmailsignatureprefix"><span style="font-size:12pt">--
</span></span><span style="font-size:12pt"><u></u><u></u></span></p>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12pt">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12pt"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12pt"><a href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzokOWx5lx$" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><u></u><u></u></span></p>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div></blockquote></div><div><br clear="all"></div><div><br></div><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!dcoz7MDn7GnxkOOW8KrkFC-3TAKZKmVlbtBSOJqC2xpb8AuzPeBKTeVUS-nWxQhzkrKp4wQF9njzokOWx5lx$" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>