Reaction integration related questions in PeleLMeX #422
Replies: 9 comments 25 replies
-
|
In PelePhysics, the state vector that is given to CVode or the other chemical integrators is the species and temperature (rather than enthalpy), which simplifies the Jacobian. See here: https://amrex-combustion.github.io/PelePhysics/CvodeInPP.html#the-different-reactors (although don't read too far because some of that documentation has fallen very out of date). The bits of code you've called out ensure that the initial temperature for chemistry integration is consistent with the initial enthalpy and that the enthalpy gets updated properly at the end of chemistry integration. |
Beta Was this translation helpful? Give feedback.
-
|
Okay. To replace the reaction integration in PeleLMeX (currently done using CVODE) with the DeepONet, I have an intermediate C++ file (with one of the headers being #include <Python.h> ) that calls the necessary functions from the DeepONet code (in python) and returns the required outputs, i.e., species MF and temperature. I have added this function in the ReactorCvode.cpp file (in the CPU region, after the box_flatten, and before the box_unflatten) and included the intermediate C++ file in the header accordingly. Even though I have added the header files and libraries in the Cmake file. I added in BuildPelePhysicsLib.cmake Can you please help with this? Thanks for your time. |
Beta Was this translation helpful? Give feedback.
-
|
As a general word of advice, I would recommend adding a new I haven't previously written software that calls python functions from C++ code (the reverse is more common) so I can't provide useful advice on how to approach that or get things to compile properly. Also, it seems likely to me that this strategy will not end up being computationally efficient. If you are training your models using pytorch or tensorflow in python, it isn't that difficult to link Pele to C++ versions of these libraries and execute a network that was exported from your python training scripts. There are some very basic examples of doing that here: https://github.nrel.gov/bperry/ML-Models-CPP, and we have previously integrated both tensorflow and pytorch models into Pele in this manner. I could provide further pointers if you're interested in this strategy. A final option would be to take advantage of the home-rolled neural network evaluator that is now in PelePhysics (https://github.com/AMReX-Combustion/PelePhysics/blob/development/Source/Utility/BlackBoxFunction/NeuralNetHomerolled.H), although this implementation is very specific and would have to be adapted for your use case. |
Beta Was this translation helpful? Give feedback.
-
|
Thank you for your advice. I will try to implement them in my work. I tried compiling the Python code and the intermediate C++ file (that calls the python functions) with PeleLMeX. It is compiling now, but there is some issue while linking Python, and that's why the executable is not created. I have added this to the "Make.PeleLMeX" file in the "Exec" folder: I have added the Python files and the intermediate C++ in the folder of the case that I am running (Exec/RegTests/EB_BackwardStepFlame) and in the Submodules/PelePhysics/Source/Reactions directory. This is the error that I get on the terminal: I will highly appreciate any help from your side in my work. If you have any idea how to resolve this, please let me know. Thank you so much. |
Beta Was this translation helpful? Give feedback.
-
|
Have you tried to compile the |
Beta Was this translation helpful? Give feedback.
-
|
Hi. After a lot of trials and errors with the makefiles, the Deep Operator Network is finally integrated with PelePhysics. I tried running ReactEval, and it is working fine. The prediction code is highly optimized through JAX implementation, and this is computationally very efficient. Any idea why this is happening? To give you a background of the work that I did, I made changes in the After flattening of the cell data, I made use of the yvec_d 1D array, and made an array Thank you for your time. |
Beta Was this translation helpful? Give feedback.
-
|
There are several runtime controls that are helpful for debugging MLMG errors. Start by increasing the verbosity: If it looks like the MLMG is converging, but not all the way to the defined |
Beta Was this translation helpful? Give feedback.
-
|
However, if you find that one or more of the solver residuals starts off with a NaN norm (or completely unreasonable) there could be a problem with the linking across languages. Here, this would show up as junk in the right hand side MultiFab for the projection(s). Look for that quantity in the initial plot file that gets generated if you run with init_proj turned off. (another very brute force quick check is to look for NaNs in the MultiFab min/max values written in the _H files in the profile folders.) |
Beta Was this translation helpful? Give feedback.
-
|
Because the failure is occurring during the initial projection, you should also visualize the state populated by pelelmex_initdata to ensure that you are initializing the simulation properly after your changes to PelePhysics. See the changes in #482 for how to do that. Often errors in the linear solvers during the initial projections/iterations are due to problems with initialization. |
Beta Was this translation helpful? Give feedback.
-
|
Hi. I tried increasing the verbosity of nodal and mac projectors, and I tried increasing the Next, I tried with Any idea what I can do to bypass this error? |
Beta Was this translation helpful? Give feedback.
-
|
That initial condition looks reasonable to me. From that output, it looks like the linear solvers are converging, although not quite to the default criterion. You can adjust the tolerances as documented here: https://amrex-combustion.github.io/PeleLMeX/manual/html/LMeXControls.html#linear-solvers I'd recommend something around 1e-9 for rtol and 1e-10 for atoll for both mac and nodal projections. |
Beta Was this translation helpful? Give feedback.
-
|
I did that, and now the linear solver is converging, but it said in the command window that dt is too small and the simulation will stop. So I added Now, the simulation is running without issues, but when i try to visualize any plotfile except the initial one, it looks like this:
Every other variable contour looks exactly like the above one. When I look at the values of variables using spreadsheet view, I see a lot of zeros and NaNs. Any idea why this is happening? |
Beta Was this translation helpful? Give feedback.
-
|
Could be lots of different things...a simple thing to try is a super small time step (1.e-10) and take one step and see if you get non-NaN results. In any case, we need to understand where in the code the NaNs are being generated first. |
Beta Was this translation helpful? Give feedback.
-
|
If a smaller timestep also fails, you can compile with |
Beta Was this translation helpful? Give feedback.
-
|
I tried giving As soon as I did |
Beta Was this translation helpful? Give feedback.
-
|
You need to rebuild the TPL (sundials) in debug mode as well. Try:
|
Beta Was this translation helpful? Give feedback.
-
|
Okay. Running the code with debug mode gave me something like this |
Beta Was this translation helpful? Give feedback.
-
|
Hi. During initial velocity projection, I can see that the velocity goes upto 4e+12. Any idea why this is happening and if that can be an issue? I am sharing the part of the code where I made the changes (CPU region of the |
Beta Was this translation helpful? Give feedback.
-
|
Can you add divu to the plot file and look to see if there is something unexpected in that field? |
Beta Was this translation helpful? Give feedback.
-
|
Okay, will look into that more. But what about the negative temperature values in the ReactEval? |
Beta Was this translation helpful? Give feedback.
-
|
Yeah the negative temperatures in ReactEval is a separate issue from the forcing term. Debugging that first is a good idea - ensure that you get results from ReactEval that match detailed chemistry before moving to PeleLMeX. My recommendation is that you implement the DeepOnet chemistry in its own reactor class (derived from ReactorBase) rather than hacking into ReactorCvode. If you're not actually using Cvode there is a ton of extra complication in ReactorCvode that you are not using and that could cause issues if not properly accounted for. Look at ReactorRK64 as an example (it is much simpler than ReactorCvode). That will help you understand the essential parts of the reactor class and potentially identify why negative temperatures may be coming in. |
Beta Was this translation helpful? Give feedback.
-
|
Yes, will try this now. Thanks. In my current approach, the DeepONet takes species mass fractions as input (not species densities). Before passing them to the model, I compute the mass fractions by dividing the species densities by the total density, which is obtained as the sum of all species densities. The DeepONet then predicts the updated mass fractions. To convert these predicted mass fractions back to species densities, I’ve been using the total density from before the reaction step, assuming it remains unchanged. However, this approach appears to be inaccurate in many cases. Is there a more appropriate way to compute the updated total density after reaction integration before converting the mass fractions back to species densities? |
Beta Was this translation helpful? Give feedback.
-
|
I have some questions about
|
Beta Was this translation helpful? Give feedback.
-
|
Beta Was this translation helpful? Give feedback.
-
Hi. I am trying to understand the algorithm of PeleLMeX, and I have some questions.
So I understand that this code uses an operator splitting technique to solve for transport and reactions separately. In the documentation (Step 2 of the algorithm), it is given that after calculating the advective and diffusive terms for reactions and enthalpy, we end up with this ODE system:
Then, these equations are integrated with the CVODE package in SUNDIALS.
In PeleLMeX_Advance.cpp file, the advanceChemistry function is called. Then that function loops over all the AMR levels and integrates chemistry in each of them. When integrating, PeleLMeX_Reactions.cpp falls the m_reactor -> react object, which takes us into the react (the one with box as one of the inputs) function in the ReactorCvode.cpp file (because our chem_integrator = "ReactorCvode") . Inside that function, before integrating, the data is flattened from multidimensional to 1D using box_flatten function, and after the reaction integration has taken place, the the data is unflattened using the box_unflatten function.
My questions are:
During flattening (using box_flatten function defined in ReactorUtils.H), why is the temperature updated using the EOS?
During unflattering (using box_unflatten function defined in ReactorUtils.H), which happens after the reaction ODEs have integrated, why is the value of energy updated again using the frcEExt term, which is essentially the advection and diffusion term related to energy? The advection and diffusion terms were already considered while integrating the reactions, right?
I understand that after energy calculations, we need to calculate temperature, which is done here using the EOS.
Thank you. Waiting for your response.
Beta Was this translation helpful? Give feedback.
All reactions