Non-unity Lewis number

[bazharz]

Hello, Regarding the Lewis number. I would like to impose a constant number, for example, Le = 0.8. I tried to consult the documentation, but I couldn't find an explanation on this matter. Is it possible to set a non-unity Lewis number?

[baperry2]

By default, PeleLMeX runs with full mixture-averaged diffusion (i.e., variable non-unity Lewis numbers). There is in option for to enforce unity Lewis number for all species, but not an option to specify a constant non-unity Lewis number. You should be able to accomplish this with a fairly trivial change in the source code though, assuming you want the same Lewis number for all species.

To run with constant unity Lewis numbers, the option is peleLM.unity_Le = 1, in which case you also need to specify a Prandtl number peleLM.Prandtl = 0.7. In that case, thermal diffusivity is computed from viscosity using that Prandtl number, and mass diffusivities are computed from the viscosity using a Schmidt number equal to the Prandtl number as implied by Le=1. For a contant non-unity Lewis number across all species, all you need to do is modify the code here where the Schmidt number is set and add a separate query so it can be set independently of the Prandtl number.

BTW, if you are interested in contributing to the code, adding this capability would be a great opportunity for a first pull request!

[bazharz]

Thank you for this clarification! I will try to integrate this into this file!

[bazharz]

I'm reaching out to you again because I've attempted to fix a non-unity Lewis number.

I've modified 2 source files: PeleLMSetup.cpp and PeleLM_K.H. In the PeleLMSetup.cpp file, I set Sc = 1 and calculated Pr for Le = 0.7, Pr = Sc/Le. Then, in the PeleLM_K.H file, I checked if the quantity lambda/(rhoDCp) = 0.7, which is the imposed value of Le.

I'm not sure how the eos.TY2Cp function calculates Cp, and I couldn't find any information about it in the documentation. Could someone please explain how this function works?

I will attempt to automate this task so that it can be directly imposed in the input files instead of the source files.

Thank you in advance.

[ThomasHowarth]

The eos.TY2Cp function takes 3 arguments, temperature, mass fraction array and a variable to store the calculated specific heat capacity. Presuming that you have compiled with Fuego, you can find a list of all the EOS functions in the PelePhysics submodule in PelePhysics/Source/Eos/Fuego.H . For this particular function, it is a essentially a renaming of the function CKCPBS which is found in the mechanism.H file.

You can also take a look at my pull request #350 to see what I've done for a fixed Lewis number.

[baperry2]

350 is now merged so fixed non-unity Lewis number capability is supported.

To extend on what Thomas said, CKCPBS is analogous to the Chemkin function of the same name. Like all functions in mechanism.H it is machine generated code that gets written by CEPTR. It evaluates heat capacity for an ideal gas based on the NASA7 or NASA9 polynomial format, with the coefficients coming from the thermo data from a mechanism file (originally in Cantera yaml format).

[bazharz]

Thank you for this!

I've reloaded the new version of PeleLMeX so that I can directly impose a Le in the input file. I followed the same steps I had taken to configure the case I'm studying for the old version.

I'm using an O'conaire mechanism; I added it to the Mechanism directory. For the Conaire mechanism, I encountered a compilation error related to Fuego.H. And for the Pele files for reading the .dat file, when I test a mechanism like grimech12 for hydrogen, I encounter a geometry error during execution.

