Set an ignition source with specific pressure and temperature in a certain area

[XUFEI-Sai]

Hello,
I want to set an ignition source with specific pressure and temperature in a certain area, such as prob.Hpoint_p=0.5e7, prob.Hpoint-T=2500. And I referred to the code in External.cpp and the previous #384 discussion "Add a heat source with a specified heat flux", trying to define problem_modify_ext_sources in my prob.H, but running the program always results in segment errors (core has been dumped). Here is my code, I'm not sure if this is appropriate. I hope to receive advice.
Thanks a lot.

`struct MyProblemSpecificFunctions : public DefaultProblemSpecificFunctions
{
static void problem_modify_ext_sources(
amrex::Real time,
amrex::Real /dt/,
const amrex::MultiFab& /state_old/,
const amrex::MultiFab& state_new,
amrex::MultiFab& ext_src,
int ng,
amrex::GeometryData const& geomdata,
ProbParmDevice const& prob_parm)
{
/* Notes: ext_src contains sources from velocity forcing coming in
This function should add to rather than overwrite ext_src.*/

auto const& Sns = state_new.const_arrays();
const amrex::Real* prob_lo = geomdata.ProbLo();
const amrex::Real* dx = geomdata.CellSize();
const amrex::IntVect ngs(ng);
auto const& Farrs = ext_src.arrays();

amrex::ParallelFor(
  ext_src, ngs, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
    const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0];
    const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1];
    const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];
    if (x <= prob_parm.Hpoint_x+prob_parm.Hpoint_r && x >= prob_parm.Hpoint_x
    && y <= prob_parm.Hpoint_r && y >= 0.0
    && time <= prob_parm.Hpoint_time_e 
    && time >= prob_parm.Hpoint_time_s) {
      
      amrex::Real rho_hpoint = 0.0;
      amrex::Real e_force = 0.0;
      amrex::Real e_hpoint = 0;
      amrex::Real massfrac_hpoint[NUM_SPECIES] = {0.0};
      for (int n = 0; n < NUM_SPECIES; n++){
        massfrac_hpoint[n] = Sns[nbx](i, j, k, UFS + n) / Sns[nbx](i, j, k, URHO);
      }
      auto eos = pele::physics::PhysicsType::eos();
      eos.PYT2RE(prob_parm.Hpoint_p,massfrac_hpoint,prob_parm.Hpoint_T,rho_hpoint,e_hpoint);
      
      for (int idir = 0; idir < AMREX_SPACEDIM; idir++) {
        
        Farrs[nbx](i, j, k, UMX + idir) = Sns[nbx](i, j, k, UMX + idir) / Sns[nbx](i, j, k, URHO) *
                                          rho_hpoint;
        e_force += Sns[nbx](i, j, k, UMX + idir) / Sns[nbx](i, j, k, URHO) *
                   Farrs[nbx](i, j, k, UMX + idir);
      }
      Farrs[nbx](i, j, k, URHO)  = rho_hpoint;
      Farrs[nbx](i, j, k, UEDEN) = rho_hpoint * e_hpoint + 0.5 * e_force;
      Farrs[nbx](i, j, k, UEINT) = rho_hpoint * e_hpoint;
      Farrs[nbx](i, j, k, UTEMP) = prob_parm.Hpoint_T;
      
    }
  });

}
}`

[marchdf]

Segfault should be easy to identify/solve. I recommend building in debug mode and then looking at the backtrace. It should point you to the exact line that fails. If you have trouble parsing it, feel free to post the trace here.

[XUFEI-Sai]

After trying, in my code, I found that in ParallelFor, it is not possible to use "if" to judge and modify values other than Farrs [nbx] (i, j, k, UEDEN); otherwise, segfault will occur. In addition, I tried calling external_forcing in the template to introduce a lot of speed when restarting the calculation, and the results did not seem to change much. May I ask if the problem_modify_ext_sources interface allows modifying other fluxes and specifying fluxes for a specific physical area? If there could be an example, it would be very helpful.

[baperry2]

As a point of clarification, problem_modify_ext_sources defines source terms acting on the cell volume, which get stored in Farrs and have dimensions of quantity/volume/time, not fluxes. It should be able to add a source term for any state variable. However, it appears that you are trying to use this capability to reset the temperature and pressure within a region of your simulation. That is not the purpose of this capability and will not work. You could instead compute a forcing that pushes the state toward your desired conditions. Or you could use different user-defined capability such as setting initial conditions with the pc_initdata function.

If you do want to add a source term, here is an example of a case that uses this capability to add a mass source term in a circular region: https://github.com/baperry2/PeleC/blob/54c8253eb53e887d45cf436d26dbb6d549d8de53/Exec/RegTests/EB-Vane/prob.H#L125. Computing the location and using a conditional should be able to work for defining a source in a particular area. Though note that if your case has EB, you should add an additional conditional to not compute the source term in covered cells. See the code linked above or the PeleC::fill_ext_source function for an example.

[XUFEI-Sai]

Thank you very much for your help. I have achieved it. I found that my problem is that the UEDEN of Farrs multiplied by the local density is on the same order of magnitude as the UEDEN of Sns, and only when there is a significant difference between the two, as in your EB_vane(where Sns [nbx] (i, j, k, UEDEN) is about 2e6, and Farrs [nbx] (i, j, k, UEDEN) is about 2e10), can it be effective.