Skip to content

Commit efba4d4

Browse files
ranochaKostaszki
ranocha
authored and
Kostaszki
committed
fix TGV test case; add computation of vorticity for Euler eqs.; closes #35 
1 parent 310c15f commit efba4d4

4 files changed

Lines changed: 150 additions & 1 deletion

File tree

‎examples/ideal_gas_Euler.m‎

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -69,17 +69,26 @@
6969
end
7070
%% Plot numerical solution
7171
num_nodes = I_Mesh('NODES_X')*I_Mesh('NODES_Y')*I_Mesh('NODES_Z');
72+
field_u2(:) = field_u1;
73+
cl_run_kernel(I_Tech('device'), 'compute_vorticity', I_BalanceLaws('g_range'), I_BalanceLaws('l_range'), field_u1, field_u2, 0);
7274
if strcmp(I_Tech('memory_layout'), 'USE_ARRAY_OF_STRUCTURES')
7375
field_u1_plot = reshape(field_u1(1:num_nodes*I_BalanceLaws('NUM_TOTAL_VARS')), [I_BalanceLaws('NUM_TOTAL_VARS'), num_nodes]);
76+
field_u2_plot = reshape(field_u2(1:num_nodes*I_BalanceLaws('NUM_TOTAL_VARS')), [I_BalanceLaws('NUM_TOTAL_VARS'), num_nodes]);
7477
elseif strcmp(I_Tech('memory_layout'), 'USE_STRUCTURE_OF_ARRAYS')
7578
field_u1_tmp = reshape(field_u1, I_Tech('NUM_NODES_PAD'), I_BalanceLaws('NUM_TOTAL_VARS'));
7679
field_u1_plot = field_u1_tmp(1:num_nodes, :)';
80+
field_u2_tmp = reshape(field_u2, I_Tech('NUM_NODES_PAD'), I_BalanceLaws('NUM_TOTAL_VARS'));
81+
field_u2_plot = field_u2_tmp(1:num_nodes, :)';
7782
else
7883
error('You must USE_ARRAY_OF_STRUCTURES or USE_STRUCTURE_OF_ARRAYS.')
7984
end
8085

8186
%Optional plots
8287
if ismember(lower(char(I_RunOps('plot_numerical_solution'))),{'x','y','z','xy', 'xz', 'yz', 'xyz'})
8388
plot_2D(field_u1_plot, I_RunOps('plot_numerical_solution'),...
84-
I_Mesh('NODES_X'), I_Mesh('NODES_Y'), I_Mesh('NODES_Z'), 'Numerical Solution', 1, 1);
89+
I_Mesh('NODES_X'), I_Mesh('NODES_Y'), I_Mesh('NODES_Z'), 'Numerical Solution', 5, 5);
90+
91+
plot_2D(field_u2_plot, I_RunOps('plot_numerical_solution'),...
92+
I_Mesh('NODES_X'), I_Mesh('NODES_Y'), I_Mesh('NODES_Z'), 'Vorticity', 6, 8);
8593
end
94+

