Auxiliary fields in PeleLMeX with MAGMA

[q86091483]

Dear community,

I'm implementing the transport equation of auxiliary variables, such as residence times (or age variable) and post-flame residence times (hereafter referred to aux ) in PeleLMeX. These variables can be considered as passive scalars and do not affect the evolution of variables state[0:NVAR]originally transported in PeleLMeX.

Our implementation of the residence times worked pretty well on CPU. Specifically, a second CVODE solver was initialized and used to advance the aux variables, with the corresponding cF_RHS_aux and cJac_aux added. The interface and implementation of utils::flattenand ùtils::unflatten are also changed accordingly. The reason we added a second CVODE, instead of using the the original CVODE for the prolonged ODE vector [$\rho Y_K$, $T$, $aux$], is that we found the addition of aux variables may change the original solution, possibly due to the addition of Jacobian of aux fields (here we use denseAJ_direct with ReactorCvode with analytical Jacobian).
Below is the transport of mixture fraction of the cross flow, residence time and post-flame residence time in a H2 JICF. These variables are added following the spirit of SDC. We also did some tests in other manufactured cases, the results agree well with some available analytical solutions.

Our problem is the GPU. We now face some issues when implementing the same equation with MAGMA on AMD HIP GPU (we need to burn through some GPU resources). An error occurs when using SUNLinSol_MagmaDense to set the linear solver of CVODE named as LS_aux. Here is how we did on GPU:

Similar to what we have done for CPU, a second initCvode_aux, cvode::cJac_aux, cF_RHS_aux and fKernelAJCHem are added for GPU. On the other hand in the GPU region of ReactorCvode::react(...), we first declare

SUNMatrix A_aux = nullptr;            // Jacobian matrix for aux
SUNNonlinearSolver NLS_aux = nullptr; // Non-linear solver for aux
SUNLinearSolver LS_aux = nullptr;     // Linear solver for aux

void* cvode_mem_aux = CVodeCreate(CV_BDF, *amrex::sundials::The_Sundials_Context());

The ODE vector for aux is then declared as

const int neq_tot_aux = NUMAUX * ncells; // say NUMAUX = 2 for a mixture fraction and mixture fraction-weighted residence time
auto y_aux = utils::setNVectorGPU(neq_tot_aux, atomic_reductions, stream);
amrex::Real* yvec_d_aux = N_VGetDeviceArrayPointer(y_aux);

After flatten(box, ncells...), we initialized the CVODE components as follows.

int flag1 = CVodeSetUserData(cvode_mem_aux, udata);
flag1 = CVodeInit(cvode_mem_aux, cF_RHS_aux, time_start, y_aux);
if (utils::check_flag(&flag1, "CVodeInit", 1)) {
  return (1);
}

// Initializey_aux, A_aux
LS_aux = SUNLinSol_MagmaDense(y_aux, A_aux, *amrex::sundials::The_Sundials_Context()); // The problem is here.
if (utils::check_flag(
  static_cast<void*>(LS_aux), "SUNLinSol_MagmaDense", 0)) {
  return (1);
}

The error occurs as
SUNDIALS_ERROR: SUNLinSol_MagmaDense() failed - returned NULL pointer
which means SUNLinSol_MagmaDense return a null pointer and didn't initliaze y_aux properly.

I want to ask:

1. Can you give me some hints on how to solve the issue caused by SUNLinSol_MagmaDense?
2. I suspect I didn't deal with the GPU streams properly, since streamis used in utils::setNVectorGPU. I don't understand how the stream works in amrex and Pele*. Sorry I looked into the online amrex doc and the past publication but I still didn't grasp the gist. Also The sundials document gives very limited introduction on N_VSetKernelExecPolicy_Hip related topics. Can you help to give me some reading direction, like how amrex::gpuStream is used in PeleLMeX?
3. What role is played by stream in ReactorCvode::react(..., stream)?
4. Is it okay to use the same streamfor both CVODE solvers (original and auxiliaray fields)?

Any comments are greatly welcomed.

thanks
Zisen

[q86091483]

Just add: the above code snippet occurs in ReactorCvode::react(...).

[baperry2]

I'm not the best person to give advice on the specific questions being asked, but I have a couple comments on overall strategy. I'd expect the aux quantities you are looking at to be non-stiff, and also with them being passive, stiff integration with CVODE using an analytical Jacobian may be overkill.

Another option would be to treat the source terms as external sources as we do with source terms from the soot and spray modules. See #427 for examples of this with non-transported auxiliary quantities. This strategy is decoupled from the SDC iterations (#429) and essentially amounts to forward Euler integration of the sources, so it definitely isn't applicable for everything, but probably would be fine for mixture fraction (no source) and residence time (constant source term).

Also keep any eye out for #467 which will provide an improved framework for supporting auxiliary scalars.

[q86091483]

Thanks Bruce.

The passive scalar we want to add does introduce some stiffness, and we found that implicit CVODE can help to solve this issue. I will also keep an eye on the threads you recommended.