Here's the error I'm facing: Initializing AMReX (97512176d8e5)... MPI initialized with 4 MPI processes MPI initialized with thread support level 0 Initializing SUNDIALS with 1 threads... SUNDIALS initialized. AMReX (97512176d8e5) initialized 0::Assertion rtmp > std::numeric_limits<ParticleReal>::lowest() && iters < maxiters' failed, file "/home/bazharz/Test_PeleLMeX/PeleLMeX/Submodules/PelePhysics/Submodules/amrex/Src/Base/AMReX_Geometry.cpp", line 586 !!! SIGABRT 1::Assertionrtmp > std::numeric_limits::lowest() && iters < maxiters' failed, file "/home/bazharz/Test_PeleLMeX/PeleLMeX/Submodules/PelePhysics/Submodules/amrex/Src/Base/AMReX_Geometry.cpp", line 586 !!! SIGABRT 2::Assertion rtmp > std::numeric_limits<ParticleReal>::lowest() && iters < maxiters' failed, file "/home/bazharz/Test_PeleLMeX/PeleLMeX/Submodules/PelePhysics/Submodules/amrex/Src/Base/AMReX_Geometry.cpp", line 586 !!! SIGABRT 3::Assertionrtmp > std::numeric_limits::lowest() && iters < maxiters' failed, file "/home/bazharz/Test_PeleLMeX/PeleLMeX/Submodules/PelePhysics/Submodules/amrex/Src/Base/AMReX_Geometry.cpp", line 586 !!! SIGABRT See Backtrace.2 file for details See Backtrace.0 file for details See Backtrace.1 file for details See Backtrace.3 file for details

MPI_ABORT was invoked on rank 2 in communicator MPI_COMM_WORLD with errorcode 6.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes. You may or may not see output from other processes, depending on exactly when Open MPI kills them.

I tried to solve it by modifying an existing case like Flamesheet or HITDecay, but I still encounter the same error. Do you have any idea what could be causing the issue? The latest issue I encountered was with fixed_Le. I tried to impose it for an existing Flamesheet case using the drm19 mechanism. When I set, for instance, peleLM.fixed_Le = 0.7, I receive the following error message: Initializing AMReX (97512176d8e5)... MPI initialized with 4 MPI processes MPI initialized with thread support level 0 Initializing SUNDIALS with 1 threads... SUNDIALS initialized. AMReX (97512176d8e5) initialized

================= Build infos ================= PeleLMeX git hash: v23.11-37-g5784bbd-dirty AMReX git hash: 9751217 PelePhysics git hash: c2d6577 AMReX-Hydro git hash: 6f77e32

ParmParse::queryval type mismatch on value number 0 of last occurrence of peleLM.fixed_Le Expected an "i" type which can't be parsed from the string "0.7" peleLM.fixed_Le(nvals = 1) :: [0.7] amrex::Abort::SIGABRT ParmParse::queryval type mismatch on value number 0 of last occurrence of peleLM.fixed_Le Expected an "i" type which can't be parsed from the string "0.7" peleLM.fixed_Le(nvals = 1) :: [0.7] amrex::Abort::SIGABRT ParmParse::queryval type mismatch on value number 0 of last occurrence of peleLM.fixed_Le Expected an "i" type which can't be parsed from the string "0.7" peleLM.fixed_Le(nvals = 1) :: [0.7] amrex::Abort::SIGABRT ParmParse::queryval type mismatch on value number 0 of last occurrence of peleLM.fixed_Le Expected an "i" type which can't be parsed from the string "0.7" peleLM.fixed_Le(nvals = 1) :: [0.7] amrex::Abort::SIGABRT See Backtrace.3 file for details See Backtrace.1 file for details See Backtrace.2 file for details See Backtrace.0 file for details

Could you please assist me in resolving these compilation and execution issues? For my setup, I'm using the O'conaire mechanism, reading temperature fields, and mass fractions of 9 species.

[ThomasHowarth]

So to use the fixed_Le feature, you turn it on with peleLM.fixed_Le= 1, and then give the Lewis number with peleLM.Lewis= 0.7 (see Transport controls). For the other issues, it's going to be tricky to debug without seeing your inputs file, how you generated the mechanism, how you created your pmf file and how you may have adjusted the problem file. Was the FlameSheet where you started from? If you encountered the problem with the HITdecay problem, even if you messed up the the pmf or the chemistry building, seeing this error suggests something is fundamentally wrong with your inputs file and geometry setup there since this problem bypasses all of the chemistry and scalar parts of the code.

[bazharz]