‎examples/ideal_gas_Euler_Taylor_Green_vortex.m‎

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
%clc
2+
clear
3+
close all
4+
5+
addpath('../matlab')
6+
7+
BalanceLaws.prepare_vars();
8+
global I_Mesh I_TI I_BalanceLaws I_Tech I_RunOps I_Results
9+
10+
N = uint32(64);
11+
I_Mesh('NODES_X') = N; I_Mesh('NODES_Y') = N; I_Mesh('NODES_Z') = N;
12+
I_Mesh('XMIN') = 0.0; I_Mesh('XMAX') = 6.283185307179586; % = 2*pi
13+
I_Mesh('YMIN') = 0.0; I_Mesh('YMAX') = 6.283185307179586;
14+
I_Mesh('ZMIN') = 0.0; I_Mesh('ZMAX') = 6.283185307179586;
15+
16+
I_TI('final_time') = 10;
17+
I_TI('cfl') = 0.85;
18+
19+
dt = I_TI('cfl') * (I_Mesh('YMAX') / double(I_Mesh('NODES_Y')-1)) / 10;
20+
num_steps = ceil(I_TI('final_time')/dt);
21+
dt = I_TI('final_time') / num_steps;
22+
23+
%Time integrator: SSPRK33 SSPRK104 CalvoFrancoRandez2R64
24+
%KennedyCarpenterLewis2R54C CarpenterKennedy2N54 ToulorgeDesmet2N84F
25+
I_TI('time_integrator') = 'CarpenterKennedy2N54';
26+
27+
I_TI('DT') = dt;
28+
I_TI('num_steps') = num_steps;
29+
30+
I_Tech('device') = 1;
31+
I_Tech('REAL') = 'double'; % float, double
32+
I_Tech('REAL4') = sprintf('%s4',I_Tech('REAL')); %Vector datatype
33+
I_Tech('memory_layout') = 'USE_ARRAY_OF_STRUCTURES'; % USE_ARRAY_OF_STRUCTURES, USE_STRUCTURE_OF_ARRAYS
34+
35+
% Use as kernel defines to keep consistency with header files?
36+
I_BalanceLaws('NUM_CONSERVED_VARS') = 5;
37+
I_BalanceLaws('NUM_AUXILIARY_VARS') = 4;
38+
I_BalanceLaws('NUM_TOTAL_VARS') = I_BalanceLaws('NUM_CONSERVED_VARS') + I_BalanceLaws('NUM_AUXILIARY_VARS');
39+
40+
%Compiler based optimizations
41+
if strcmp(I_Tech('REAL'),'float')
42+
I_Tech('optimizations') = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only -cl-single-precision-constant';
43+
else
44+
I_Tech('optimizations') = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only';
45+
end
46+
47+
I_RunOps('periodic') = 'USE_PERIODIC'; % 'USE_PERIODIC', 'USE_PERIODIC'; must be set to 'USE_PERIODIC'
48+
% if periodic boundary conditions should be used
49+
50+
I_RunOps('order') = 6; I_RunOps('operator_form') = 'classical'; % order: 2, 4, 6; operator_form: classical, extended
51+
I_RunOps('conservation_laws') = 'ideal_gas_Euler';
52+
I_RunOps('testcase') = 'Taylor_Green_vortex';
53+
I_RunOps('plot_numerical_solution') = 'z';
54+
I_RunOps('save_integrals_over_time') = false;
55+
% Choose between L2 and LInfinity norm for error calculation
56+
I_RunOps('norm') = 'LInf'; % L2, LInf
57+
%% Initialize variables
58+
[field_u1, field_u2] = BalanceLaws.initialize();
59+
60+
fprintf('Testcase: %s \nOrder: %d \nTime integrator: %s\nDT: %.16e N_STEPS: %5d FINAL_TIME: %.16e\nDX: %.16e NODES_X: %5d\nDY: %.16e NODES_Y: %5d\nDZ: %.16e NODES_Z: %5d \nREAL: %s\n\n',...
61+
I_RunOps('testcase'), I_RunOps('order'), I_TI('time_integrator'), I_TI('DT'), I_TI('num_steps'), I_TI('final_time'), I_Mesh('DX'), I_Mesh('NODES_X'), I_Mesh('DY'), I_Mesh('NODES_Y'), I_Mesh('DZ'), I_Mesh('NODES_Z'), I_Tech('REAL'));
62+
%% Compute numerical solution
63+
BalanceLaws.compute_numerical_solution(field_u1, field_u2);
64+
fprintf('Total runtime: %.3f seconds Kernel runtime: %d\n', I_Results('runtime'), I_Results('kernel_runtime'));
65+
66+
rel_err = I_Results('rel_err');
67+
for comp=0:I_BalanceLaws('NUM_CONSERVED_VARS') - 1
68+
fprintf('Relative Error of Field Component %d: %.15f %%\n', comp, 100*rel_err(comp + 1))
69+
end
70+
71+
%% Plot numerical solution
72+
num_nodes = I_Mesh('NODES_X')*I_Mesh('NODES_Y')*I_Mesh('NODES_Z');
73+
field_u2(:) = field_u1;
74+
cl_run_kernel(I_Tech('device'), 'compute_vorticity', I_BalanceLaws('g_range'), I_BalanceLaws('l_range'), field_u1, field_u2, 0);
75+
if strcmp(I_Tech('memory_layout'), 'USE_ARRAY_OF_STRUCTURES')
76+
field_u1_plot = reshape(field_u1(1:num_nodes*I_BalanceLaws('NUM_TOTAL_VARS')), [I_BalanceLaws('NUM_TOTAL_VARS'), num_nodes]);
77+
field_u2_plot = reshape(field_u2(1:num_nodes*I_BalanceLaws('NUM_TOTAL_VARS')), [I_BalanceLaws('NUM_TOTAL_VARS'), num_nodes]);
78+
elseif strcmp(I_Tech('memory_layout'), 'USE_STRUCTURE_OF_ARRAYS')
79+
field_u1_tmp = reshape(field_u1, I_Tech('NUM_NODES_PAD'), I_BalanceLaws('NUM_TOTAL_VARS'));
80+
field_u1_plot = field_u1_tmp(1:num_nodes, :)';
81+
field_u2_tmp = reshape(field_u2, I_Tech('NUM_NODES_PAD'), I_BalanceLaws('NUM_TOTAL_VARS'));
82+
field_u2_plot = field_u2_tmp(1:num_nodes, :)';
83+
else
84+
error('You must USE_ARRAY_OF_STRUCTURES or USE_STRUCTURE_OF_ARRAYS.')
85+
end
86+
87+
%Optional plots
88+
if ismember(lower(char(I_RunOps('plot_numerical_solution'))),{'x','y','z','xy', 'xz', 'yz', 'xyz'})
89+
plot_2D(field_u1_plot, I_RunOps('plot_numerical_solution'),...
90+
I_Mesh('NODES_X'), I_Mesh('NODES_Y'), I_Mesh('NODES_Z'), 'Numerical Solution', 5, 5);
91+
92+
plot_2D(field_u2_plot, I_RunOps('plot_numerical_solution'),...
93+
I_Mesh('NODES_X'), I_Mesh('NODES_Y'), I_Mesh('NODES_Z'), 'Vorticity', 6, 8);
94+
end
95+

