# Incompressible BCs not fully enforced

I’m lookig at your ini file now, there are quite a few things wrong with it. Somethings to bear in mind:

1. A pressure of zero will always give issues
2. You didn’t increase the max pseudo iteration number as I recommended
3. You changed the pseudo-dt-fact to 50. This will definitely cause issues, this menas that as you do from the p-multigrid level p to p-1, the pseudo-dt in increasing by a factor of 50. Which will easily make things unstable.
4. You reduced nu to ~1e-6, this will make the Reynolds number very high and, unless you have enough mesh, will lead to instabilities.
5. In the integrate plugin I recommended you calculate the intgral of the absolute divergence, this is an important difference as the integral of divergence often equals zero whereas the pointwise values might not be zero.

I’ll let you know if I get a working inifile for your case

Sorry, I misundersood what you tought me about pseudo iteration number.
I appliciate if you give me ini file.
Thank you.

Thank you so much your supports. I could calculate my case 1 step without NAN.(Since the second and subsequent steps fall no matter how I do it, I may need to devise a way, so keep thinking of ways to make it work.) I attach ini file since I want you to review it in order to calculate it next the second and subsequent steps.)
By the way, in terms of velocity, pressure and nu(kinematic viscosity), What are the units? This is because I only know the no-dimensional NS equation and the units are not opened at User Guide.
Moreover , I would like to know how and where in codes Reynolds number and Mach number are defined.

[backend]
precision = double

[constants]
nu = 15e-3
Uin = 10.0
Vin = 0.0
Win = 0.0
Pc = 1.0
ac-zeta = 2.5

[solver]
system = ac-navier-stokes
order = 3

[solver-time-integrator]
formulation = dual
scheme = sdirk33
pseudo-scheme = rk45
controller = none
pseudo-controller = local-pi
tstart = 0.0
tend = 0.1
dt = 0.0005
pseudo-dt = 0.00005
pseudo-niters-min = 100
pseudo-niters-max = 200
pseudo-resid-norm = l2
pseudo-resid-tol = 5e-4
pseudo-resid-tol-p = 2.5e-2
atol = 1e-1
pseudo-dt-max-mult = 1.0
[solver-dual-time-integrator-multip]
pseudo-dt-fact = 1.75
cycle = [(3, 1), (2, 2), (1, 4), (0, 8), (1, 4), (2, 2), (3, 1)]
[solver-interfaces]
riemann-solver = rusanov
ldg-beta = 0.5
ldg-tau = 0.1

[solver-interfaces-line]
flux-pts = gauss-legendre

[solver-elements-tri]
soln-pts = williams-shunn

[solver-elements-tet]
soln-pts = shunn-ham

soln-pts = gauss-legendre

[solver-interfaces-tri]
flux-pts = williams-shunn

[soln-plugin-nancheck]
nsteps = 1

[soln-plugin-pseudostats]
flushsteps = 1
file = ./results1/csv/residual_merge2.csv

[soln-plugin-writer]
dt-out = 0.001
basedir = ./results1
basename = target_with_fluid_merge2-{t:.2f}

[soln-plugin-integrate]
nsteps = 1
file = ./results1/csv/integral_merge2.csv
int-dive = abs(%(div1)s + %(div2)s + %(div3)s)
[soln-bcs-wall]
type = no-slp-wall

[soln-bcs-outlet]
type = ac-out-fp
p = Pc

[soln-bcs-inlet]
type = ac-in-fv
u = Uin
v = Vin
w = Win

[soln-ics]
u = Uin
v = Vin
w = Win
p = Pc


Best regard,

Hi @Yujif1Aero, worth checking if your Uin value of 10 means your inflow velocity is 10 times the speed of sound. Can you describe the system you are trying to solve (i.e., mesh scale, input velocity, fluid viscosity). Also, might be worth getting this to run using the plain Navier-Stokes solver first.

I have tried to simulate the vortex around buldings.
inpul velocity is more than 10m/s, viscosity is air ( nu = 1.48e-5 )
In terms of mesh, I am not used to make mesh using Gmesh (I do not have commercial mesher eg pointwise) I also want to know how check the mesh scalse at wall but I only know automatic meshing by gmesh.( I am learning how to use gmesh now) If possible I want to calaculate it with mesh Y+=1 ~ 10 at wall.

I attach mesh file at target_with_fluid_merge2.msh - Google Drive

Below I attach ini file which generates NAN at some step,

[backend]
precision = double
rank-allocator = linear

[backend-cuda]
device-id = local-rank
mpi-type = standard

[constants]
nu = 1.48e-5
#nu = 0.1
Uin = 10.0
Vin = 0.0
Win = 0.0
Pc = 101325
ac-zeta = 2.5

[solver]
system = ac-navier-stokes
order = 3

