Varying inlet composition with PMF

[asmunder]

I'm running simulations using PMF to set the initial and boundary conditions. This works well, but now I'd like to vary the composition at the inlet as a function of time during the simulation.

I can see a way to implement this functionality: basically just copying the PMF stuff and modifying it to read a second file that describes inlet composition as a function of time.

Just wanted to check if this sounds like a good approach, and that I'm not reinventing the wheel etc?

If it sounds good, I'll submit a PR when I have something that works.

[drummerdoc]

Seems like you’re doing too much. The boundary functions only need the ii flow value and you have the value of time to work with. So you can define any function of time you’d like. Every user and every case will probably be different so I would get too fancy. Maybe I misunderstand what you wanna do?

[emotheau]

You can simply adapt the the bc_fill_nd.F90 file to impose a time varying BC. For example to introduce a pulse of temperature of 100Hz you can do: eos_state % T = T_ref + 2.0d0 sin(2.0d0M_PI100.0d0time)

You can do a similar thing with eos_state % molefrac(1:nspecies) to change the composition over time.

[drummerdoc]

Oh wait...I re-read the comment and think I understand the question now. The answer, I believe, is a variation of what we both said above, but a clarification might help.

The Pele codes (like all AMReX codes) have been written in a way that divides up "general purpose" and "case-specific" bits. Within the case-specific bits, it's still the case that many problems people do look very similar in a generic sense. Over the years, I have found it useful to have a "bcfunction" type thing, for example to manage self-consistently specifying all the states for each boundary. While that's handy, it's a bit restrictive, because it shoves all the data for each boundary into module data as constants. Users will undoubtedly have needs that stretch beyond what I planned for, such as needing a time-dependent function rather than a constant, which requires more data.

The general boundary function interface to the code should cover most other cases - that's how turbulent inlet data is managed, for example. Users will occasionally find that they have to rig up something using the more general stuff (or bypass the "general" stuff with even more brute force hacking). Sounds like you are somewhere in that intermediate area. The real questions are:

1. Has someone done this particular variant before, and
2. Would anyone else find it useful

I'm am not aware of anyone needing that particular setup (yet), so you're on your own to design what best suites your needs. But second, it is probably a bit too specifically geared toward your case to add to the demos.

More generally, the question you appear to really be asking is associated with where "user" code should be stored. That's interesting to think about. Note that I have moved to putting much of my work into a separate repository called "PeleProduction", and this accomplishes a couple of things:

1. I can manage the set of repositories to all be consistent with each other for a particular study. For example, as we hack out the EB stuff in PeleLM, I can work with the correct set of commits across all the repos as they evolve, but can also easily switch over to reproduce a user's issues with some other branches.
2. I can manage my build/run folders outside of the main code for my own work.

Maybe you might consider a similar strategy? Maybe you can try it, and if it works well, we can suggest it for people in general?

Let me know if this became too rambling, and I lost you.

[asmunder]

Thank you both for the answers. As you say, having temperature/pressure/density/velocity vary with time is straight-forward. This works as it should, see gif below (hopefully working) with oscillating temperature at the inlet at 100 Hz, 750 +/- 50 K.

To be more specific on what I am trying to achieve next, is that I want to vary inlet equivalence ratio as a function of time. Currently I work with Cantera in Python to generate my PMF files, and this is very nice (I have it in a Jupyter Notebook, that plots everything and saves it to a file).

I could of course implement varying equivalence ratio myself just as functions modifying molefrac, but especially if my fuel is a mixture I'd have to recreate a lot of the logic already in Cantera, so I'd like to generate a file with Cantera and just read that.

I do already work with my own "user" code in a separate repo, but I don't have it linked to specific versions of the PeleC repos. I will look into how you do that, thanks for the tip Marc.

[drummerdoc]

those are cool little movies.
Check the README in the PeleProduction repo for a discussion about git submodules.

[emotheau]

Ok I understand your first question. So basically if I'm correct, you want to read through multiple PMF files that you have generated (and that contain the oscillation), and interpolate the solution in time to correspond to the physical time at t^n. If I'm correct I think it is similar to what we do to inject turbulence.

If you have developed a feature that does that, I'm definitely interested to put that in a more generic way in the NSCBC stuff. This is also a feature that could be useful in PeleLM.

[asmunder]

I have been working on this for a bit in the past couple of months. I have something that works now, which is strongly based off the PMF stuff. I'll outline the details below.

If you think this is sufficiently useful to be integrated in PeleC and PeleLM, what is the best approach for proceeding? Should I submit a PR for PeleC with an example case in Exec/RegTests/PMF that demonstrates this? The functionality is mainly implemented in pmf_generic.F90 and bc_fill_nd.F90 , there are no changes to the "core" PeleC code.

Basically in this approach there is only one additional input file, which contains the inlet variation in temperature, density, velocity, and compositions as a function of time. (As opposed to a PMF file, that contains the same variables but as a function of space, for setting the initial condition.)

