Suggest to me for calculation channel flow, please

[Natky]

Hi,
I'm trying to test channel flow case on HiFiLES-code, when I calculated periodic channel flow 3D for lamina flow follow set up my input file, eles.cpp and my domain , at calculation round massflux = -60.7985xxx , Ubulk = -9.67758xxx and bodyforce(1) = -1, the code have show “Fatal error”
There is a Fatal error:

Fatal error 'NaN in residual, exiting.' at ../src/eles.cpp:1638

For set up and edit “eles.cpp” at line

6401 // case-specific parameters
6402 // periodic channel:
6403 area = 2_pi;
6404 vol = 4_pi*pi;

6406 // periodic hill (HIOCFD3 Case 3.4):
6407 //area = 9.162;
6408 //vol = 114.34;
6409 //mdot0 = 9.126; initial mass flux ********** i dont know that 's set up comment "mdot0" , you know? or suggest to me, ****************

and

6438 // set body force for streamwise momentum and energy
6439 body_force(0) = 0.;
6440 //body_force(1) = alpha/area/dt_(mdot0 - mass_flux); // modified SD3D version
6441 //body_force(1) = 1.0/area/dt_(mdot0 - 2.0_mass_flux + mdot_old); // HIOCFD C3.4 version
6442 body_force(1) = -1.0; //Hack for bodyforce constant
6443 body_force(2) = 0.;
6444 body_force(3) = 0.;
6445 body_force(4) = body_force(1)_ubulk; // energy forcing

Here is my set up domain periodic channel flow and grid element
Lx = 0 : 2*pi
Ly = 0 : pi
Lz = -1 : 1
and cell element : nonuniform grid = ( 24 * 12 * 28 )

and here is my input file

---

Solver parameters

// 0: Euler/Navier-Stokes, 1:Advection/Adv-Diffusion
equation 0
viscous 1

// 0: Rusanov, 1: Lax-Friedrich, 2: Roe
riemann_solve_type 0
vis_riemann_solve_type 0

// 0: Isentropic Vortex, 1: Uniform flow, 2: Sine Wave
ic_form 6

// 0: Normal (doesn't have an analytical solution), 1:Isentropic Vortex, 2: Advection-Equation
test_case 0
order 5 // Order of basis polynomials
dt_type 0 // 0: User-supplied, 1: Global, 2: Local
dt 0.0001
CFL 1
n_steps 1000000
adv_type 3 // 0: Forward Euler, 3: RK45
tau 1.0
pen_fact 0.5

---

Mesh options

mesh_file ch3d_pi.msh
dx_cyclic 6.283185307179586
dy_cyclic 3.141592653589793
dz_cyclic 2.0

---

Monitoring, plotting parameters

p_res 4 // Plotting resolution, # of nodes per direction
write_type 0 // 0: Paraview, 1: Tecplot
// Choose extra fields to be written to file: u v w energy pressure mach vorticity q_criterion.
// Set to 0 or comment out for no diagnostic fields
n_diagnostic_fields 5 u v w pressure mach
inters_cub_order 4 // Order of cubature rule for integrating over element interfaces
volume_cub_order 4 // Order of cubature rule for integrating over element volumes

monitor_res_freq 100
plot_freq 500
data_file_name ch3d_v01zy_
restart_dump_freq 10000
res_norm_type 2 // 1:L1 norm, 2:L2 norm, 0:infinity norm
error_norm_type 2 // 1:L1 norm, 2:L2 norm, 0:infinity norm
res_norm_field 0 // Density

---

Element parameters

==== Tris ====
upts_type_tri 0
fpts_type_tri 1
vcjh_scheme_tri 3
c_tri 0.0
sparse_tri 0

==== Quads ====
upts_type_quad 0 // 0: Gauss, 1: Gauss-Lobatto
vcjh_scheme_quad 1 // 0: VCJH, 1: DG, 2: SD, 3: Hu, 4: c_+
eta_quad 0.0
sparse_quad 0

==== Hexs ====
upts_type_hexa 0
vcjh_scheme_hexa 1
eta_hexa 0.
sparse_hexa 0

==== Tets ====
upts_type_tet 1
fpts_type_tet 0
vcjh_scheme_tet 1
eta_tet 0.0
sparse_tet 0

