Hi all,
I was able to make simple mesh myself. However, it’s result is not right…
Could you teach me how to fix it?
What I want to do is that the external fluid calculation around 3D cude.
I attach the image fig and the result I have got bellow.
The files of mesh, geo, ini and so on are attached in mesh1_test.zip - Google Drive
Thank you for your always helping.
Best,
Yuji
From a cursory look at the msh file your mesh seems reasonable, although I’m not as familar with the gmsh v4 format as v2.
Some things you might want to consider are:
abs(grad_u_x + grad_v_y) ie the integral of the absolute value of the divergence. This should be identically zero for incompressible flows so if it’sparticularly large then the system isn’t converging. (Note that your impulsive IC and BCs will probably cause some non-zero divergnce initially as the BC get enforced)What is the figure supposed to represent? Is it a flow snapshot after some number of time steps, or just the initial condition? I’ll note that the scale is extremely misleading, as the velocity is basically uniform throughout the entire domain despite what one may think from the figure.
Regards, Freddie.
Dear Dr. Will, Prof Freddie,
Thank your for replying.
The figure and below figure is the initial condition. I put the results at google drive mesh1_test.zip - Google Drive
I have got bellow error messages I cannot fix…
$pyfr run -b cuda -p tartget_with_fluid2.pyfrm tartget_with_fluid2.ini
0.7% [> ] 0.49/75.00 ela: 00:00:11 rem: 00:28:36Traceback (most recent call last):
File "/home/ubuntu2204/libxsmm/pyfr-venv/bin/pyfr", line 8, in <module>
sys.exit(main())
File "/home/ubuntu2204/libxsmm/pyfr-venv/lib/python3.10/site-packages/pyfr/__main__.py", line 118, in main
args.process(args)
File "/home/ubuntu2204/libxsmm/pyfr-venv/lib/python3.10/site-packages/pyfr/__main__.py", line 251, in process_run
_process_common(
File "/home/ubuntu2204/libxsmm/pyfr-venv/lib/python3.10/site-packages/pyfr/__main__.py", line 247, in _process_common
solver.run()
File "/home/ubuntu2204/libxsmm/pyfr-venv/lib/python3.10/site-packages/pyfr/integrators/base.py", line 115, in run
self.advance_to(t)
File "/home/ubuntu2204/libxsmm/pyfr-venv/lib/python3.10/site-packages/pyfr/integrators/dual/phys/controllers.py", line 56, in advance_to
self._accept_step(self.pseudointegrator._idxcurr)
File "/home/ubuntu2204/libxsmm/pyfr-venv/lib/python3.10/site-packages/pyfr/integrators/dual/phys/controllers.py", line 35, in _accept_step
csh(self)
File "/home/ubuntu2204/libxsmm/pyfr-venv/lib/python3.10/site-packages/pyfr/plugins/nancheck.py", line 21, in __call__
raise RuntimeError(f'NaNs detected at t = {intg.tcurr}')
RuntimeError: NaNs detected at t = 0.5000000000000002
Ok so the boundary conditions won’t be enforced on the initial condition as the IC will be exactly what you defined in the ini file. On the first time step they’ll be enforced, but asa with all simulations it will take some time for these to develop.
This error message means that your simulation has diverged, there are a number of things that can cause this, but normally its just that your time step is too big. Have a look at the things I listed above, the integral and the pseudostats can be helpfully in working out what needs to be done.
Thank you for replying.
To start with the ini modify, I changed velocity and presure inlet and outlet condition, and the max number of pseudo iterations (I attach ini file). The initial velocity seems correct, but puresure is as same as before velocity fig.(I attach the preview fig) . I cannot solve it still now.
Now I am looking for the way to set up the integration plugin and calculating abs(grad_u_x + grad_v_y) .
[backend]
precision = double
[constants]
nu = 0.005
Uin = 0.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 = 75.0
dt = 0.01
pseudo-dt = 0.005
pseudo-niters-min = 3
pseudo-niters-max = 3
pseudo-resid-norm = l2
pseudo-resid-tol = 1e-12
atol = 1e-6
pseudo-dt-max-mult = 2.5
[solver-dual-time-integrator-multip]
pseudo-dt-fact = 50
cycle = [(3, 1), (2, 1), (1, 1), (0, 2), (1, 1), (2, 1), (3, 3)]
[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 = 50
[soln-plugin-pseudostats]
flushsteps = 20
file = residual.csv
header = true
[soln-plugin-writer]
dt-out = 5.0
basedir = .
basename = tartget_with_fluid2-{t:.2f}
[soln-bcs-wall]
type = no-slp-wall
[soln-bcs-inlet]
type = ac-char-riem-inv
ac-zeta = 180
p = Pc
u = Uin
v = Vin
w = Win
[soln-bcs-inlete]
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 = 2.0
u = Uin
v = Vin
w = Win
[soln-ics]
u = Uin
v = Vin
w = Win
p = Pc
pressure fig
velocity fig
In the above you have:
Uin = 0.0
Vin = 0.0
Win = 0.0
Hence why the velocity is zero. Also If those screen shots are of the solution at t=0, that is perfectly normal.
I appliciate your comments.
I have checked initial condition is applied correctly. However when boundary condition is applied, NAN is produced(grad u seems nearby 0(the relult is bellow)(I think it is correct.).
I wrote in ini file
[soln-plugin-integrate]
nsteps = 1
file = ./results/csv/integral.csv
header = true
quad-deg = 9
div1 = grad_u_x
div2 = grad_v_y
div3 = grad_w_z
int-dive = %(div1)s + %(div2)s + %(div3)s
I got result
t,int-dive
0.0,-1.7856327527709936e-11
Due to the fact, I tried 2 cases in order to find the resone why NAN is produced,
No. 1
The U=10 P=0 initial condition is applied in the domein. The inlet BC is U=10 in ac-in-fv and outlet BC is P=0 ac-out-fp.
The ini file and simulation file are in case1_INLET_OUTLET_WALL_BC - Google Drive
No. 2
The U=10 P=0 initial condition is applied in the domein. Outlet BC is P=0 ac-out-fp.(NO INLET BC)
The ini file and simulation file are in case2_onlyOUTLET_WALL_BC - Google Drive
Both cases are calculated initial condition, but produced NAN into the next step.
I cannot find why these NAN have been produced. In my opinion, I think BC is not working well.
Could you teach me how to fix it?
Thank you.
I’m lookig at your ini file now, there are quite a few things wrong with it. Somethings to bear in mind:
p to p-1, the pseudo-dt in increasing by a factor of 50. Which will easily make things unstable.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 fix it along your comments.
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
[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 = ./results1/csv/residual_merge2.csv
header = true
[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
header = true
quad-deg = 9
div1 = grad_u_x
div2 = grad_v_y
div3 = grad_w_z
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.
Thank you for replyig.
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
[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 = ./results3/csv/residual.csv
header = true
[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
#errors header = true
#int-div = (grad_u_x + grad_v_y + grad_w_z)
[soln-plugin-integrate]
nsteps = 1
file = ./results3/csv/integral.csv
header = true
quad-deg = 9
div1 = grad_u_x
div2 = grad_v_y
div3 = grad_w_z
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.
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
[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 = 100
[soln-plugin-pseudostats]
flushsteps = 100
file = ./csv/pseudostats_results5.csv
header = true
[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
header = true
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
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-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.