The way I see it, the outlet doesn't really need this kind of functionality since only pressure (at most) is actually set at the outlet. Thus my code is only for applying variation to the inlet. And it only supports one inlet, but I guess it's pretty easy to extend if one is interested in multiple inlets in some complex geometry.

Below is an animation where I have taken the output from a Cantera computation of an autoignition flame, started the simulation with the entire domain filled with the Cantera-predicted combustion products and predicted final temperature as initial condition, and then slowly change the inlet from the products to the reactants. There is no perturbation of the inlet here, the flame is just very noisy. I've had trouble initializing these properly in the past, but this approach is gentle enough to work.

[drummerdoc]

So, to make sure I understand, what you've done is create a pmf-like database using time as the independent coordinate, and then step through that for the BC using the simulation time?

Yes, this seems handy. We have been throwing around ideas for creating a set of tools to generate boundary conditions. For example, in 2D or 3D it might be nice to place an inlet jet wherever convenient and be able to specify things like jet composition, transition thickness, radius, etc as part of that. I suppose the turbulent inlet stuff would be another such thing. And, your function could be something available to compose inlet mixtures on top of this (that is, what if we want a round jet issuing fuel whose mixture varies over time? This seems really useful although I have to admit that I haven't thought of a good organization for that yet.

I wonder if we can stuff these things somewhere for now and collect them up when we think of a clever way to organize them. Let me think about it. If you have ideas, post them here.

[asmunder]

Yes, that is exactly what I have done here.

Just to provide my two cents on the general case (and apologies for the long-winded-ness):
For "organization" of boundary conditions, there are several aspects of OpenFOAM's approach that I like. You have one tool to generate a set of patches that together cover the entire boundary of the domain. Then you have a set of routines for different types of boundary conditions, each taking different parameters and/or input files. Finally you can combine these two sets in an arbitrary way. With this approach, you benefit from the "explosion of possibilities" that arises when you take the product of two sets.

In practice, I could imagine a solution where the first tool to generate a set of patches takes an input file that specifies a set of geometric shapes (boxes, cylinders, spheres, etc.) and boolean operators on these, and then creates a .cpp file analogous to PeleC_init_eb.cpp. The input file would specify names for each patch, like circular_inlet, outlet, outer_wall, bluff_body etc. that the .cpp file then uses in its definitions, such that they become accessible from inputs.3d after compilation.

Then you would have a generic bc_fill_nd.F90 and related routines that define the generic Interior, Symmetry, SlipWall, TimeVaryingInletComposition, TurbulenceFeed, IsothermalNoSlipWall etc. And you would specify in inputs.3d something like

pelec.use_bc_patches = 1
circular_inlet.bc_type = "TimeVaryingInletComposition"
circular_inlet.input_file = "inlet_data.dat"
outer_wall.bc_type = "IsothermalNoSlipWall"
outer_wall.Temp = 750.
....

Of course there are many open questions here. For instance, should you have just one bc_type for each patch, and then you need to write in your code the full set of combination of inlet types like IsothermalNoSlipWall and TurbulentVelocityTimeVaryingCompositionInlet? Or should you do the alternative which is less coding but more verbose and redundant for the input file (which is like what OpenFOAM actually does), and have in your inputs.3d be something like the following?

pelec.use_bc_patches = 1
circular_inlet.bc_type.T = "TimeVaryingInlet"
circular_inlet.bc_type.U = "TimeVaryingInlet"
circular_inlet.bc_type.rho = "TimeVaryingInlet"
circular_inlet.bc_type.composition = "TimeVaryingInlet"
circular_inlet.input_file = "inlet_data.dat"
outer_wall.bc_type.T = "Isothermal"
outer_wall.bc_type.U = "NoSlip"
outer_wall.bc_type.rho = "ZeroGradient"
outer_wall.bc_type.composition = "ZeroGradient"
outer_wall.isothermal.Temp = 750.
....

I don't know if this fits with the overall vision, whether it would be hard to make it play nicely with EB, etc. But I think it could be an approach that would be nice and flexible from the user perspective, at least.

[drummerdoc]

Thanks for taking the time to think about this and write it out. This topic is an example of how a code becomes hugely more complex when the domains start looking less like idealized boxes - what I have been working with for 25 years now! Let me digest this, and maybe get some feedback from others on my team here, and anyone else with ideas, and try to decide how to move forward.

Again, I do appreciate the effort here. And rest assured that we'll need to do something along these lines eventually, so the input is valuable.

[emotheau]

Currently, the procedure is that there is a multifab of integers named bcMask. This is a mask containing an integer value pointing to the type of boundary condition that you want to apply (wall, inflow, outflow, etc...). When a routine impose_NSCBC_xd is called, the bc_type provided by the user is copied into the bcMask. Thus, if the user puts a wall, this information will propagate in other place of the code, such as the Riemann solver.

The bc_type is provided point wise by the user via the routine bcnormal. This allows the user to construct a combination of bcs on the physical boundaries of the domain.

If I understand correctly your idea, you would like to construct the bcMask via generic objects similar to what is done to create EB. This is a good idea, but the bcMask should not only contains integers pointing to the imposed BC, because it should embed everything. Actually I think it should be a c++ object with multiple members that would contain the bc_params and bc_target provided by the current bcnormal routine. Then in the actual NSCBC routine, these data would be extracted from the object rather than passed through a fortran routine.

That's a good idea, but a significant piece of work in perspective !

[asmunder]

So coming back to this issue, I have written the code to vary inlet composition based on a text file, and I've been running simulations with this for a while.

It seems to work fine, but I have a couple of questions regarding the pressure.

Regarding which pressure is actually set at the outlet This is just to clarify my understanding, but is the pressure we set at the outlet the static pressure p, or the total pressure P = p + 0.5 rho*U^2 ?

Regarding the thermodynamic state at the inlet Currently with NSCBC, the thermodynamic state at the inlet is set using temperature and pressure. The temperature is whatever I set at the inlet, be it constant or varying. But the pressure is currently set based on the outlet pressure, as far as I understand the code. I think this becomes inconsistent if there is a significant pressure drop along the domain - even in the animations I posted above, there is about 0.1 bar pressure drop from inlet to outlet.

I have tried to figure out how to get the pressure just inside the domain at the inlet, but so far I've not been able to accomplish this. The pressure seems to be in the q array, and I can't trace the call sequence of the NSCBC routines back so that I can figure out how to add this array to the argument list.

Is my thinking here correct? If yes, how do I get access to the pressure just inside of the inlet in bc_fill_nd.f90 ?

[drummerdoc]

My guess is that you want to call something like ctoprim on the state data passed in the valid box to the bc filler routine. You can call it pointwise, or over a slab at just inside the exit plane. That would get you p at the cell centers anyway. If you want values extrapolated to faces, it's a much more nuanced question.

I still am interested in the stuff above. I was just talking with a couple of the postdocs about the concept of a Gas object, as in Cantera. The idea is that it holds the entire state of the gas in a zeroD kind of way. You could, for example, use it to set the fuel and oxidizer tank mixtures when setting a mixture fraction definition, or it could be the result of a pmf (or time interp like yours above) result. It should also be able to easily convert between mole and mass fraction, etc. It seems useful and could be a part of what you've done.

What form is all your stuff in, and where is it?

[emotheau]

Hi Asmund, I'm glad that you are still working on the NSCBC stuff and made progress ! For your question, NSCBC is based on characteristic variables, and the pressure is the one coming from the equation of state through the routine ctoprim.

However I'm not sure to understand your second question. When you are imposing an inlet, only the target temperature and velocities are employed to compute the characteristic waves, and there is no target pressure involved in the formulation to impose an inflow. Moreover, once the waves are created and the derivative of the characteristic variables are computed, the values of the characteristic variables in the ghost-cells are reconstructed with the help of the characteristic values that already exist in the cells inside the domain, the ones close to the physical boundary. So I'm confused with what you ask, because I don't see where an outlet pressure comes here? Also the pressure inside the domain is already accessible and used here, so I don't understand why you have difficulty to get it?

By the way, is it possible that you share some results ? I'm curious of the performance of this new inlet.

[asmunder]

Hi both, thanks for the responses.

My stuff is all in files that live in the "run directory", specifically pmf_generic.f90, bc_fill_nd.F90, Prob_nd.F90. Having a Gas object would certainly make this a whole lot easier, but I guess this would also be quite a bit of work.

I'll do a little cleanup and then submit a pull request with this stuff, make it a new case in RegTests probably?

@emotheau : it could very well be that I'm misunderstanding the code here. If you look at line 241 of Exec/RegTests/PMF_NSCBC/bc_fill_nd.F90 in the current development branch, at the low BC with which_bc_type = Inflow it sets eos_state % p = pamb. Then it does call eos_xty(eos_state) and call eos_tp(eos_state), obtaining the density etc. at the pressure pamb that comes from the probin file, and at the composition and temperature set a few lines further up, possibly varying with time. The density thus obtained is then used when setting values in u_ext(:) . Won't this use the density computed at outlet pressure? Also, to clarify, it's here in bc_fill_nd.F90 that I didn't see how to get the pressure inside the domain - but the suggestion of calling ctoprim should probably work.

[emotheau]

Ok I understand the confusion. The bcnormal routine is a "user-defined" defined routine used to fill the ghost-cells during a fill-patch operation. The keep consistency, and to provide data from only one place whenever the nscbc flag is turned on or off, this routine is also called in the GC-NSCBC routine. However U_ext values are not used in the context of NSCBC. Actually, the purpose of GC-NSCBC is to recompute the ghost-cells values that were filled by a fillpatch operation.

If you need to pass values to bcnormal from inside the domain, there is the U_in array, which is currently not used in the GC-NSCBC context (we pass an empty U_dummy array). But I think you could fill it with values from "Q" that contains the characteristic variables in the domain.

I'm sorry that all of this stuff seem confusing, but keep in mind that GC-NSCBC were developed in a separate "experimental" code, and then ported directly into PeleC while trying to conserve as much as possible the original structure of the code.

[marchdf]

Closing for staleness. Please reopen if need be.