[FLASH-USERS] Heatexchange_computeDt problems and possible solutions

[Julian Bouffard]

Hello all,

I am running a calculation using the Heatexchange module and have encountered a case where dt_HeatXc becomes extremely small in cold cells.  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.

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.

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 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.

The relevant portion of the current code, omitting checks, second component, etc., is, from Heatexchange_computeDt.F90, approximately,

q2dot = cvele / eqtime * t12diff
q1dot = -q2dot
...
denom = max(abs(t12diff)*eint1/temp1, hx_relTol*eint1)
energy1LogRate = abs(q1dot) / denom
...
dt_temp = 1.0 / energyLogRate
dt_HeatXc = dt_temp*hx_dtFactor

I find this method to be problematic for ambient, near-equilibrium cells.  For example, I encounter cells ahead of a shock with diagnostics,

blockID       = 101
temp1         = 297.731
temp2         = 296.865
t12diff       = 8.659E-01
eqtime        = 6.257E-17
energyLogRate = 9.323E+14
dt_HeatXc     = 1.073E-15

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.
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.


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:

           absTempThreshold = 1000
           relTempThreshold = 0.01
           if( doTempFilter ) then
              relTDiff = abs(t12diff) / max(temp1, temp2)
              if( max(temp1, temp2) < absTempThreshold .AND. relTDiff < relTempThreshold) then
                 ! This cell will not restrict dt
                 cycle
              end if
           end if

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.
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

Qmax = (CvEle*CvIon)/(CvEle + CvIon)*abs( dT )

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),

Qallowed = dT/T*Eint

If Qmax < Qallowed, then even complete equilibration won't exceed the allowed energy change, so this cell doesn't need to restrict timestep.
Otherwise, if Qmax exceeds Qallowed, calculate the timestep that would result in Qallowed heat exchange:

dtAllowed = -tauEq/(1 + CvEle/CvIon)*ln(1 - Qallow/Qmax)

Something schematically like (missing checks if dT is zero, energies positive, etc):

           if( doEnergyFilter ) then

              dT = abs( t12diff )
              Qmax = (cvele*cvion)/(cvele + cvion)*dT

              Qallow1 = max( dT*eint1/temp1, hx_relTol*eint1)
              Qallow2 = max( dT*eint2/temp2, hx_relTol*eint2)

              Qallow = min(Qallow1, Qallow2)

              if(Qmax > Qallow) then
                 eqtime = eqtime * hx_ieTimeCoef
                 dt_temp_cell = -( eqtime/( 1.0 + cvele/cvion )) * LOG(1.0 - Qallow/Qmax)

                 if( dt_temp_cell < dt_temp ) then
                    dt_temp = dt_temp_cell
                 end if
              end if
           end if

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?

It is worth noting that in the limit Qallow << Qmax, the proposed methodology reduces to the current initial-rate-extrapolation timestep estimate.

Thank you!

Julian
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://flash.rochester.edu/pipermail/flash-users/attachments/20260819/00aed8ba/attachment.htm>