[FLASH-USERS] MatRad3 case for MultiTemp

[Klaus Weide]

On Wed, 11 Nov 2015, Mark Richardson wrote:

>    I'm curious if someone with experience with the MultiTemp aspect of
> FLASH's EOS module, in particular (but not necessary) the MatRad3 case
> which sets the ion component to zero, could help me understand it better.

Hello Mark,

The Users Guide does need some updating regarding the Eos unit, in 
particular the MatRad3 variants of the multiTemp implementations.
In the meantime I'll try to answer your questions as much as possible.

> I am trying to implement a new EOS that works within MatRad3. I have this EOS
> working with Multispecies, but for the MatRad3 scenario I'm a bit confused
> about when the NR solver needs to be called (so Table 15.1 in the user
> guide).
>
>   Why does the Newton-Raphson solver need to be called in the MODE_DENS_EI
> cases in Gamma and Multigamma, given that eint is solved analytically in
> the case of MODE_DENS_TEMP. Why isn't temp set just with the inverse? 

First, you probably should not be using the MODE_DENS_EI or MODE_DENS_TEMP
variants in a 3T (including MatRad3) scenario. You should be using 
MODE_DENS_EI_GATHER (or similar) and MODE_DENS_TEMP_<something> instead.
While the 1T modes (MODE_DENS_EI,MODE_DENS_EI,MODE_DENS_PRES) have been 
left in the code of the 3T implementations, for some possible use in the 
future (according to table 15.1, when one really does not care at all for 
the component-specific output variables), those code paths have not been 
recently used or tested by us.

Let me try to answer the question "Why does the Newton-Raphson solver need 
to be called in the MODE_DENS_EI_GATHER case?" (which is what you may have 
meant). 

Calling Eos with a *_GATHER mode variant is taken to imply that you want 
to update the component-specific output variables (like 
Tion,Tele,Trad,pion,pele,prad, in the case of MODE_DENS_EI_GATHER),
but in addition also are interested in having them "gathered" into
"combined" output variables (like T,p,ei). For pressures and energies,
p and ei are just sums over the three components. For a combined 
temperature, T is understood to be the temperature that the three 
components in a cell WOULD equilibrate to IF you let them.

To compute this as-if temperature T, the various multiTemp Eos 
implementations use a Newton-Raphson(-like) root search (RS).  When 
getting this combined T is the only reason for doing the RS, this is 
indicated in Table 15.1 with a parenthesized "(y)" in the "RS" column.
Because the output variable T is rarely needed, several multiTemp Eos 
implementations have been modified to check a runtime parameter in order 
to determine whether to actually compute it. See code like
   	if(eos_combinedTempRule==0) then
           call eos_newtonRaphson(...)
        endif
in various places. (This is not reflected in the UG test around Table 15.1 
yet.)

Why use a Newton-Raphson RS instead of just inverting the 
energy(temperature) dependence analytically? This would be trivial for a 
purely ideal-gas-like temperature dependence of energy, but that is not 
what we have in the presence of radiation. In general we have a 
temperature dependence of the form
   ei(T) = a*T + b*T^4 .
Finding T for a given ei(T) amounts to solving a quartic equation, 
which COULD be done analytically but we have chosen to use the 
pre-existing RS code instead.


> My
> confusion makes understanding the more complex MODE_DENS_EI_{ } submodes
> quite difficult.
> 
> Also, do I have much control beyond Simulation_initBlock in how EOS will be
> called throughout the simulation, where we are including RadTrans and
> Diffusion. It seems like for the most part it's limited to MODE_DENS_EI,
> but there are some other cases, including _EI_MAT_GATHER_PRAD_SCALE
> (hy_uhd_unsplit), and _PRES (hy_uhd_getFaceFlux).

There are now four runtime parameters that specify Eos modes to be used in 
various situations by the FLASH code:

eosMode     - Generic mode used by the Grid unit when it is asked to make
              sure that cells are thermodynamically consistent,
              and not given more specific instructions
              (as possibly after guard cell filling, and after block
              creation or destruction [refinement or derefinement])
eosModeInit - as eosMode but used during initialization.
              During initial construction of the Grid, each call of
              Simulation_initBlock will be automatically followed by
              a call of Eos_wrapped(eosModeInit,...) for that block.
hy_eosModeGc  - Used by unsplit Hydro (including USM MHD) to update
                guard cells whose variable values may be the result
                of interpolation or averaging at fine-coarse boundaries,
                before those values are used by the Hydro code as inputs
                for reconstruction etc.
                (Default setting is to use the mode given by eosMode.)
hy_eosModeAfter - Used by unsplit Hydro when it calls Eos_wrapped to 
                finish up its work on a block.

In addition, there are various modes hardwired in various code units...

The _PRADSCALE variants are (currently) used only in connection with 
the radiation flux-limiter aware Hydro implementation (under 
source/physics/Hydro/HydroMain/unsplit_rad). They are not implemented at 
the Eos() level or lower, but are caught by Eos_wrapped() [which calls 
Eos() with a "real" mode and then modifies some of the output variables to 
effectively "scale down" the radiation pressure using the radiation flux 
limiter].

The calls with MODE_DENS_PRES in hy_uhd_getFaceFlux are only for use of 
Hydro with EOSforRiemann=TRUE. I believe this has not been used or tested 
in 3T setups. To do that, one would have to change hy_uhd_getFaceFlux
to call Eos with a different mode in the 3T case (and provide the 
appropriate input variables for that mode).



HTH,

Klaus