Reaction integration related questions in PeleLMeX #422
-
|
Hi. I am trying to understand the algorithm of PeleLMeX, and I have some questions. 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:
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.
Replies: 11 comments 44 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.
-
|
Yes sure. Will consider, although pretty sure that I have to maintain the confidentiality of the DeepONet code. But will ask my advisor. Thanks again. |
Beta Was this translation helpful? Give feedback.
-
Hello, [Y, T]_n+1 = CVODE([Y, T]_n) + F_AD * dt Thanks. |
Beta Was this translation helpful? Give feedback.
-
Hello, |
Beta Was this translation helpful? Give feedback.
-
|
Not all of this is in the PeleLMeX docs, but if you look at the PeleLM docs there's some more details on the temporal integration approach, which is mostly the same between the two codes: https://amrex-combustion.github.io/PeleLM/manual/html/Model.html#the-pelelm-temporal-integration In CVODE, we integrate the system |
Beta Was this translation helpful? Give feedback.
-
Thankyou for your reply, (referenced in the final code update as shown below // Convert mass fractions back to species densities (with forcing) May you please let me know this as it is like explicit treatment of advection diffusion terms. |
Beta Was this translation helpful? Give feedback.
-
|
Hello, Thank you very much for your help! |
Beta Was this translation helpful? Give feedback.
-
|
I believe that your equation above is correct. The F_AD term is treated as temporally constant during the ODE evolution step. That is then consistent with how the effects of chemistry are extracted from the solution after the chemistry update (that is, it is assumed that the contribution of advection and diffusion during that update is F_AD * dt. During this step, density can indeed change. This part of the algorithm is fundamentally different than most split integrators that I am aware of, but it allows for a reasonable steady balance between A+D and chemistry in the cases where one should exist. With traditionally split approaches a quasi steady system typically is found when the solution is bouncing back and forth between two states, one after the AD step, and one after the chemistry. Also, I see no reason why your DeepONet solution wouldn't be able to serve as an ODE replacement (assuming you get all the pieces in place correctly). Another area to focus on carefully is how the chemical production terms appear in the right hand side to the DivU constraint. There are subtleties in algorithm performance that are different when using an instantaneous evaluation of the production rates vs the integrated average of the production over time. Try to keep a distinction in evaluation capability between the two and use them consistent with the default strategy in PeleLMeX, at least at first. Once you get a baseline functioning, you can explore algorithmic options to understand how they perform, especially when using chemical models of vastly different sizes and numerical stiffness. |
Beta Was this translation helpful? Give feedback.
-
|
Hello, Here RY2RRinvY line should be there after the species density update using deeponets and forcing term. Then the updated densities and mass fractions will be used for further analysis. One more doubt I have. First we have RY2RRinvY(rhoY ^n, rho^n, rho_inv^n, y^n). Then we update the states as shown below, |
Beta Was this translation helpful? Give feedback.
-
|
Very glad to hear about your success, and your ability to sort through the very subtle issues! |
Beta Was this translation helpful? Give feedback.
-
|
Hi. I am trying to run the backward-step flame problem (with Cvode as the chemical integrator) on my university's HPC cluster with GPUs. I am working on a node, which is equipped with an NVIDIA RTX A4000 GPU, and I have compiled PeleLMeX with CUDA and MPI support using the recommended CMake flags. This is how I run the simulation, and the initial output that I see related to MPI and GPU. Is there anything that I am doing wrong? Please let me know. Also, I could not find any tutorial that explains how to compile and run a problem with MPI and CUDA support, so I am not sure if I am doing it correctly. Thanks for your time. |
Beta Was this translation helpful? Give feedback.
-
|
A few more notes on running on GPU:
|
Beta Was this translation helpful? Give feedback.
-
|
Okay got it. Thank you. Also, when using DeepONet as a surrogate chemical integrator, sometimes, the model does not perform well, and it leads to this error: Is there any changes that I can make in the input file to maybe increase the tolerance of these iterations (or something else) which may help in avoiding this error? |
Beta Was this translation helpful? Give feedback.
-
|
From the documentation: First step is to turn up the verbosity. If iterations are converging but slowly, add more. If it converges to a small value, but higher than the tolerance, increase the tolerance. If iterations are diverging, debugging is needed (check that the state looks reasonable, heat capacities are reasonable, etc). |
Beta Was this translation helpful? Give feedback.
-
|
Thank you, and apologies for not checking the documentation first. |
Beta Was this translation helpful? Give feedback.
-
|
That said, we are seeing cases where the deltaT iterations are problematic, and we are trying to think of approaches to make this more robust. So, in the irritating case that the deltaT iterations work fine for a longish time then start to fail, one reason might be that the time step is increasing and eventually crosses a threshold where the deltaT iterations fail to reach the convergence criteria. If this is happening, you should see a correlation with step size and solver robustness. We have seen in that case that you can limit the time step by decreasing the CFL, however we have also seen that this may just delay the onset of the problem. As I said, we are thinking about ways to robustlify this part of the algorithm, so if you can experiment a bit with these settings, and if you find that you are still failing in the modes discussed, let us know and try to create a reproducer that we can use to test any improvements we make. Thanks for your patience and communication on this. |
Beta Was this translation helpful? Give feedback.
Not all of this is in the PeleLMeX docs, but if you look at the PeleLM docs there's some more details on the temporal integration approach, which is mostly the same between the two codes: https://amrex-combustion.github.io/PeleLM/manual/html/Model.html#the-pelelm-temporal-integration
In CVODE, we integrate the system$\frac{d\xi}{dt} = \dot{\omega}_{\xi}(T) + F$ from $t$ to $t+ \Delta t$ . That means the state after CVODE integration includes the effects of reaction, diffusion, and advection, so F_AD does not need to be added separately for the updated state. In fact, to compute I_R, the forcing just due to reactions, we have to subtract F_AD from the update. The density can change during …