‎include_physics/ideal_gas_Euler.h‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1933,4 +1933,24 @@ inline void analytical_solution(uint ix, uint iy, uint iz, global REAL *u, REAL
19331933
set_field_component(ix, iy, iz, Field_E, u, E_ana);
19341934
}
19351935

1936+
1937+
1938+
1939+
/*
1940+
Compute the vorticity of the velocity in u1 and store it as the velocity of u2.
1941+
*/
1942+
kernel void compute_vorticity(global REAL const *u1, global REAL *u2) {
1943+
1944+
uint ix = get_global_id(0);
1945+
uint iy = get_global_id(1);
1946+
uint iz = get_global_id(2);
1947+
1948+
REAL4 vorticity = curl(ix, iy, iz, u1, Field_ux);
1949+
set_field_component(ix, iy, iz, Field_ux, u2, vorticity.x);
1950+
set_field_component(ix, iy, iz, Field_ux, u2, vorticity.y);
1951+
set_field_component(ix, iy, iz, Field_ux, u2, vorticity.z);
1952+
1953+
return;
1954+
}
1955+
19361956
#endif // IDEAL_GAS_EULER

‎include_testcases/ideal_gas_Euler_Taylor_Green_vortex.h‎

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,3 +47,28 @@ inline REAL p_init(uint ix, uint iy, uint iz) {
4747

4848
return (REAL)(100.0/GAMMA + 0.0625 * (cos(2*x)*cos(2*z) + 2*cos(2*y) + 2*cos(2*x) + cos(2*y)*cos(2*z)));
4949
}
50+
51+
52+
/*
53+
There is no analytical solution of the density.
54+
*/
55+
inline REAL rho_analytical(uint ix, uint iy, uint iz, REAL time) {
56+
57+
return rho_init(ix, iy, iz);
58+
}
59+
60+
/*
61+
There is no analytical solution of the velocity.
62+
*/
63+
inline REAL4 u_analytical(uint ix, uint iy, uint iz, REAL time) {
64+
65+
return u_init(ix, iy, iz);
66+
}
67+
68+
/*
69+
There is no analytical solution of the pressure.
70+
*/
71+
inline REAL p_analytical(uint ix, uint iy, uint iz, REAL time) {
72+
73+
return p_init(ix, iy, iz);
74+
}

0 commit comments

Comments
 (0)