Thank you for the clarification regarding Lewis. Regarding the execution issue, I had actually modified the 3 files: pelelmex_prob.cpp, pelelmex_prob.H, pelelmex_prob_parm.H In order to define new arrays for reading and initializing data from an input .dat file. Normally, in the old version, compilation and execution worked fine. When I mentioned Flamesheet and HITdecay, it's because they were reference files for me to understand how to define my test case. I could share the 3 files with you. As for the mechanism, I'm having trouble identifying the source of the problem. The approach I followed for the previous version of PeleLMeX was to create a directory in Mechanism, and then I placed the .yaml, therm.dat, and tran.dat files in the directory, and it worked fine. Thank you in advance.

[ThomasHowarth]

Did you run the convert script that is required to create the mechanism.H file? Pele does not read the yaml files directly, PelePhysics has a converter called CEPTR (CEPTR) that will take your yaml to create a mechanism.H file.

To impose a flamelet you've made with a particular mechanism, you should use the PMF approach that is used in FlameSheet. To create a PMF file that can be read by the code, use the Utility also given in PelePhysics. See Utility and look at the premixed flame initialisation section.

[bazharz]

Yes, I've indeed run the convert script. I did it for a second time, and I still encounter the same compilation error when using this O'Conaire mechanism! Here's an example of the compilation error I'm facing: /home/bazharz/Test_PeleLMeX/PeleLMeX/Submodules/PelePhysics/Source/Eos/Fuego.H: In static member function ‘static void pele::physics::eos::Fuego::RTY2WDOT(amrex::Real, amrex::Real, const Real, amrex::Real)’: /home/bazharz/Test_PeleLMeX/PeleLMeX/Submodules/PelePhysics/Source/Eos/Fuego.H:213:18: error: ‘mw’ was not declared in this scope 213 | WDOT[n] = mw(n); | ^~ /home/bazharz/Test_PeleLMeX/PeleLMeX/Submodules/PelePhysics/Source/Eos/Fuego.H: In static member function ‘static void pele::physics::eos::Fuego::Y2WBAR(const Real, amrex::Real&)’: /home/bazharz/Test_PeleLMeX/PeleLMeX/Submodules/PelePhysics/Source/Eos/Fuego.H:426:22: error: ‘imw’ was not declared in this scope 426 | summ += Y[i] imw(i); | ^~~ In file included from /home/bazharz/Test_PeleLMeX/PeleLMeX/Submodules/PelePhysics/Source/Eos/EOS.H:14, from /home/bazharz/Test_PeleLMeX/PeleLMeX/Submodules/PelePhysics/Source/PelePhysics.H:8, from /home/bazharz/Test_PeleLMeX/PeleLMeX/Source/PeleLMeX.H:18, from /home/bazharz/Test_PeleLMeX/PeleLMeX/Source/main.cpp:1: /home/bazharz/Test_PeleLMeX/PeleLMeX/Submodules/PelePhysics/Source/Eos/SRK.H: In member function ‘void pele::physics::eos::SRK::RTY2C(amrex::Real, amrex::Real, const Real, amrex::Real)’: /home/bazharz/Test_PeleLMeX/PeleLMeX/Submodules/PelePhysics/Source/Eos/SRK.H:582:13: error: ‘mw’ was not declared in this scope 582 | (mw(ii) InvRuT) * | ^~

I've used PMF for reading initial fields such as temperature and mass fractions. Below is the initial file, pelelmex_prob.H

ifndef PELELM_PROB_H

define PELELM_PROB_H

include

include

include

include

include

include

include

include

// ----------------------------------------------------------- // Search for the closest index in an array to a given value // using the bisection technique. // INPUTS/OUTPUTS: // xtable(0:n-1) => array to search in (ascending order) // n => array size // x => x location // idxlo <=> output st. xtable(idxlo) <= x < xtable(idxlo+1) // ----------------------------------------------------------- AMREX_GPU_DEVICE AMREX_FORCE_INLINE void pelelmex_initdata( int i, int j, int k, int is_incompressible, amrex::Array4 const& state, amrex::Array4 const& aux, amrex::GeometryData const& geomdata, ProbParm const& prob_parm, pele::physics::PMF::PmfData::DataContainer const pmf_data) { const amrex::Real prob_lo = geomdata.ProbLo(); const amrex::Real prob_hi = geomdata.ProbHi(); const amrex::Real dx = geomdata.CellSize();