[solver-time-integrator]
formulation = dual
scheme = sdirk33
pseudo-scheme = rk45
controller = none
pseudo-controller = local-pi
tstart = 0.0
tend = 2.0e-3
#dt = 0.01
#pseudo-dt = 0.005
#dt = 0.001
#pseudo-dt = 0.0001
dt = 1.0e-3
pseudo-dt = 1.0e-4
pseudo-niters-min = 50
pseudo-niters-max = 50
pseudo-resid-norm = l2
#pseudo-resid-tol = 5e-4
#pseudo-resid-tol-p = 2.5e-2
#atol = 1e-1
#pseudo-dt-max-mult = 1.0
pseudo-resid-tol = 5.0e-4
pseudo-resid-tol-p = 2.5e-2
atol = 1e-1
pseudo-dt-max-mult = 1.0
#pseudo-dt-max-mult = 2.5

[solver-dual-time-integrator-multip]
pseudo-dt-fact = 1.75
#cycle = [(3, 1), (2, 1), (1, 1), (0, 2), (1, 1), (2, 1), (3, 3)]
cycle = [(3, 1), (2, 2), (1, 4), (0, 8), (1, 4), (2, 2), (3, 1)]
[solver-interfaces]
riemann-solver = rusanov
ldg-beta = 0.5
ldg-tau = 0.1

[solver-interfaces-line]
flux-pts = gauss-legendre

[solver-elements-tri]
soln-pts = williams-shunn

[solver-elements-tet]
soln-pts = shunn-ham

soln-pts = gauss-legendre

[solver-interfaces-tri]
flux-pts = williams-shunn

[soln-plugin-nancheck]
nsteps = 1

[soln-plugin-pseudostats]
flushsteps = 1
file = ./results3/csv/residual.csv

[soln-plugin-writer]
dt-out = 1.0e-3
basedir = ./results3
basename = target_with_fluid_merge3-{t:-2f}

#[soln-plugin-integrate]
#nsteps = 1
#file = ./results/csv/integral.csv

[soln-plugin-integrate]
nsteps = 1
file = ./results3/csv/integral.csv
int-dive = abs(%(div1)s + %(div2)s + %(div3)s)

#[soln-bcs-inlet]
#type = ac-char-riem-inv
#ac-zeta = 180
#p = Pc
#u = Uin
#v = Vin
#w = Win

#[soln-bcs-outlet]
#type = ac-char-riem-inv
#ac-zeta = 180
#p = Pc
#u = Uin
#v = Vin
#w = Win
##p = Pc

[soln-bcs-wall]
type = no-slp-wall

[soln-bcs-outlet]
type = ac-out-fp
ac-zeta = 120
p = Pc

[soln-bcs-inlet]
type = ac-in-fv
ac-zeta = 120
u = Uin
v = Vin
w = Win

[soln-ics]
u = 3.3
v = 3.3
w = 3.3
p = Pc


Two points, you probably want to reduce the pseudo-niters-min, you don’t want to force the system to take 100 steps if it has already converged, 4 is probably more appropriate. For pseudo-niters-max, 100 is probably better too.

PyFR doesn’t assume units, see one of the many other posts asking similar questions.

For example:

Finally on you question about GMSH, this might be better asked to there support team/community.

1 Like

I got results finally. But presure is oscillation and inlet pressure is gradually becoming bigger.
Do you know why it is happed ?

Pressure t = 1.0

Pressure t = 1.9

Velocity = 1.0

Velocity = 1.9

[backend]
precision = double
rank-allocator = linear

[backend-cuda]
device-id = local-rank
mpi-type = standard

[constants]
nu = 0.5e-5

Uin = 1.0
Vin = 0.0
Win = 0.0
Pc = 783.642
ac-zeta = 2.5
rhom = 1.0

R = 287.1

n = 7.0

bt = 7.5
[solver]
system = ac-navier-stokes
order = 3

[solver-time-integrator]
formulation = dual
scheme = sdirk33
pseudo-scheme = rk45
controller = none
pseudo-controller = local-pi
tstart = 0.0
tend = 2.0
dt = 0.01
pseudo-dt = 0.0005
pseudo-niters-min = 4
pseudo-niters-max = 50
pseudo-resid-norm = l2
pseudo-resid-tol = 5e-4
pseudo-resid-tol-p = 2.5e-2
atol = 1e-1
pseudo-dt-max-mult = 1.0

[solver-dual-time-integrator-multip]
pseudo-dt-fact = 1.75
cycle = [(3, 1), (2, 2), (1, 4), (0, 8), (1, 4), (2, 2), (3, 1)]
[solver-interfaces]
riemann-solver = rusanov
ldg-beta = 0.5
ldg-tau = 0.1

[solver-interfaces-line]
flux-pts = gauss-legendre

[solver-elements-tri]
soln-pts = williams-shunn

[solver-elements-tet]
soln-pts = shunn-ham

soln-pts = gauss-legendre

[solver-interfaces-tri]
flux-pts = williams-shunn

[soln-plugin-nancheck]
nsteps = 100

[soln-plugin-pseudostats]
flushsteps = 100
file = ./csv/pseudostats_results5.csv

[soln-plugin-writer]
dt-out = 0.1
basedir = ./results5
basename = target_with_fluid_results5-{t:.2f}
post-action = echo "Wrote file {soln} at time {t} for mesh {mesh}."