==== Prisms ====
upts_type_pri_tri 0
upts_type_pri_1d 0
vcjh_scheme_pri_1d 1
eta_pri 0.0
sparse_pri 0

---

Fluid Parameters

gamma 1.4
prandtl 0.72
S_gas 120.
T_gas 291.15
R_gas 286.9
mu_gas 1.827E-05

---

Boundary conditions

==== Viscous ====
fix_vis 0 // 0: Sutherland's law, 1: Constant viscosity
Mach_free_stream 0.2
nx_free_stream 1.
ny_free_stream 0.
nz_free_stream 0.
Re_free_stream 180.
L_free_stream 1.
T_free_stream 300
T_wall 300
//Mach_wall 0.
//nx_wall 1.
//ny_wall 0.
//nz_wall 0.
//T_wall 300

==== Inviscid ====
rho_bound 1.
u_bound 1.0
v_bound 1.0
w_bound 0.0
p_bound 17.857142857142854098

---

Forcing

body_forcing 1 // Used for the periodic channel test case. 0 = no forcing, 1 = forcing

Initial Profile: If ic_form=6, use these parameters to describe polynomial initial condition. Useful for cyclic problems. Specify coeffs of poly of form a + bx + cx^2 + dx^3 + ex^4 + by + cy^2 + dy^3 + ey^4 + bz + cz^2 + dz^3 + ez^4

// Coeffs of x-momentum [need the 13 so reader know how many quantities]
x_coeffs 13 1.875 0 0 0 0 0 0 0 0 0 -3.75 0 1.875
// Coeffs of y-momentum
y_coeffs 13 0 0 0 0 0 0 0 0 0 0 0 0 0
// Coeffs of y-momentum
z_coeffs 13 0 0 0 0 0 0 0 0 0 0 0 0 0
perturb_ic 1 // If 1, add velocity perturbation to kickstart turbulence in periodic channel. 0 otherwise.

Please, suggest to me
thank you,

suggest to me or contact to me : [email protected]

[mlopez14]

The timestep (1e-4) seems too large for the order (5) of the scheme. I would suggest you try a lower order and a smaller timestep. Let us know if this error happens at order = 1 and dt less than 1e-6.

[Natky]

Thank you for your suggestion and I have set timestep = 1e-6 and set to use the order = 1 of the scheme, but The code can't calculate periodic channel flow. Could you please check the setting in eles.cpp and my input file for me again?

And there is a fatal error:

---

iter, mdot0, mdot_old, mass_flux, body_force(1): 4538, 0.00000000, -60.39439814, -60.42285522, -1.00000000
iter, mdot0, mdot_old, mass_flux, body_force(1): 4539, 0.00000000, -60.42285522, -60.45131858, -1.00000000
iter, mdot0, mdot_old, mass_flux, body_force(1): 4540, 0.00000000, -60.45131858, -60.47978820, -1.00000000

Fatal error 'NaN in residual, exiting.' at ../src/eles.cpp:1638

---

[mlopez14]

Correct me if I'm wrong, but it seems the flow is being initialized going in the +x direction (nx_free_stream 1.) and the forcing function is returning a negative value, so the forcing is going in the -x direction. I would try initializing the flow to go in the -x direction (nx_free_stream -1.).

[Natky]

Thank you for suggestion. But i found the problem during code calculation same as before. I think the problem causes from uncomplete set up at 'eles.cpp' file.Could you please give me a suggestion about setting up the code for calculation periodic channel flow and then i will use the code to develop the application for more calulation on periodic turbulent channel flow domain. if you have any manual or the code user guide, could you please let me know or give it to me?

Thank you very much in advance.

[mlopez14]

The code user guide is the wiki we have here in github. The forcing function is experimental and was used only for one instance that is not published yet. In what part of the domain is the code blowing up? Is it close to the boundary, right at the middle? As long as the boundary is labeled periodic in the mesh and the correct dx_cyclic, etc. variables are in the input file (they seem to be correct), the code would let you know if there is a mismatch.

Do the elements on the periodic boundaries match as if they were adjacent? If they do not, the code would be blowing up at the boundaries.

[Natky]

Thank you for suggestion.