<html xmlns:v="urn:schemas-microsoft-com:vml" 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=us-ascii">
<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:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
font-size:12.0pt;
font-family:"Calibri",sans-serif;
mso-ligatures:standardcontextual;}
p.MsoListParagraph, li.MsoListParagraph, div.MsoListParagraph
{mso-style-priority:34;
margin-top:0in;
margin-right:0in;
margin-bottom:0in;
margin-left:.5in;
font-size:12.0pt;
font-family:"Calibri",sans-serif;
mso-ligatures:standardcontextual;}
span.EmailStyle17
{mso-style-type:personal-compose;
font-family:"Calibri",sans-serif;
color:windowtext;}
.MsoChpDefault
{mso-style-type:export-only;}
@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:1870145129;
mso-list-type:hybrid;
mso-list-template-ids:-118971206 67698703 67698713 67698715 67698703 67698713 67698715 67698703 67698713 67698715;}
@list l0:level1
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level2
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level3
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
@list l0:level4
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level5
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level6
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
@list l0:level7
{mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level8
{mso-level-number-format:alpha-lower;
mso-level-tab-stop:none;
mso-level-number-position:left;
text-indent:-.25in;}
@list l0:level9
{mso-level-number-format:roman-lower;
mso-level-tab-stop:none;
mso-level-number-position:right;
text-indent:-9.0pt;}
ol
{margin-bottom:0in;}
ul
{margin-bottom:0in;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang="EN-US" link="#0563C1" vlink="#954F72" style="word-wrap:break-word">
<div class="WordSection1">
<p class="MsoNormal"><span style="font-size:11.0pt">Hello all,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">I am running a calculation using the Heatexchange module and have encountered a case where dt_HeatXc becomes extremely small in cold cells.
</span> <span style="font-size:11.0pt">I am aware of hx_relTol and hx_dtFactor; hx_relTol can relax this behavior and hx_dtFactor scales the resulting HeatExchange timestep, but neither seems to directly address the underlying behavior I am seeing. I have
made a local modification and wanted to ask whether I have overlooked something, or whether a slightly different timestep estimate might be useful more generally.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">As I understand it, the timestep estimate effectively linearly extrapolates the initial exchange rate over the proposed timestep. This appears unnecessarily restrictive when tauEq is much shorter than the
timestep, because the temperature difference itself decays exponentially. Once the electron and ion temperatures equilibrate, the driving term disappears, so Qdot*dt substantially overestimates the energy that can actually be exchanged.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">Specifically, dt_HeatXc is the timestep constraint supplied by the HeatExchange unit and is compared with timestep constraints from the other units in determining the global timestep. It is of the form D*Eint/Qdot,
where D is the max of either dT/T or hx_relTol, and dT = abs(Tion – Tele). However, Qdot has the form CvEle*dT/tauEq, where tauEq is the electron-ion equilibration time. Thus, when dT/T is greater than hx_relTol, this gives dt_heatXc = tauEq*Eint/(CvEle*T).
In this case, dt_heatXc is proportional to the equilibration time</span> <span style="font-size:11.0pt">
and the explicit dependence on the electron-ion temperature difference cancels. This allows a cell that is arbitrarily close to equilibrium to still drive the scale of the global timestep.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">The relevant portion of the current code, omitting checks, second component, etc., is, from Heatexchange_computeDt.F90, approximately,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">q2dot = cvele / eqtime * t12diff<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">q1dot = -q2dot<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">…<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">denom = max(abs(t12diff)*eint1/temp1, hx_relTol*eint1)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">energy1LogRate = abs(q1dot) / denom<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">…<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">dt_temp = 1.0 / energyLogRate<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">dt_HeatXc = dt_temp*hx_dtFactor<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">I find this method to be problematic for ambient, near-equilibrium cells. For example, I encounter cells ahead of a shock with diagnostics,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">blockID = 101<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">temp1 = 297.731<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">temp2 = 296.865<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">t12diff = 8.659E-01<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">eqtime = 6.257E-17<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">energyLogRate = 9.323E+14<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">dt_HeatXc = 1.073E-15<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">Here the temperatures differ by only about 0.3%, but the small equilibration time produces a large instantaneous exchange rate and consequently a very small timestep.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">This appears to be the regime in which the linear extrapolation becomes overly restrictive. The actual Heatexchange update relaxes the two temperatures exponentially, so the exchange rate goes to zero as equilibrium
is approached and does not remain equal to its initial value over the entire timestep.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoListParagraph" style="margin-left:0in"><span style="font-size:11.0pt">As an immediate and minimal solution I considered excluding cells that are relatively cold and close to equilibrium from participating in computation of dt_HeatXc. In Heatexchange_computeDt.F90,
something like:<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> absTempThreshold = 1000<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> relTempThreshold = 0.01<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> if( doTempFilter ) then<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> relTDiff = abs(t12diff) / max(temp1, temp2)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> if( max(temp1, temp2) < absTempThreshold .AND. relTDiff < relTempThreshold) then<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> ! This cell will not restrict dt<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> cycle<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> end if<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> end if<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">This ameliorates the immediate problem, but I think there may be a more general way to formulate the timestep criterion that does not require an arbitrary cold-cell threshold.
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">If I interpret the existing denom as the energy-change scale that the timestep estimator is intended to resolve, it seems possible to retain that criterion but calculate the time required to exchange that
energy using the same exponential relaxation employed by Heatexchange.F90, rather than the initial-rate approximation. For the existing relaxation scheme, with CvIon and CvEle held constant during the step, the maximum energy that can possibly be exchanged
before the two components reach equilibrium is<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">Qmax = (CvEle*CvIon)/(CvEle + CvIon)*abs( dT )<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">This can then be compared directly to the existing energy-change scale used by the timestep estimator. Keeping the current limit of dT/T*Eint (when above the hx_relTol threshold),<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">Qallowed = dT/T*Eint<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">If Qmax < Qallowed, then even complete equilibration won’t exceed the allowed energy change, so this cell doesn’t need to restrict timestep.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">Otherwise, if Qmax exceeds Qallowed, calculate the timestep that would result in Qallowed heat exchange:<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">dtAllowed = -tauEq/(1 + CvEle/CvIon)*ln(1 – Qallow/Qmax)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">Something schematically like (missing checks if dT is zero, energies positive, etc):<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> if( doEnergyFilter ) then<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> dT = abs( t12diff )<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Qmax = (cvele*cvion)/(cvele + cvion)*dT
<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> <o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Qallow1 = max( dT*eint1/temp1, hx_relTol*eint1)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Qallow2 = max( dT*eint2/temp2, hx_relTol*eint2)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> <o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> Qallow = min(Qallow1, Qallow2)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> <o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> if(Qmax > Qallow) then<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> eqtime = eqtime * hx_ieTimeCoef<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> dt_temp_cell = -( eqtime/( 1.0 + cvele/cvion )) * LOG(1.0 - Qallow/Qmax)<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> if( dt_temp_cell < dt_temp ) then<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> dt_temp = dt_temp_cell<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> end if<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> end if<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> end if <o:p>
</o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"> <o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">Does this interpretation of the purpose of denom and the timestep calculation look correct? In particular, is the intent of dt_HeatXc to bound an energy change of this scale, or is it deliberately intended
to resolve the electron-ion relaxation timescale itself even when the finite energy available to exchange is small? If the former, is there a reason not to use the finite exponential exchange rather than its initial-rate approximation in the stiff limit?<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">It is worth noting that in the limit Qallow << Qmax, the proposed methodology reduces to the current initial-rate-extrapolation timestep estimate.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt">Thank you!<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:11.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal">Julian<o:p></o:p></p>
</div>
</body>
</html>