<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=utf-8">
<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;}
@font-face
{font-family:"Apple Color Emoji";
panose-1:0 0 0 0 0 0 0 0 0 0;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
font-size:10.0pt;
font-family:"Aptos",sans-serif;}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
p.m-1124721174010954366m-6386679822550478849m3355120538688373402msolistparagraph, li.m-1124721174010954366m-6386679822550478849m3355120538688373402msolistparagraph, div.m-1124721174010954366m-6386679822550478849m3355120538688373402msolistparagraph
{mso-style-name:m_-1124721174010954366m-6386679822550478849m3355120538688373402msolistparagraph;
mso-margin-top-alt:auto;
margin-right:0in;
mso-margin-bottom-alt:auto;
margin-left:0in;
font-size:12.0pt;
font-family:"Aptos",sans-serif;}
span.m-1124721174010954366m-6386679822550478849m3355120538688373402gmailsignatureprefix
{mso-style-name:m_-1124721174010954366m-6386679822550478849m3355120538688373402gmailsignatureprefix;}
span.m-1124721174010954366m-6386679822550478849gmailsignatureprefix
{mso-style-name:m_-1124721174010954366m-6386679822550478849gmailsignatureprefix;}
span.m-1124721174010954366gmailsignatureprefix
{mso-style-name:m_-1124721174010954366gmailsignatureprefix;}
span.gmailsignatureprefix
{mso-style-name:gmail_signature_prefix;}
span.EmailStyle23
{mso-style-type:personal-reply;
font-family:"Aptos",sans-serif;
color:windowtext;}
.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;}
/* List Definitions */
@list l0
{mso-list-id:1351103408;
mso-list-template-ids:741610242;}
@list l1
{mso-list-id:1406412846;
mso-list-template-ids:1653488824;}
ol
{margin-bottom:0in;}
ul
{margin-bottom:0in;}
--></style>
</head>
<body lang="EN-US" link="blue" vlink="purple" style="word-wrap:break-word">
<div class="WordSection1">
<p class="MsoNormal"><span style="font-size:11.0pt">Hi,<br>
<br>
Here is the output when I run with “lu” and the settings you suggested:<br>
<br>
0 TS dt 0.1 time 0.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> 0 SNES Function norm 2.982668991189e-01<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> PC failed due to FACTOR_NUMERIC_ZEROPIVOT
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.000e-01 retrying with dt=2.500e-02
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> 0 SNES Function norm 7.456672477972e-02<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> PC failed due to FACTOR_NUMERIC_ZEROPIVOT
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 2.500e-02 retrying with dt=6.250e-03
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> 0 SNES Function norm 1.864168119493e-02<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> PC failed due to FACTOR_NUMERIC_ZEROPIVOT
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 6.250e-03 retrying with dt=1.563e-03
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> 0 SNES Function norm 4.660420298733e-03<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> PC failed due to FACTOR_NUMERIC_ZEROPIVOT
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 1.563e-03 retrying with dt=3.906e-04
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> 0 SNES Function norm 1.165105074683e-03<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> PC failed due to FACTOR_NUMERIC_ZEROPIVOT
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 3.906e-04 retrying with dt=9.766e-05
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> 0 SNES Function norm 2.912762686708e-04<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> PC failed due to FACTOR_NUMERIC_ZEROPIVOT
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 9.766e-05 retrying with dt=2.441e-05
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> 0 SNES Function norm 7.281906716770e-05<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> PC failed due to FACTOR_NUMERIC_ZEROPIVOT
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 2.441e-05 retrying with dt=6.104e-06
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> 0 SNES Function norm 1.820476679192e-05<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> PC failed due to FACTOR_NUMERIC_ZEROPIVOT
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none step 0 stage rejected (SNES reason DIVERGED_LINEAR_SOLVE) t=0 + 6.104e-06 retrying with dt=1.526e-06
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> 0 SNES Function norm 4.551191697981e-06<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> TSAdapt none beuler 0: step 0 accepted t=0 + 1.526e-06 dt=1.526e-06
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">1 TS dt 1.52588e-06 time 1.52588e-06<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<div id="mail-editor-reference-message-container">
<div>
<div>
<div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0in 0in 0in">
<p class="MsoNormal" style="margin-bottom:12.0pt"><b><span style="font-size:12.0pt;color:black">From:
</span></b><span style="font-size:12.0pt;color:black">Matthew Knepley <knepley@gmail.com><br>
<b>Date: </b>Tuesday, February 25, 2025 at 20:48<br>
<b>To: </b>Eirik Jaccheri Høydalsvik <eirik.hoydalsvik@sintef.no><br>
<b>Cc: </b>petsc-users@mcs.anl.gov <petsc-users@mcs.anl.gov><br>
<b>Subject: </b>Re: [petsc-users] TS Solver stops working when including ts.setDM<o:p></o:p></span></p>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">On Tue, Feb 25, 2025 at 9:51</span><span style="font-size:12.0pt;font-family:"Arial",sans-serif"> </span><span style="font-size:12.0pt">AM Eirik Jaccheri Høydalsvik <<a href="mailto:eirik.hoydalsvik@sintef.no">eirik.hoydalsvik@sintef.no</a>>
wrote:<o:p></o:p></span></p>
</div>
<div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-right:0in">
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;margin-bottom:12.0pt"><span style="font-size:11.0pt">Hi,<br>
<br>
After sending you the email, I rescaled the residual function and got the two jacobians to agree down to e-7.<br>
<br>
I have tried with “lu” and ”ilu” as preconditioners, and this does not work. However, I just tried to use “sor” as a preconditioner, and using sor using the da jacobian works just fine!</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">Why should it work with sor and not with ilu or lu?</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
</div>
</div>
</div>
</blockquote>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">ILU fails all the time, so that is not surprising. However, I do not understand why SOR would succeed and LU would fail,<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">except that SOR is functioning as a kind of globalization by solving very inexactly. Can you run with<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> -snes_monitor -snes_converged_reason -ksp_monitor_true_solution -ksp_converged_reason -snes_linesearch_monitor<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">and send the output?<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> Thanks,<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> Matt<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-right:0in">
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">Eirik<br>
<br>
Jacobians:<br>
row 0: (0, 1.1012) (1, -51356.3) (2, 0.258649) (3, -0.0644364) (4, -6402.63) (5, 6.19796e-08)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 1: (0, -0.445291) (1, 901708.) (2, 0.) (3, 0.44529) (4, -901708.) (5, 3.63946e-07)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 2: (0, 1.10139) (1, -40239.6) (2, 0.258157) (3, -0.0642761) (4, -6985.51) (5, 6.19796e-08)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 3: (0, -0.101197) (1, 51356.3) (2, -0.258649) (3, 1.16563) (4, -44953.7) (5, 0.258649) (6, -0.0644364) (7, -6402.63) (8, 8.23293e-08)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 4: (0, 0.) (1, 0.) (2, 0.) (3, -0.44529) (4, 901708.) (5, -3.63946e-07) (6, 0.44529) (7, -901708.) (8, -2.8832e-07)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 5: (0, -0.101388) (1, 51394.3) (2, -0.258157) (3, 1.16566) (4, -33254.1) (5, 0.258157) (6, -0.0642762) (7, -6985.51) (8, 8.27566e-08)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 6: (3, -0.101197) (4, 51356.3) (5, -0.258649) (6, 1.06444) (7, 6402.63) (8, -5.80354e-08)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 7: (3, 0.) (4, 0.) (5, 0.) (6, 0.) (7, 0.) (8, 1.)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 8: (3, -0.101388) (4, 51394.3) (5, -0.258157) (6, 1.06428) (7, 18140.2) (8, -5.88806e-08)<br>
<br>
1.1011966721737030e+00 -5.1356338141763350e+04 2.5864916418200712e-01 -6.4436434390614972e-02 -6.4026259743175688e+03 7.0149583230222760e-08 -1.2114185055821600e-08 7.4938912205780135e-08 -1.2114185055821600e-08
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">-4.4529045442108389e-01 9.0170829570587131e+05 -3.6394551278146074e-07 4.4529038352260736e-01 -9.0170829570607911e+05 -2.8832047116453381e-07 2.8832047116453381e-07
-2.8832047116453381e-07 2.8832047116453381e-07 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">1.1013882301195626e+00 -4.0239613965471210e+04 2.5815705499526392e-01 -6.4276195275552644e-02 -6.9855091616075770e+03 7.0994758931791700e-08 -1.2114185055821600e-08
7.5502362673492770e-08 -1.2114185055821600e-08 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">-1.0119660581208668e-01 5.1356338141756765e+04 -2.5864909514315410e-01 1.1656331102401210e+00 -4.4953712167476042e+04 2.5864905556513557e-01 -6.4436357255260146e-02
-6.4026259743998889e+03 8.6207799859128616e-08 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -4.4529067184307852e-01 9.0170829570636747e+05 0.0000000000000000e+00 4.4529045442108389e-01
-9.0170829570587131e+05 3.6394551278146074e-07 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">-1.0138834089854361e-01 5.1394313808568040e+04 -2.5815698520133573e-01 1.1656640380374783e+00 -3.3254086573823952e+04 2.5815685372183378e-01 -6.4276095304361430e-02
-6.9855068889361701e+03 8.6288440103988163e-08 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">-6.9304407528653807e-08 0.0000000000000000e+00 -6.9304407528653807e-08 -1.0119667848877271e-01 5.1356338141793494e+04 -2.5864912586737532e-01 1.0644363666594443e+00
6.4026259743251758e+03 -7.4093736504211182e-08 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
0.0000000000000000e+00 1.0000000000000000e+00 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;margin-bottom:12.0pt"><span style="font-size:11.0pt">-6.9586132762510124e-08 0.0000000000000000e+00 -6.9586132762510124e-08 -1.0138837786069978e-01 5.1394295578534489e+04 -2.5815692427475545e-01 1.0642752170084262e+00
1.8140206731963925e+04 -7.4375461738067500e-08<br>
<br>
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<div id="m_-1124721174010954366mail-editor-reference-message-container">
<div>
<div>
<div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0in 0in 0in">
<p class="MsoNormal" style="mso-margin-top-alt:auto;margin-bottom:12.0pt"><b><span style="font-size:12.0pt;color:black">From:
</span></b><span style="font-size:12.0pt;color:black">Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>><br>
<b>Date: </b>Tuesday, February 25, 2025 at 15:27<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</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
</div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt">On Tue, Feb 25, 2025 at 3:19</span><span style="font-size:12.0pt;font-family:"Arial",sans-serif"> </span><span style="font-size:12.0pt">AM Eirik Jaccheri
Høydalsvik <<a href="mailto:eirik.hoydalsvik@sintef.no" target="_blank">eirik.hoydalsvik@sintef.no</a>> wrote:<o:p></o:p></span></p>
</div>
<div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-top:5.0pt;margin-right:0in;margin-bottom:5.0pt">
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">Thanks again for the quick response,<br>
<br>
I tried prining the jacobians with -ksp_view_mat as you suggested, with a system of only 3 cells (I am studying at a 1d problem). Printing the jacobian in the first timestep I got the two matrices attached at the end of this email. The jacobians are in general
agreement, with some small diviations, like the final element of the matrix being 1.6e-5 in the sparse case and 3.7 In the full case.</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
</div>
</div>
</div>
</blockquote>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt">We usually expect to see single precision accuracy (1e-7), so this indicates that your condition number is high.<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt">If you use LU (-pc_type lu) to solve the linear system, do you get similar results?<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> Thanks,<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> Matt<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-top:5.0pt;margin-right:0in;margin-bottom:5.0pt">
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">Questions:<br>
<br>
1. Are differences on the order of 1e-5 expected when computing the jacobians in different ways?<br>
<br>
2. Do you think these differences can be the cause of my problems? Any suggestions for furtner debugging strategies?<br>
<br>
Eirik<br>
<br>
! sparse jacobian</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 0: (0, 1.1012) (1, -104.568) (2, 0.258649) (3, -0.0644364) (4, -13.1186) (5, 1.3237e-08)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 1: (0, -0.44489) (1, 1846.04) (2, 2.12629e-07) (3, 0.445291) (4, -1846.04) (5, 7.08762e-08)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 2: (0, 540.692) (1, -40219.1) (2, 126.734) (3, -31.5544) (4, -7023.46) (5, 6.48896e-06)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 3: (0, -0.101197) (1, 104.568) (2, -0.258649) (3, 1.16563) (4, -91.4489) (5, 0.258649) (6, -0.0644365) (7, -13.1186) (8, -4.4809e-08)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 4: (0, 0.) (1, 0.) (2, 0.) (3, -0.44489) (4, 1846.04) (5, -2.17357e-07) (6, 0.445291) (7, -1846.04) (8, 2.17355e-07)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 5: (0, -49.7734) (1, 51373.8) (2, -126.734) (3, 572.246) (4, -33195.6) (5, 126.734) (6, -31.5544) (7, -7023.46) (8, -2.19026e-05)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 6: (3, -0.101197) (4, 104.568) (5, -0.258649) (6, 1.06444) (7, 13.1186) (8, 3.32334e-08)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 7: (3, 0.) (4, 0.) (5, 0.) (6, 0.) (7, 0.) (8, 1.)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">row 8: (3, -49.7734) (4, 51373.8) (5, -126.734) (6, 522.472) (7, 18178.2) (8, 1.61503e-05)
</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">! full jacobian</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">1.1011966827009450e+00 -1.0456754702270389e+02 2.5864915220241336e-01 -6.4436436239323838e-02 -1.3118626729240630e+01 7.6042484344402957e-08 2.9290438414140398e-08
2.5347494781467651e-08 7.2381179542635411e-08 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">-4.4488897562431995e-01 1.8460406897256150e+03 -5.0558790784552242e-07 4.4529051871980491e-01 -1.8460406898012175e+03 0.0000000000000000e+00 0.0000000000000000e+00
0.0000000000000000e+00 0.0000000000000000e+00 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">5.4069168701069862e+02 -4.0219094660763396e+04 1.2673402711669499e+02 -3.1554364136687809e+01 -7.0234605760797031e+03 3.7635960251523168e-05 1.4708306305192963e-05
1.2833718246687978e-05 3.5617173111594722e-05 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">-1.0119659898285956e-01 1.0456754705230254e+02 -2.5864912672499502e-01 1.1656331184446040e+00 -9.1448920317937109e+01 2.5864905777109459e-01 -6.4436443447843730e-02
-1.3118626754633008e+01 3.8772800266974086e-09 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 -4.4488904650049199e-01 1.8460406899429699e+03 -2.8823314009443733e-07 4.4529051871980491e-01
-1.8460406898012175e+03 0.0000000000000000e+00 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">-4.9773392795905657e+01 5.1373794518018862e+04 -1.2673401444613337e+02 5.7224585844509852e+02 -3.3195615874594827e+04 1.2673393520690603e+02 -3.1554356503250116e+01
-7.0234583029005144e+03 1.9105655626471588e-06 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">-8.1675260962506883e-08 -2.9290438414140398e-08 -2.5347494781467651e-08 -1.0119667997558363e-01 1.0456754704720647e+02 -2.5864913361425051e-01 1.0644364161519400e+00
1.3118626729240630e+01 -7.6042484344402957e-08 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
0.0000000000000000e+00 1.0000000000000000e+00 </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt">-4.0087344635721997e-05 -1.4564107223769502e-05 -1.2401121002417596e-05 -4.9773414819747863e+01 5.1373776293551586e+04 -1.2673397275364130e+02 5.2247224727851687e+02
1.8178158133060850e+04 -3.7347562088676249e-05</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<div id="m_-1124721174010954366m_-6386679822550478849mail-editor-reference-message-container">
<div>
<div>
<div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0in 0in 0in">
<p class="MsoNormal" style="mso-margin-top-alt:auto;margin-bottom:12.0pt"><b><span style="font-size:12.0pt;color:black">From:
</span></b><span style="font-size:12.0pt;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 15:00<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</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
</div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt">On Mon, Feb 24, 2025 at 8:56</span><span style="font-size:12.0pt;font-family:"Arial",sans-serif"> </span><span style="font-size:12.0pt">AM Eirik Jaccheri
Høydalsvik <<a href="mailto:eirik.hoydalsvik@sintef.no" target="_blank">eirik.hoydalsvik@sintef.no</a>> wrote:<o:p></o:p></span></p>
</div>
<div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-top:5.0pt;margin-right:0in;margin-bottom:5.0pt">
<div>
<div>
<div>
<ol start="1" type="1">
<li class="m-1124721174010954366m-6386679822550478849m3355120538688373402msolistparagraph" style="mso-list:l1 level1 lfo1">
<span style="font-size:11.0pt">Thank you for the quick answer, I think this sounds reasonable</span><span style="font-size:11.0pt;font-family:"Apple Color Emoji"">😊</span><span style="font-size:11.0pt"> Is there any way to compare the brute-force jacobian to
the one computed using the coloring information?</span><o:p></o:p></li></ol>
</div>
</div>
</div>
</blockquote>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt">The easiest way we have is to print them both out:<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> -ksp_view_mat<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt">on both runs. We have a way to compare the analytic and FD Jacobians (-snes_test_jacobian), but<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt">not two different FDs.<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> Thanks,<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> Matt<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-top:5.0pt;margin-right:0in;margin-bottom:5.0pt">
<div>
<div>
<div>
<ol start="1" type="1">
<li class="m-1124721174010954366m-6386679822550478849m3355120538688373402msolistparagraph" style="mso-list:l0 level1 lfo2">
<o:p></o:p></li></ol>
<div id="m_-1124721174010954366m_-6386679822550478849m_3355120538688373402mail-editor-reference-message-container">
<div>
<div>
<div style="border:none;border-top:solid #B5C4DF 1.0pt;padding:3.0pt 0in 0in 0in">
<p class="MsoNormal" style="mso-margin-top-alt:auto;margin-bottom:12.0pt"><b><span style="font-size:12.0pt;color:black">From:
</span></b><span style="font-size:12.0pt;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</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
</div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt">On Mon, Feb 24, 2025 at 8:41</span><span style="font-size:12.0pt;font-family:"Arial",sans-serif"> </span><span style="font-size:12.0pt">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:<o:p></o:p></span></p>
</div>
<div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-top:5.0pt;margin-right:0in;margin-bottom:5.0pt">
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
</div>
</div>
</div>
</blockquote>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt">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.<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-top:5.0pt;margin-right:0in;margin-bottom:5.0pt">
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
</div>
</div>
</div>
</blockquote>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt">Yes, it speeds up the computation of the FD Jacobian.<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> Thanks,<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> Matt<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<blockquote style="border:none;border-left:solid #CCCCCC 1.0pt;padding:0in 0in 0in 6.0pt;margin-left:4.8pt;margin-top:5.0pt;margin-right:0in;margin-bottom:5.0pt">
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black">Best regards,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black">Eirik Høydalsvik</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black">t 0 of 1 with dt = 0.2</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black">0 TS dt 0.2 time 0.</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black">Traceback (most recent call last):</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> File "/Users/iept1445/Code/tank_model/closed_tank.py", line 200, in <module></span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> return_dict1d = get_tank_composition_1d(tank_params)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> File "/Users/iept1445/Code/tank_model/src/tank_model1d.py", line 223, in get_tank_composition_1d</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ts.solve(u=x)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> File "petsc4py/PETSc/TS.pyx", line 2478, in petsc4py.PETSc.TS.solve</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black">petsc4py.PETSc.Error: error code 91</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> da = PETSc.DMDA().create(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> dim=(N_vertical,),</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> dof=3,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> stencil_type=PETSc.DMDA().StencilType.STAR,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> stencil_width=1,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # boundary_type=PETSc.DMDA.BoundaryType.GHOSTED,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> )</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> x = da.createGlobalVec()</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> x_old = da.createGlobalVec()</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f = da.createGlobalVec()</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> J = da.createMat()</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> rho_ref = rho_m[0] # kg/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> e_ref = e_m[0] # J/mol</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> p_ref = p0 # Pa</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> x.setArray(np.array([rho_m / rho_ref, e_m / e_ref, ux_m]).T.flatten())</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> x_old.setArray(np.array([rho_m / rho_ref, e_m / e_ref, ux_m]).T.flatten())</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB = PETSc.Options()</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["snes_lag_preconditioner_persists"] = False</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["snes_lag_jacobian"] = 1</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["snes_lag_jacobian_persists"] = False</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["snes_lag_preconditioner"] = 1</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black">optsDB["ksp_type"] = "gmres" # "gmres" # gmres"</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="NO-BOK" style="font-size:11.0pt;color:black"> </span><span style="font-size:11.0pt;color:black">optsDB["pc_type"] = "ilu" # "lu" # "ilu"</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["snes_type"] = "newtonls"</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["ksp_rtol"] = 1e-7</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["ksp_atol"] = 1e-7</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["ksp_max_it"] = 100</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["snes_rtol"] = 1e-5</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["snes_atol"] = 1e-5</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["snes_stol"] = 1e-5</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["snes_max_it"] = 100</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["snes_mf"] = False</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["ts_max_time"] = t_end</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black">optsDB["ts_type"] = "beuler" # "bdf" #</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="NO-BOK" style="font-size:11.0pt;color:black"> </span><span style="font-size:11.0pt;color:black">optsDB["ts_max_snes_failures"] = -1</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["ts_monitor"] = ""</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["ts_adapt_monitor"] = ""</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # optsDB["snes_monitor"] = ""</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # optsDB["ksp_monitor"] = ""</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> optsDB["ts_atol"] = 1e-4</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> x0 = x_old</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> residual_wrap = residual_ts(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> eos,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> x0,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> N_vertical,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> g,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> pos,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> z,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> mw,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> dt,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> dx,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> p_amb,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> A_nozzle,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black">r_tank_inner,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="NO-BOK" style="font-size:11.0pt;color:black"> mph_uv_flsh_L,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="NO-BOK" style="font-size:11.0pt;color:black"> </span><span style="font-size:11.0pt;color:black">rho_ref,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> e_ref,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> p_ref,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> closed_tank,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> J,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> da,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> drift_func,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> T_wall,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> tank_params,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> )</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # optsDB["ts_adapt_type"] = "none"</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ts = PETSc.TS().create(comm=COMM)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # TODO: Figure out why DM crashes the code</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # ts.setDM(residual_wrap.da)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ts.setIFunction(residual_wrap.residual_ts, None)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ts.setTimeStep(dt)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ts.setMaxSteps(-1)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ts.setTime(t_start) # s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ts.setMaxTime(t_end) # s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ts.setMaxSteps(1e5)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ts.setStepLimits(1e-3, 1e5)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ts.setFromOptions()</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ts.solve(u=x)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"><br>
<br>
Residual function:<br>
<br>
class residual_ts:</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> def __init__(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> eos,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> x0,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> N,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> g,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> pos,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> z,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> mw,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> dt,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> dx,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> p_amb,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> A_nozzle,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black">r_tank_inner,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="NO-BOK" style="font-size:11.0pt;color:black"> mph_uv_flsh_l,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="NO-BOK" style="font-size:11.0pt;color:black"> </span><span style="font-size:11.0pt;color:black">rho_ref,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> e_ref,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> p_ref,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> closed_tank,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> J,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> da,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> drift_func,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> T_wall,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> tank_params,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ):</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.eos = eos</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.x0 = x0</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.N = N</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.g = g</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.pos = pos</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.z = z</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-_73AJ_I$" target="_blank">self.mw</a> = mw</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black">self.dt = dt</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="NO-BOK" style="font-size:11.0pt;color:black"> self.dx = dx</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="NO-BOK" style="font-size:11.0pt;color:black"> </span><span style="font-size:11.0pt;color:black">self.p_amb = p_amb</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.A_nozzle = A_nozzle</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black">self.r_tank_inner = r_tank_inner</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="NO-BOK" style="font-size:11.0pt;color:black"> self.mph_uv_flsh_L = mph_uv_flsh_l</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="NO-BOK" style="font-size:11.0pt;color:black"> </span><span style="font-size:11.0pt;color:black">self.rho_ref = rho_ref</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.e_ref = e_ref</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.p_ref = p_ref</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.closed_tank = closed_tank</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span lang="NO-BOK" style="font-size:11.0pt;color:black">self.J = J</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="NO-BOK" style="font-size:11.0pt;color:black"> self.f = f</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="NO-BOK" style="font-size:11.0pt;color:black"> self.da = da</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span lang="NO-BOK" style="font-size:11.0pt;color:black"> </span><span style="font-size:11.0pt;color:black">self.drift_func = drift_func</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.T_wall = T_wall</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.tank_params = tank_params</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.Q_wall = np.zeros(N)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.n_iter = 0</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.t_current = [0.0]</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.s_top = [0.0]</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.p_choke = [0.0]</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # setting interp func # TODO figure out how to generalize this method</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self._interp_func = _jalla_upwind</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # allocate space for new params</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.p = np.zeros(N) # Pa</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.T = np.zeros(N) # K</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.alpha = np.zeros((2, N))</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho = np.zeros((2, N))</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.e = np.zeros((2, N))</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # allocate space for ghost cells</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.alpha_ghost = np.zeros((2, N + 2))</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_ghost = np.zeros((2, N + 2))</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_m_ghost = np.zeros(N + 2)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.u_m_ghost = np.zeros(N + 1)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.u_ghost = np.zeros((2, N + 1))</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.e_ghost = np.zeros((2, N + 2))</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.pos_ghost = np.zeros(N + 2)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.h_ghost = np.zeros((2, N + 2))</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # allocate soace for local X and Xdot</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.X_LOCAL = da.createLocalVec()</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.XDOT_LOCAL = da.createLocalVec()</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> def residual_ts(self, ts, t, X, XDOT, F):</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.n_iter += 1</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # TODO: Estimate time use</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> """</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> Caculate residuals for equations</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> (rho, rho (e + gx))_t + (rho u, rho u (h + gx))_x = 0</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> P_x = - g \rho</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> """</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> n_phase = 2</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.da.globalToLocal(X, self.X_LOCAL)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.da.globalToLocal(XDOT, self.XDOT_LOCAL)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> x = self.da.getVecArray(self.X_LOCAL)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> xdot = self.da.getVecArray(self.XDOT_LOCAL)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f = self.da.getVecArray(F)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> T_c, v_c, p_c = self.eos.critical(self.z) # K, m3/mol, Pa</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> rho_m = x[:, 0] * self.rho_ref # kg/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> e_m = x[:, 1] * self.e_ref # J/mol</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> u_m = x[:-1, 2] # m/s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # derivatives</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> rho_m_dot = xdot[:, 0] * self.rho_ref # kg/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> e_m_dot = xdot[:, 1] * self.e_ref # kg/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> dt = ts.getTimeStep() # s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> for i in range(self.N):</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # get new parameters</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.mph_uv_flsh_L[i] = self.eos.multi_phase_uvflash(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.z, e_m[i], <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-_73AJ_I$" target="_blank">self.mw</a> / rho_m[i], self.mph_uv_flsh_L[i]</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> )</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> betaL, betaV, betaS = _get_betas(self.mph_uv_flsh_L[i]) # mol/mol</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> beta = [betaL, betaV]</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> if betaS != 0.0:</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> print("there is a solid phase which is not accounted for")</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.T[i], self.p[i] = _get_tank_temperature_pressure(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.mph_uv_flsh_L[i]</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ) # K, Pa)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> for j, phase in enumerate([self.eos.LIQPH, self.eos.VAPPH]):</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # new parameters</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_ghost[:, 1:-1][j][i] = (</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-_73AJ_I$" target="_blank">self.mw</a></span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> / self.eos.specific_volume(self.T[i], self.p[i], self.z, phase)[0]</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ) # kg/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.e_ghost[:, 1:-1][j][i] = self.eos.internal_energy_tv(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.T[i], <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-_73AJ_I$" target="_blank">self.mw</a> / self.rho_ghost[:, 1:-1][j][i], self.z, phase</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> )[</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> 0</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ] # J/mol</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.h_ghost[:, 1:-1][j][i] = (</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.e_ghost[:, 1:-1][j][i]</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> + self.p[i] * <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-_73AJ_I$" target="_blank">self.mw</a> / self.rho_ghost[:, 1:-1][j][i]</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ) # J/mol</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.alpha_ghost[:, 1:-1][j][i] = (</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> beta[j] / self.rho_ghost[:, 1:-1][j][i] * rho_m[i]</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ) # m3/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # calculate drift velocity</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> for i in range(self.N - 1):</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.u_ghost[:, 1:-1][0][i], self.u_ghost[:, 1:-1][1][i] = (</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> calc_drift_velocity(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> u_m[i],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self._interp_func(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_ghost[:, 1:-1][0][i],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_ghost[:, 1:-1][0][i + 1],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> u_m[i],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ),</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self._interp_func(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_ghost[:, 1:-1][1][i],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_ghost[:, 1:-1][1][i + 1],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> u_m[i],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ),</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.g,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self._interp_func(self.T[i], self.T[i + 1], u_m[i]),</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> T_c,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.r_tank_inner,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self._interp_func(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.alpha_ghost[:, 1:-1][0][i],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.alpha_ghost[:, 1:-1][0][i + 1],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> u_m[i],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ),</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self._interp_func(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.alpha_ghost[:, 1:-1][1][i],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.alpha_ghost[:, 1:-1][1][i + 1],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> u_m[i],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ),</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.drift_func,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> )</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ) # liq m / s , vapour m / s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> u_bottom = 0</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> if self.closed_tank:</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> u_top = 0.0 # m/s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> else:</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # calc phase to skip env_isentrope_cross</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> if (</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.mph_uv_flsh_L[-1].liquid != None</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> and self.mph_uv_flsh_L[-1].vapour == None</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> and self.mph_uv_flsh_L[-1].solid == None</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ):</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> phase_env = self.eos.LIQPH</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> else:</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> phase_env = self.eos.TWOPH</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.h_m = e_m + self.p * <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-_73AJ_I$" target="_blank">self.mw</a> / rho_m # J / mol</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> mdot, self.p_choke[0] = calc_mass_outflow(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.eos,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.z,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.h_m[-1],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.s_top[0],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.p[-1],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.p_amb,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.A_nozzle,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-_73AJ_I$" target="_blank">self.mw</a>,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> phase_env,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> debug_plot=False,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ) # mol / s , Pa</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> u_top = -mdot * <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-_73AJ_I$" target="_blank">self.mw</a> / rho_m[-1] / (np.pi * self.r_tank_inner**2) # m/s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # assemble vectors with ghost cells</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.alpha_ghost[:, 0] = self.alpha_ghost[:, 1:-1][:, 0] # m3/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.alpha_ghost[:, -1] = self.alpha_ghost[:, 1:-1][:, -1] # m3/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_ghost[:, 0] = self.rho_ghost[:, 1:-1][:, 0] # kg/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_ghost[:, -1] = self.rho_ghost[:, 1:-1][:, -1] # kg/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_m_ghost[0] = rho_m[0] # kg/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_m_ghost[1:-1] = rho_m # kg/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_m_ghost[-1] = rho_m[-1] # kg/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # u_ghost[:, 1:-1] = u # m/s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.u_ghost[:, 0] = u_bottom # m/s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.u_ghost[:, -1] = u_top # m/s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.u_m_ghost[0] = u_bottom # m/s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.u_m_ghost[1:-1] = u_m # m/s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.u_m_ghost[-1] = u_top # m/s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.e_ghost[:, 0] = self.e_ghost[:, 1:-1][:, 0] # J/mol</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.e_ghost[:, -1] = self.e_ghost[:, 1:-1][:, -1] # J/mol</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.pos_ghost[1:-1] = self.pos # m</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.pos_ghost[0] = self.pos[0] # m</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.pos_ghost[-1] = self.pos[-1] # m</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.h_ghost[:, 0] = self.h_ghost[:, 1] # J/mol</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.h_ghost[:, -1] = self.h_ghost[:, -2] # J/mol</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # recalculate wall temperature and heat flux</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # TODO ARE WE DOING THE STAGGERING CORRECTLY?</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> lz = self.tank_params["lz_tank"] / self.N # m</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> if ts.getTime() != self.t_current[0] and self.tank_params["heat_transfer"]:</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.t_current[0] = ts.getTime()</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> for i in range(self.N):</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.T_wall[i], self.Q_wall[i], h_ht = (</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> solve_radial_heat_conduction_implicit(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.tank_params,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.T[i],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.T_wall[i],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> (self.u_m_ghost[i] + self.u_m_ghost[i + 1]) / 2,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_m_ghost[i + 1],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.mph_uv_flsh_L[i],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> lz,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> dt,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> )</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ) # K, J/s, W/m2K</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # Calculate residuals</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f[:, :] = 0.0</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f[:, 0] = dt * rho_m_dot # kg/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;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:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f[:, 2] = (</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> dt</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> * (</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> rho_m_dot * (e_m / <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-_73AJ_I$" target="_blank">self.mw</a> + self.g * self.pos)</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> + rho_m * e_m_dot / <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-_73AJ_I$" target="_blank">self.mw</a></span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> )</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> - rho_m_dot * e_m_dot / <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-_73AJ_I$" target="_blank">self.mw</a> * dt**2</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> - self.Q_wall / (np.pi * self.r_tank_inner**2 * lz) * dt</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> ) # J / m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # add contribution from space</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> for i in range(n_phase):</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> e_flux_i = np.zeros_like(self.u_ghost[i]) # J/m3 m/s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> rho_flux_i = np.zeros_like(self.u_ghost[i]) # kg/m2/s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> for j in range(1, self.N + 1):</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> if self.u_ghost[i][j] >= 0.0:</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> rho_flux_new = _rho_flux(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.alpha_ghost[i][j], self.rho_ghost[i][j], self.u_ghost[i][j]</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> )</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> e_flux_new = _e_flux(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.alpha_ghost[i][j],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_ghost[i][j],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.h_ghost[i][j],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-_73AJ_I$" target="_blank">self.mw</a>,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.g,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.pos_ghost[j],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.u_ghost[i][j],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> )</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # backward euler</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> rho_flux_i[j] = rho_flux_new # kg/m2/s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> e_flux_i[j] = e_flux_new # J/m3 m/s</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> else:</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> rho_flux_new = _rho_flux(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.alpha_ghost[i][j + 1],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_ghost[i][j + 1],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.u_ghost[i][j],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> )</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> e_flux_new = _e_flux(</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.alpha_ghost[i][j + 1],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.rho_ghost[i][j + 1],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.h_ghost[i][j + 1],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> <a href="https://urldefense.us/v3/__http://self.mw/__;!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-_73AJ_I$" target="_blank">self.mw</a>,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.g,</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.pos_ghost[j + 1],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> self.u_ghost[i][j],</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> )</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # backward euler</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> rho_flux_i[j] = rho_flux_new</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> e_flux_i[j] = e_flux_new</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # mass eq</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f[:, 0] += (dt / self.dx) * (rho_flux_i[1:] - rho_flux_i[:-1]) # kg/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # energy eq</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f[:, 2] += (dt / self.dx) * (e_flux_i[1:] - e_flux_i[:-1]) # J/m3</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f1_ref, f2_ref, f3_ref = self.rho_ref, self.p_ref, self.e_ref</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f[:, 0] /= f1_ref</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f[:-1, 1] /= f2_ref</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f[:, 2] /= f3_ref</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> # dummy eq</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> f[-1, 1] = x[-1, 2]</span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt;color:black"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:11.0pt"> </span><span style="font-size:12.0pt"><o:p></o:p></span></p>
</div>
</div>
</div>
</blockquote>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"><br clear="all">
<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span class="m-1124721174010954366m-6386679822550478849m3355120538688373402gmailsignatureprefix"><span style="font-size:12.0pt">--
</span></span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt">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<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"><a href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-EmFBV7E$" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><o:p></o:p></span></p>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"><br clear="all">
<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span class="m-1124721174010954366m-6386679822550478849gmailsignatureprefix"><span style="font-size:12.0pt">--
</span></span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt">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<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"><a href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-EmFBV7E$" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><o:p></o:p></span></p>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"><br clear="all">
<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span class="m-1124721174010954366gmailsignatureprefix"><span style="font-size:12.0pt">--
</span></span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt">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<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"> <o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal" style="mso-margin-top-alt:auto;mso-margin-bottom-alt:auto"><span style="font-size:12.0pt"><a href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-EmFBV7E$" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><o:p></o:p></span></p>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"><br clear="all">
<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
</div>
<p class="MsoNormal"><span class="gmailsignatureprefix"><span style="font-size:12.0pt">--
</span></span><span style="font-size:12.0pt"><o:p></o:p></span></p>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">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<o:p></o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"><a href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!YLsa1SRpJkLsV8gFq3kXXxcNwmZ3lx9qJMebtTjTnxhU0-43s0nvM3UbT2FOmlqyxa9bU2er_jIx3OYk_Ds4kW542MX-EmFBV7E$" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><o:p></o:p></span></p>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</body>
</html>