Pressure driven flows

[Miguel Beneitez]

Hi,

I am new to Dedalus and I think that it is a great tool! I wanted to ask the list since I am trying a very basic flow configuration, but I have a few questions about my implementation:

I am trying to simulate a 2D plane channel flow in primitive variables. It is possible to have Poiseuille flow (and equally pipe flow) considering a fixed mass flux (or bulk velocity). In most codes such as NEK5000, (https://github.com/Nek5000/Nek5000)  this condition is enforced in a correction.

Since in Dedalus the equations are introduced in symbolic form. The way I have implemented this can be found below. The configuration works and the mass flux is indeed kept constant, but I feel that it does not really fit the way Dedalus works.

In particular regarding the variable 'Fdp', the additional eq.

problem.add_equation("Fdp = (integ(u,'y')*0.5-cmflux)/timestep")

mixes the symbolic formulation in Dedalus with the numerical considerations of the body force.

Are you aware of any way of implementing this better? I have tested adding a constraint of the form

problem.add_equation("(integ(u,'y')*0.5-cmflux)=0")

but I get problems regarding the amount of equations/bc.

Any help is appreciated,

Best,

Miguel

# ----------------- SET UP ------------------
# Parameters
Lx, Ly = (2.0*np.pi, 2.0)
nx, ny = (128, 128)
Re = 6000.0

timestep = 0.01

# Constant mass flux
cmflux = 2/3

# Create bases and domain
start_init_time = time.time()

# Basis
x_basis = de.Fourier('x', nx, interval=(0, Lx), dealias=3/2)
y_basis = de.Chebyshev('y', ny, interval=(-Ly/2, Ly/2), dealias=3/2)
domain = de.Domain([x_basis, y_basis], grid_dtype=np.float64)

problem = de.IVP(domain, variables=['p','u','v','uy','vy','Fdp'], time='t')
problem.parameters['Re'] = Re

problem.parameters['cmflux'] = cmflux
problem.parameters['timestep'] = timestep

problem.add_equation('dt(u) + dx(p) - 1/Re * Lap(u, uy) + Fdp = - Adv(u, uy)
problem.add_equation('dt(v) + dy(p) - 1/Re * Lap(v, vy) = - Adv(v, vy)

problem.add_equation("Fdp = (integ(u,'y')*0.5-cmflux)/timestep")

problem.add_bc('dx(u) + vy = 0')

problem.add_equation('uy - dy(u) = 0')
problem.add_equation('vy - dy(v) = 0')

problem.add_bc('left(u)  = 0')
problem.add_bc('right(u) = 0')
problem.add_bc('left(v)  = 0')
problem.add_bc('right(v) = 0', condition='(nx != 0)')
problem.add_bc('left(p)  = 0', condition='(nx == 0)')

solver =  problem.build_solver(de.timesteppers.RK443)
solver.stop_wall_time = np.inf
solver.stop_iteration = np.inf
solver.stop_sim_time = 100.0
initial_dt = 0.01
logger.info('Solver built')

[Miguel Beneitez]

Sorry if I was not too concrete, but the specific question is: can one access the timestepper in Dedalus to make a corrections to ensure constant mass flux or is

problem.add_equation("(integ(u,'y')*0.5-cmflux)=0")

the only way to enforce such. Thanks a lot!

Best,

Miguel

[Jeffrey S. Oishi]

Hi Miguel,

Apologies for the delay. I think adding an additional scalar variable (cmflux) is the most Dedalus way to do whjat you want. I have suggested this to others wanting to do constant mass flux in pipes in Dedalus, though I've never done it myself. This adds the cost of one scalar variable (even though you only need the zero mode of that variable) but eliminates splitting errors. It's essentially what we do with the pressure in incompressible flows.

However, I don't quite know what you mean by "access the timestepper". If I read the paper you posted earlier, it suggested that the entire timestep be taken uncorrected and then the correction applied. You can certainly do this in Dedalus by simply performing the correction step in the timestep loop after completion of the "normal" timestep without U (the required mean flow). If you wanted to do something more complicated, that may also be possible, but perhaps you could say a bit more about what it is you need from the timestepper itself.

Jeff

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/67445a5a-31b0-4b9a-9c1d-a00a7dba6a96n%40googlegroups.com.

[Miguel Beneitez]

Hi Jeff,

Thanks a lot for the answer! After given it a bit more thought I think that this is the way I did it. It reduced to add the following lines:

problem.substitutions['lhs_q'] = "integ( beta/Re * Lap(u, uy) + dx(p),'y')"

problem.add_equation('dt(u) + dx(p) - 1/Re * Lap(u, uy) + Fdp = - Adv(u, uy)')

---------------------------

problem.add_equation("Fdp - lhs_q/Ly = 0")

This equation is obtained by integrating the mass and then considering that d/dt of the integral quantity is 0 and that the nonlinear terms do not contribute to the flow rate.

----------------------------

problem.add_bc('left(dt(Fdp))  = 0')

where Fdp is the body force to keep the constant mass flux/bulk velocity. As you say the equation only needs to be considered for the wavenumber 0, so maybe it would be more efficient to write it as (can one do that?):

problem.add_equation("Fdp - lhs_q/Ly = 0", condition='nx==0')

Also the boundary condition below is a dummy one since the actual value at the wall does not really matter. Maybe there is a better choice?

problem.add_bc('left(dt(Fdp))  = 0')

I also about adding the correction in the time looping itself, but it would imply having to timestep the whole system twice, which might not be the most efficient.

Thanks again for the helpful answer and the great job developing Dedalus!

Best,

Miguel

[Jeffrey S. Oishi]

Hi Miguel,

There's probably no point in limiting only to kx == 0, since you're carrying around the whole field anyway. Let me know how this works!

Jeff

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/a7fbc6b6-c94e-4148-aa3a-9ac6f160c66en%40googlegroups.com.

[Miguel Beneitez]

Hi Jeff,

it worked just fine! The flow rate was kept constant. Two minor details: 1. the value that is kept constant is that of the initial condition, so that is the imposed value. 2. adding the forcing explicitly as

problem.add_equation("Fdp = lhs_q/Ly ")

blows up after a few iterations, probably due to the lag in the value of the forcing.

Best,

Miguel

[Keaton Burns]

Hi Miguel,

I think you should be able to have the mean pressure drop/gradient (Fdp) determining completely implicitly by having a boundary condition directly simply set the flux (integ(u,'y') = Q), similar to what you mentioned in your first email, rather than needing the lhs_q variable. This should be very performant since you could set the Fdp field to be constant in y (so it only adds a single number to each solve).  Then as Jeff says you shouldn't have to do any splitting or adjustment during timestepping, or calculate additional complex substitutions.  So all the ingredients would be:

1) Add Fdp to the solver as a variable.

2) Set "problem.meta['Fdp']['y']['constant'] = True" so only one Chebyshev mode is added to each pencil

3) Add "problem.substitutions['Q'] = Q" to pass in the value of the mass flux you want.