[soln-plugin-fluidforce-bwall]
nsteps = 100
file = ./csv/bluding-forces_results5.csv
morigin = (0.0, 0.0, 0.0)

[soln-plugin-tavg]
nsteps = 100
dt-out = 0.1
mode = continuous
basedir = ./results5
basename = files-avg-{t:.2f}

avg-p = p

avg-u = u
avg-v = v
avg-w = w
avg-uu = u*u
avg-vv = v*v
avg-ww = w*w
avg-uv = u*v
avg-uw = u*w
avg-vw = v*w
fun-avg-upup = uu - u*u
fun-avg-vpvp = vv - v*v
fun-avg-wpwp = ww - w*w
fun-avg-upvp = uv - u*v
fun-avg-upwp = uw - u*w
fun-avg-vpwp = vw - v*w
avg-vel = sqrt(u*u + v*v + w*w)
fun-avg-urms = sqrt(uu - u*u + vv - v*v + ww -w*w)

[soln-plugin-integrate]
nsteps = 100
file = ./csv/integral_results5.csv

int-E = rhom*(u*u + v*v + w*w)
int-enst = rhom*(%(vor1)s*%(vor1)s + %(vor2)s*%(vor2)s + %(vor3)s*%(vor3)s)
int-dive = abs(%(div1)s + %(div2)s + %(div3)s)

[soln-bcs-wall]
type = no-slp-wall

[soln-bcs-bwall]
type = no-slp-wall

[soln-bcs-outlet]
type = ac-out-fp
ac-zeta = 180
p = Pc

[soln-bcs-inlet]
type = ac-in-fv
ac-zeta = 180
u = Uin * pow((z / bt), (1.0/n))
v = Vin
w = Win

[soln-ics]
ac-zeta = 180
u = Uin * pow((z / bt), (1.0/n))
v = Vin
w = Win
p = Pc


Inlet surface presure are
t = 1

t = 1.9

t = 0

t = 0.1

In your case, I thought you were simulating a channel with an obstruction? I don’t see the obstuction in the pictures, rather, you seem to have mesh where it should be.

Edit: I looked at the mesh you sent previously and I now understand what your case look like.

The pressure at the inlet isn’t set by the ac-in-fv boundary and is free to change. Pressure waves will occur in incompressible fluids as the flow accelerates, initially you’ll see larger pressure waves as the solution settles down from the unrealistic IC to a more realistic state. You just have to power though this period and it normally takes a couple flow throughs of the domain. For example, given the size of the domain ([-10,10]^2\times[0,20]) and the inlet velocity (u=1), it may take until t=20 to 80 or longer. What you can do initially is reduce the order to order = 1 and increase nu. This does two things, the reduced order speeds up the simualtion and adds diffusion to damp out the waves. Increasing nu also add diffusion. Once the solution has developed you can then return the order and nu to the desired values.

Sorry, I was lack of my explanation. I attach detale figure.
I do not understand why puresure gets unsatable at Inlet. Could you teach me?

Thank you for replyng so much.
I understood I needed to caluclate many period to get stable.
However I could not understand that if nu is increased, the diffusion is increasing.

Oh nu = 1 /Re .
So that if nu is increased, Re is decreased. Is it right?

Yep, changing nu will change the Retynolds number, but this is only for the initialisation. Once you’ve done some flow throughs, you then change it to the target order and nu.

Thank you. Is the order solver’s order? If mesh is made in 3 order, is it impossible?

Yep, it should be fine. The pyfr import doesn’t use the ini file so the pyfrm file is independent of the mesh order. Also looking at your domain, it should have entirely flat surfaces anyway.

1 Like

When 3-order mesh is made, Solver should be 3-order, isn’t it?

If the order of mesh and solve don’t match, will be “ValueError: The multigrid level orders cannot exceed the solution order” error message shown?

I want to confirm it because I was stuggling this error message.

That error is coming from p-multigrid, this isn’t to do with the mesh.

P-multigrid is a convergence accelleration method that can be simply decsribed as down sampling the solution to a lower order, solving that system, and using that low order data to help solve the high order system.

Try changing the solver-dual-time-integrator-multip section in the ini the following for order 1:

[solver-dual-time-integrator-multip]
pseudo-dt-fact = 1.75
cycle = [(1, 1), (0, 4), (1, 4)]

1 Like

When I use 2nd-order ac-NS, could you tell me the best pseudo-dt-fact and cylcle of solver-dual-time-integrator-multip?

This is very often problem specific. You will need to do some experimentation here to find out what works best for your problem.

In general I would start with a V cycle and go from there.

Regards, Freddie.

My case is working well at 1st-order of [solver-dual-time-integrator-multip] pseudo-dt-fact = 1.75 cycle = [(1, 1), (0, 4), (1, 4)] I am not sure what works best for my problem at 2nd and 3rd order.
I would try something like the following a go for order=2:
[(2,1),(1,1),(0,5),(1,2),(2,4)]

But you could probably go straight to order=3 after your order=1 run.