Incompressible BCs not fully enforced

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.

Thank your reply.
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.
Moreover I am not familiar with V cycle. What is V cycle?

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.

1 Like

I noticed that I cannot understand pressure out-let condition.
In pyfr Pressure, what paramiters are used to do no-dimensional?
I think there are 2 types:(upper- letters are no-dimensional,lower-case letters are dimensional.
No1. is p statics pressure in P=p/rho * u * u ?
No2. is p dynamic pressure in P=p/rho * u * u?
If No2 is applyed, I think reference P is 1. However No1 is applyed, refecenre P must be calculated by the pressure and dencity of the atmosphere

Pyfr does no explicit non-dimensionalisation, this is all left to the user in how they define the initial and boundary conditions.

1 Like

Thank you replying.
Sorry I couldn’t understand exactly .
If reference velocity 10m/s ,reference deniers 1kg/m3 and pressure 1000hPa, what pressure should be put on outlet boundary?
I define velocity 10 m/s = 1 and density = 1kg/m3 as non-dimensionaly.
And in terms of outlet condition, is it Neumann boundary?

At the outlet you can put whatever pressure you like at the outlet boundary, but this is normally depermined by the case you are investigating.

In terms of is it a Neumann oundary condition, no, all the conditions are Dirichlet boundary conditions.

Thank you.

I understood it. It is because Pressure “GRADIENT” is used in NS equations, it isnt?
Best

Sorry again.My simulation is working with 1st order and low Re. But I have a problem at this oder and Re before calculation with 3rd order and high Re.
At Inlet condition, the pressure at the edge of wall is getting rasing. Why dose this pressure raising occurr? Is it possible to fix it? I apply the power low model (Uin = U0 * z ^(1/7),at z= 0 Uin =0.0) . So, near wall velocity close to 0.

 [backend]
precision = double
rank-allocator = linear

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

[constants]
nu = 1.48e-3
;nu = 0.1
Uin = 1.0
Vin = 0.0
Win = 0.0
Pc = 783.642
ac-zeta = 2.5
rhom = 1.0
;gas constance
R = 287.1
;input flow n
n = 7.0
;input hight 
bt = 7.5
[solver]
system = ac-navier-stokes
order = 1

[solver-time-integrator]
formulation = dual
scheme = sdirk33
pseudo-scheme = rk45
controller = none
pseudo-controller = local-pi
tstart = 0.0
tend = 80.0
dt = 0.001
pseudo-dt = 0.0001
pseudo-niters-min = 3
pseudo-niters-max = 10
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 = [(1, 1), (0, 4), (1, 4)]

[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

[solver-elements-quad]
soln-pts = gauss-legendre

[solver-interfaces-tri]
flux-pts = williams-shunn
quad-deg = 10
quad-pts = williams-shunn

[soln-plugin-nancheck]
nsteps = 10

[soln-plugin-pseudostats]
flushsteps = 100
file = ./csv/pseudostats_results1.csv
header = true

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

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



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

avg-p = p
;avg-rho = rho
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)
;fun-avg-p = p
;fun-avg-rho = rho

[soln-plugin-integrate]
nsteps = 100
file = ./csv/integral_results1.csv
header = true
quad-deg = 9
div1 = grad_u_x
div2 = grad_v_y
div3 = grad_w_z
vor1 = (grad_w_y - grad_v_z)
vor2 = (grad_u_z - grad_w_x)
vor3 = (grad_v_x - grad_u_y)

;int-E = rho*(u*u + v*v + w*w)
;int-enst = rho*(%(vor1)s*%(vor1)s + %(vor2)s*%(vor2)s + %(vor3)s*%(vor3)s)
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

velocity

So in answer to your question about Neumann boundary conditions, for conventional incompressible Navier–Stokes simulation you would probably apply a pressure BC via a Neumann condition. However, PyFR uses artifical compressibility and as a result it is more natural to apply a pressure BC as a Dirichlet condition.

As to your inlet, I think the pressure you are seeing is correct. At the inlet you are applying a BC of constant velocity and then adjacent to this there is a wall boundary condition. To enforce this condition the flow will need to decelerate and so a pressure gradient arises. You could get around this by changing you inlet velocity profile, but this depends on what you are trying to simulate.

Thank you. I looked up the velocity U at z = 0 at Inlet gets approx 0.4. So it is the reason why pressure is increasing at z= 0 at Inlet.
I set that “U =0” at z=0 at Inlet using u = Uin * pow((z / bt), (1.0/n)). But something error 0.4 is remaning. Why?

I’m not completely sure, but I would guess that the profile you’re imposing doesn’t quite match the boundary layer profile for this case. And so the flow adjusts which causes a pressure spike. I think this is probably fine.

Thank you for replying. I undestood.
I did 3 test cases which is for inlet.
Case No1 at [soln-bcs-inlet] and [soln-ics] Inlet u is u = Uin * pow((z / 7.5), (1.0/8.0))
Case No2 at [soln-bcs-inlet] and [soln-ics] Inlet u is u = Uin * sqrt((z / 7.5))* sqrt((z / 7.5))* sqrt((z / 7.5))sqrt((z/7.5))
Case No3 at [soln-bcs-inlet] and [soln-ics] Inlet u is u = Uin * (z / 7.5)
(z / 7.5)

Case No2 and Case No3 of [soln-bcs-inlet] and [soln-ics] must be same.

No1.test case1
I added

[soln-plugin-sampler]
nsteps = 1
samp-pts = [(0.0, 0.0, 0.0), (0.0, 0.0, 0.01),(0.0, 0.0, 15.0)]
format = primitive
file = point-data.csv
header = true

in a ini file and I cheked (0.0, 0.0, 0.0) which is the point at z= 0 of inlet (It is suporsed to be U = 0 at z=0) .
However the U data at (0.0, 0.0, 0.0) is 0.518395365 at z=0.
I attach the No1 test case at pyfr_test1 - Google Drive

No2. test case 2
the U data at (0.0, 0.0, 0.0) is -0.0003160222634275517.
I attach No2 test case at pyfr_test2 - Google Drive

No3. test case3
the U data at (0.0, 0.0, 0.0) is -0.00031602226342755175.
test case3 matches test case2. However u is negative.
I attach No3 test case at pyfr_test3 - Google Drive

I do not know why this is.
Could you look up my test cases if possible?

Is the mesh z extent what you think it should be? You can see this in paraview

So sorry, the [soln-plugin-sampler] should be at [(-20.0, 0.0, 0.0) which is on the Inlet surface at z=0.
However I got about inlet u =0.7 at z=0.
Your advise " I would guess that the profile you’re imposing doesn’t quite match the boundary layer profile for this case." is so useful that I think correct Inlet velocity. But I am not sure which is correct… I would like to apply figure expression.
image

I attach files at pyfr_test4 - Google Drive

and configulation file is

[backend]
precision = double
rank-allocator = linear

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

[constants]
nu = 1.48e-3
#nu = 0.1
Uin = 1.0e0
Vin = 0.0
Win = 0.0
Pc = 783.642
ac-zeta = 2.5
rhom = 1.0

R = 287.1

n = 7.0

bt = 1.5
[solver]
system = ac-navier-stokes
order = 1

[solver-time-integrator]
formulation = dual
scheme = sdirk33
pseudo-scheme = rk45
controller = none
pseudo-controller = local-pi
tstart = 0.0
tend = 0.01
dt = 0.01
pseudo-dt = 0.0001
pseudo-niters-min = 3
pseudo-niters-max = 10
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 = [(1, 1), (0, 4), (1, 4)]

[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

[solver-elements-quad]
soln-pts = gauss-legendre

[solver-interfaces-tri]
flux-pts = williams-shunn
quad-deg = 10
quad-pts = williams-shunn

[soln-plugin-nancheck]
nsteps = 1

[soln-plugin-pseudostats]
flushsteps = 1
file = ./csv_test/pseudostats_results1.csv
header = true

[soln-plugin-writer]
dt-out = 0.01
basedir = ./results_test
basename = target_with_fluid_results1-{t:.2f}

[soln-plugin-fluidforce-bwall]
nsteps = 1
file = ./csv_test/bluding-forces_results1.csv
header = true
morigin = (0.0, 0.0, 0.0)

[soln-plugin-tavg]
nsteps = 1
dt-out = 0.01
mode = continuous
basedir = ./results_test
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 = 1
file = ./csv_test/integral_results1.csv
header = true
quad-deg = 9
div1 = grad_u_x
div2 = grad_v_y
div3 = grad_w_z
vor1 = (grad_w_y - grad_v_z)
vor2 = (grad_u_z - grad_w_x)
vor3 = (grad_v_x - grad_u_y)

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-plugin-sampler]
nsteps = 1
samp-pts = [(-20.0, 0, 0.0), (-20.0, 0, 0.01),(-20.0, 5.0, 15.0)]
format = primitive
file = ./csv_test/point-data.csv
header = true

[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 = 0
u = Uin *pow(z,1.0/n)
v = Vin
w = Win

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