4) Add "+ Fdp" to the u-momentum equation.

5) Add "problem.add_equation("Fdp = 0", condition="nx != 0")" to eliminate the horizontally-varying modes of Fdp.

6) Add "problem.add_equation("integ(u,'y') = Q", condition="nx == 0")" to fix the mass flux and implicitly determine Fdp.  Note Q has to be on the RHS.

-Keaton

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/3ae45314-7eff-4dd7-9e80-fc824d6bd55an%40googlegroups.com.

[Miguel Beneitez]

Hi Keaton,

Thank you for the detailes indications. It works perfectly. It was in line with my first version of the implementation, but I did not know about how to make it that flexible with Dedalus. Your suggestion for the implementation was about 20% faster than the version where I was carrying the whole Fdp field.

Best,

Miguel

[Miguel Beneitez]

Hi again,

Should the implementation suggested work no matter the basis? When I change the Chebyshev basis of the problem from

y_basis = de.Chebyshev('y1', ny, interval=(-Ly/2, Ly/2), dealias=3/2)

to

yb1 = de.Chebyshev('y1', int(ny/2), interval=(-Ly/2, 0), dealias=3/2)
yb2 = de.Chebyshev('y2', int(ny/2), interval=(0, Ly/2), dealias=3/2)
y_basis = de.Compound('y', (yb1, yb2))

I get the following error which is not shown when running the flow on a single basis. I have checked in the group, but I could not find a similar issue being reported. Have you experienced anything similar? Thanks a lot for the time!

  File "fullprob_prim.py", line 355, in <module>
    solver.step(dt)
  File "/home/mikibene/work/Numerics/codes/anaconda3/envs/dedalus/lib/python3.8/site-packages/dedalus/core/solvers.py", line 507, in step
    self.timestepper.step(self, dt)
  File "/home/mikibene/work/Numerics/codes/anaconda3/envs/dedalus/lib/python3.8/site-packages/dedalus/core/timesteppers.py", line 605, in step
    pX = p.LHS_solvers[i].solve(pRHS)  # CREATES TEMPORARY
  File "/home/mikibene/work/Numerics/codes/anaconda3/envs/dedalus/lib/python3.8/site-packages/dedalus/libraries/matsolvers.py", line 121, in solve
    return self.LU.solve(vector, trans='T')
ValueError: b is of incompatible size

Best,

Miguel

onsdag 11 augusti 2021 kl. 18:47:12 UTC+1 skrev keaton...@gmail.com:

[Miguel Beneitez]

I attach here a minimal example of the problem. The code with just 1 domain works fine, but the compound domain causes the error above in Dedalus, could this be related to the property "problem.meta['Fdp']['y']['constant'] = True" in the y direction and the way it is treated in compound domains?

Best,

Miguel

[Keaton Burns]

Hi Miguel,

Yes it should work with any basis, but it looks like there may be a bug in how the constant flag is applied to the compound basis. I’ll try to debug this later this week, but in the meantime does it work if you leave out the constant flag and add the equation “dy(Fdp) = 0”?

Keaton

On Aug 17, 2021, at 5:47 AM, Miguel Beneitez <miki...@gmail.com> wrote:



To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/8431e6c7-2790-4d7b-bc07-f9344ec5f809n%40googlegroups.com.

<channel_compound.py>

<channel_1domain.py>

[Miguel Beneitez]

Hi Keaton,

Great! Thanks  for looking into it. Removing the flag and adding the equation does not work, since it gets imbalanced as follows:

# Equation to ensure constant bulk velocity

problem.add_equation('dy(Fdp) = 0')
problem.add_equation('Fdp = 0', condition = 'nx != 0')
problem.add_equation("integ(u,'y') = Ub*Ly", condition = 'nx==0')

without the constant flag gives the error

ValueError: Pencil (1,) has 7 non-constant equations for 6 non-constant variables.

Best,

Miguel

[Keaton Burns]

Hi MIguel,

Sorry that equation will also need the condition nx==0.

-Keaton

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/da58bba5-de1e-430b-8e9b-32d51f223f70n%40googlegroups.com.

[Miguel Beneitez]

Hi Keaton,

It works fine adding the condition. Thanks a lot. I will check the performance a bit and get back. Thanks again for all the support.

Best,

Miguel