AMREX_D_TERM(const amrex::Real x = prob_lo[0] + (i+0.5)dx[0];, //x = prob_lo[0] + (i+0.5)dx[0];, const amrex::Real y = prob_lo[1] + (j+0.5)dx[1];, //y = prob_lo[1] + (j+0.5)dx[1];, const amrex::Real z = prob_lo[2] + (k+0.5)dx[2];); //z = prob_lo[2] + (k+0.5)dx[2];);

AMREX_D_TERM( const amrex::Real Lx = prob_hi[0] - prob_lo[0];, const amrex::Real Ly = prob_hi[1] - prob_lo[1];, const amrex::Real Lz = prob_hi[2] - prob_lo[2];)

// Fill in the velocities

int Nx = prob_parm.input_nx; int Ny = prob_parm.input_ny; int Nz = prob_parm.input_nz;

// // Fill Velocity // state(i,j,k,VELX) = 0.0; //0.0; state(i,j,k,VELY) = 0.0; state(i,j,k,VELZ) = 0.0;

if ( is_incompressible != 0) { return; }

state(i,j,k,TEMP) = prob_parm.d_Tinput[i + Nx j + Nx Ny * k]; //amrex::Print() << "i =" << i << "j =" << j << "k =" << k << "state(i,j,k,TEMP) =" << state(i,j,k,TEMP) << "\n";

amrex::Real Yt[NUM_SPECIES] = {0.0}; amrex::Real massfrac[NUM_SPECIES] = {0.0}; Yt[H_ID] = prob_parm.d_Y_H[i + Nx j + Nx Ny k];
Yt[H2_ID] = prob_parm.d_Y_H2[i + Nx j + Nx Ny k]; Yt[O_ID] = prob_parm.d_Y_O[i + Nx j + Nx Ny k]; Yt[O2_ID] = prob_parm.d_Y_O2[i + Nx j + Nx Ny k]; Yt[OH_ID] = prob_parm.d_Y_OH[i + Nx j + Nx Ny k]; Yt[H2O_ID] = prob_parm.d_Y_H2O[i + Nx j + Nx Ny k]; Yt[HO2_ID] = prob_parm.d_Y_HO2[i + Nx j + Nx Ny k]; Yt[H2O2_ID] = prob_parm.d_Y_H2O2[i + Nx j + Nx Ny k]; Yt[N2_ID] = prob_parm.d_Y_N2[i + Nx j + Nx Ny * k]; amrex::Print() << "Yt[H2] = " << Yt[H2_ID] << "\n";

auto eos = pele::physics::PhysicsType::eos();
eos.X2Y(Yt, massfrac); amrex::Real rho_cgs, P_cgs; P_cgs = prob_parm.P_mean * 10.0;

eos.PYT2R(P_cgs, Yt, state(i,j,k,TEMP), rho_cgs); //Yt, state(i,j,k,TEMP), rho_cgs); state(i,j,k,DENSITY) = rho_cgs * 1.0e3; // CGS -> MKS conversion

eos.TY2H(state(i,j,k,TEMP), Yt, state(i,j,k,RHOH)); //Yt, state(i,j,k,RHOH)); state(i,j,k,RHOH) = 1.0e-4 state(i,j,k,DENSITY); // CGS -> MKS conversion

for (int n = 0; n < NUM_SPECIES; n++) { state(i,j,k,FIRSTSPEC+n) = Yt[n] state(i,j,k,DENSITY); //Yt[n] state(i,j,k,DENSITY); } }

AMREX_GPU_DEVICE AMREX_FORCE_INLINE void bcnormal( const amrex::Real /x[AMREX_SPACEDIM]/, const int m_nAux, amrex::Real s_ext[NVAR], const int idir, const int sgn, const amrex::Real time, amrex::GeometryData const& geomdata, ProbParm const& prob_parm, pele::physics::PMF::PmfData::DataContainer const pmf_data) { const amrex::Real* prob_lo = geomdata.ProbLo(); }

AMREX_GPU_DEVICE AMREX_FORCE_INLINE void zero_visc (int i, int j, int k, amrex::Array4 const& beta, amrex::GeometryData const& geomdata, amrex::Box const& domainBox, const int dir, const int beta_comp, const int nComp) { amrex::ignore_unused(i,j,k,beta,geomdata,domainBox,dir,beta_comp,nComp); // We treat species when beta_comp == 0 and nComp == NUM_SPECIES // otherwise this routine could be called for other face diffusivity (Temp, velocity, ...) }

endif

[ThomasHowarth]

I have attached a version of the O'Conaire mechanism that I created using CEPTR - I am not entirely sure why the version you have produced code that doesn't work. Maybe compare your yaml and mechanism.H file to what I have attached.

Additionally, there is a PMF folder that contains a python script to create the appropriate dat files (and an example dat is also there) to read in using the PMF workflow. Take a look at the FlameSheet example to see how to use this file - I would recommend doing this over writing your own code to read in flamelet profiles. oconaire.zip

[bazharz]

Thank you, now it's working with the new version of the mechanism. The Conaire mechanism I'm using does not include Ar and only consists of 9 species. Attached, you will find the conaire.yaml file used for generating the mechanism.H and mechanism.cpp. Unfortunately, I still encounter the same execution error.

Regarding the initial data, I'm not generating a temperature field from a 1D profile; instead, I'm reading the field similar to the HITDecay case. Attached, you will find the 3 pelelm files. case.zip

[drummerdoc]

random question....any chance your mech file has cr/lf's from DOS?

[ThomasHowarth]

Also, could you attach your inputs file?

[baperry2]

@bazharz rather than attaching or copy/pasting files here, it might be easier to just provide a link to a git branch with all the files for your case setup. In the above thread, it's a bit difficult to follow and understand the issues you are having due to the formatting.

[baperry2]

Regarding the issue you were having with the mechanism, are you using the latest versions of PeleLMeX and PelePhysics/CEPTR? I can generate mechanism.{H,cpp} files from your yaml file and compile just fine using the latest versions of PeleLMeX and PelePhysics/CEPTR. Some stuff changed with how we store/compute the molecular weights in https://github.com/AMReX-Combustion/PelePhysics/pull/385 and that appears to be what is tripping your compilation up.

[bazharz]

@drummerdoc I actually have no idea! You can find a new branch created for the case I'm trying to fix in the provided link. @baperry2, yes, I am using the latest version of PeleLMeX, which is why I'm encountering this execution problem in particular.

[ThomasHowarth]

In your inputs file you have amr.max_level = 100. You don't have any amr criteria set, but if I do the same thing for one of my problems I have the error you see. If you don't want amr, set amr.max_level = 0.

[bazharz]

Thank you! It's working well now ! However, I still encounter a compilation error when I choose the mechanism from the Conaire model, which I'm employing with only 9 species instead of the usual 10. I obtained the mechanism.yaml file from the .cti file using the script provided by Cantera, cti2yaml.py. After that, I executed the convert.sh script.

[baperry2]

Not sure why the mechanism conversion is not working for you. Can you provide the git hash of PelePhysics and details of the computing system you are using to convert the mechanism, as wall as the mechanism.{H,cpp} files that you obtain?

In any event, this version I converted from your yaml file should work: https://github.com/baperry2/PelePhysics/tree/bazharz-mech/Mechanisms/oconaire

[bazharz]

Thanks, it's working now! Usually, I copied the YAML file and then ran the convert.sh script to generate the mechanism{.H,.cpp} files. You can find here the old version the mechanism{.H,.cpp} files.