diff --git a/examples/ideal_MHD_alfven_fluxo.m b/examples/ideal_MHD_alfven_fluxo.m new file mode 100644 index 0000000..525bad7 --- /dev/null +++ b/examples/ideal_MHD_alfven_fluxo.m @@ -0,0 +1,91 @@ +%clc +clear +close all + +addpath('../matlab') + +BalanceLaws.prepare_vars(); +global I_Mesh I_TI I_BalanceLaws I_Tech I_RunOps I_Results + +% Alfven wave testcase used in [Fluxo](https://github.com/project-fluxo/fluxo), domain [-1,1]^3. +N = uint32(40); +I_Mesh('NODES_X') = N; I_Mesh('NODES_Y') = N; I_Mesh('NODES_Z') = uint32(5); +I_Mesh('XMIN') = -1.0; I_Mesh('XMAX') = 1.0; +I_Mesh('YMIN') = -1.0; I_Mesh('YMAX') = 1.0; +I_Mesh('ZMIN') = -1.0; I_Mesh('ZMAX') = 1.0; + +I_TI('final_time') = 25; %125 +I_TI('cfl') = 0.4; + +dt = I_TI('cfl') * 2.0 / double(I_Mesh('NODES_Y')); +num_steps = ceil(I_TI('final_time')/dt); +dt = I_TI('final_time') / num_steps; + +%Time integrator: SSPRK33 SSPRK104 CalvoFrancoRandez2R64 +%KennedyCarpenterLewis2R54C CarpenterKennedy2N54 ToulorgeDesmet2N84F +I_TI('time_integrator') = 'CarpenterKennedy2N54'; + +I_TI('DT') = dt; +I_TI('num_steps') = num_steps; + +I_Tech('device') = 1; +I_Tech('REAL') = 'double'; % float, double +I_Tech('REAL4') = sprintf('%s4',I_Tech('REAL')); %Vector datatype +I_Tech('memory_layout') = 'USE_STRUCTURE_OF_ARRAYS'; % USE_ARRAY_OF_STRUCTURES, USE_STRUCTURE_OF_ARRAYS + +% Use as kernel defines to keep consistency with header files? +I_BalanceLaws('NUM_CONSERVED_VARS') = 8; +I_BalanceLaws('NUM_AUXILIARY_VARS') = 4; +I_BalanceLaws('NUM_TOTAL_VARS') = I_BalanceLaws('NUM_CONSERVED_VARS') + I_BalanceLaws('NUM_AUXILIARY_VARS'); + +%Compiler based optimizations +if strcmp(I_Tech('REAL'),'float') + I_Tech('optimizations') = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only -cl-single-precision-constant'; +else + I_Tech('optimizations') = ' -cl-mad-enable -cl-no-signed-zeros -cl-finite-math-only'; +end + +I_RunOps('periodic_x') = 'NONE'; % 'NONE', 'USE_PERIODIC_X'; must be set to 'USE_PERIODIC_X' + % if periodic boundary conditions in x-direction should be used +I_RunOps('periodic_y') = 'NONE'; % 'NONE', 'USE_PERIODIC_Y'; must be set to 'USE_PERIODIC_Y' + % if periodic boundary conditions in y-direction should be used +I_RunOps('periodic_z') = 'USE_PERIODIC_Z'; % 'NONE', 'USE_PERIODIC_Z'; must be set to 'USE_PERIODIC_Z' + % if periodic boundary conditions in z-direction should be used + +I_RunOps('order') = 6; I_RunOps('operator_form') = 'classical'; % order: 2, 4, 6; operator_form: classical, extended +I_RunOps('conservation_laws') = 'ideal_MHD'; +I_RunOps('testcase') = 'alfven_fluxo'; +I_RunOps('plot_numerical_solution') = 'z'; +I_RunOps('save_integrals_over_time') = false; +% Choose between L2 and LInfinity norm for error calculation +I_RunOps('norm') = 'LInf'; % L2, LInf +%% Initialize variables +[field_u1, field_u2] = BalanceLaws.initialize(); + +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',... + 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')); +%% Compute numerical solution +BalanceLaws.compute_numerical_solution(field_u1, field_u2); +fprintf('Total runtime: %.3f seconds Kernel runtime: %d\n', I_Results('runtime'), I_Results('kernel_runtime')); + +rel_err = I_Results('rel_err'); +for comp=0:I_BalanceLaws('NUM_CONSERVED_VARS') - 1 + fprintf('Relative Error of Field Component %d: %.15f %%\n', comp, 100*rel_err(comp + 1)) +end + +%% Plot numerical solution +num_nodes = I_Mesh('NODES_X')*I_Mesh('NODES_Y')*I_Mesh('NODES_Z'); +if strcmp(I_Tech('memory_layout'), 'USE_ARRAY_OF_STRUCTURES') + field_u1_plot = reshape(field_u1(1:num_nodes*I_BalanceLaws('NUM_TOTAL_VARS')), [I_BalanceLaws('NUM_TOTAL_VARS'), num_nodes]); +elseif strcmp(I_Tech('memory_layout'), 'USE_STRUCTURE_OF_ARRAYS') + field_u1_tmp = reshape(field_u1, I_Tech('NUM_NODES_PAD'), I_BalanceLaws('NUM_TOTAL_VARS')); + field_u1_plot = field_u1_tmp(1:num_nodes, :)'; +else + error('You must USE_ARRAY_OF_STRUCTURES or USE_STRUCTURE_OF_ARRAYS.') +end + +%Optional plots +if ismember(lower(char(I_RunOps('plot_numerical_solution'))),{'x','y','z','xy', 'xz', 'yz', 'xyz'}) + plot_2D(field_u1_plot, I_RunOps('plot_numerical_solution'),... + I_Mesh('NODES_X'), I_Mesh('NODES_Y'), I_Mesh('NODES_Z'), 'Numerical Solution', 6, 8); +end diff --git a/include_testcases/ideal_MHD_alfven_fluxo.h b/include_testcases/ideal_MHD_alfven_fluxo.h new file mode 100644 index 0000000..329b066 --- /dev/null +++ b/include_testcases/ideal_MHD_alfven_fluxo.h @@ -0,0 +1,110 @@ +// This project is licensed under the terms of the Creative Commons CC BY-NC-ND 4.0 license. + +//-------------------------------------------------------------------------------------------------- +// Functions for field initialisation +//-------------------------------------------------------------------------------------------------- + +/* +Alfven wave testcase used in [Fluxo](https://github.com/project-fluxo/fluxo), domain [-1,1]^3. +*/ + +// paramters of the Alfven wave test case +#define CONST_FREQUENCY 1 +REAL constant CONST_omega = 6.283185307179586 * CONST_FREQUENCY; +REAL constant CONST_r = 2; // lenght-variable = lenght of computational domain +REAL constant CONST_e = 0.2; // epsilon = 0.2 + +/* +Analytical solution of the velocity field. +*/ +inline REAL4 u_analytical(uint ix, uint iy, uint iz, REAL time) { + + REAL x = (REAL)XMIN + ix*(REAL)DX; + REAL y = (REAL)YMIN + iy*(REAL)DY; + REAL z = (REAL)ZMIN + iz*(REAL)DZ; + + REAL nx = rsqrt(CONST_r*CONST_r + 1); + REAL ny = CONST_r * rsqrt(CONST_r*CONST_r + 1); + REAL sqr = 1; + REAL Va = CONST_omega / (ny * sqr); + REAL phi_alv = CONST_omega / ny * (nx * (x - 0.5*CONST_r) + ny * (y - 0.5*CONST_r)) - Va * time; + + REAL u1 = - CONST_e * ny * cos(phi_alv); + REAL u2 = CONST_e * nx * cos(phi_alv); + REAL u3 = CONST_e * sin(phi_alv); + + return (REAL4) {u1, u2, u3, (REAL)(0)}; +} + +/* +Analytical solution of the magnetic field. +*/ +inline REAL4 b_analytical(uint ix, uint iy, uint iz, REAL time) { + + REAL x = (REAL)XMIN + ix*(REAL)DX; + REAL y = (REAL)YMIN + iy*(REAL)DY; + REAL z = (REAL)ZMIN + iz*(REAL)DZ; + + REAL nx = rsqrt(CONST_r*CONST_r + 1); + REAL ny = CONST_r * rsqrt(CONST_r*CONST_r + 1); + REAL sqr = 1; + REAL Va = CONST_omega / (ny * sqr); + REAL phi_alv = CONST_omega / ny * (nx * (x - 0.5*CONST_r) + ny * (y - 0.5*CONST_r)) - Va * time; + + REAL B1 = nx + CONST_e * ny * cos(phi_alv) * sqr; + REAL B2 = ny - CONST_e * nx * cos(phi_alv) * sqr; + REAL B3 = - CONST_e * sin(phi_alv) * sqr; + + return (REAL4) {B1, B2, B3, (REAL)(0)}; +} + +/* +Boundary condition of the magnetic field. +*/ +inline REAL4 b_boundary(uint ix, uint iy, uint iz, REAL time) { + + return b_analytical(ix, iy, iz, time); +} + +/* +Boundary condition of the velocity field. +*/ +inline REAL4 u_boundary(uint ix, uint iy, uint iz, REAL time) { + + return u_analytical(ix, iy, iz, time); +} + + +/* Boundary condition of the pressure and density */ + +inline REAL p_boundary(uint ix, uint iy, uint iz, REAL time) { + + return (REAL)(1); +} + +inline REAL rho_boundary(uint ix, uint iy, uint iz, REAL time) { + + return (REAL)(1); +} + + +/* initial conditions for the simulation */ + +inline REAL4 b_initial(uint ix, uint iy, uint iz, REAL time){ + + return b_analytical(ix, iy, iz, time); +} + +inline REAL4 u_initial(uint ix, uint iy, uint iz, REAL time){ + + return u_analytical(ix, iy, iz, time); +} + +inline REAL p_initial(uint idx, uint iy, uint iz, REAL time){ + return (REAL)(1); +} + +inline REAL rho_initial(uint idx, uint iy, uint iz, REAL time){ + + return (REAL)(1); +}