A Dipole as a new example (#31)

* boundary conditions for MHD

* applied requested changes

* used a Macro instead of pow(..., 2)

* Added speeds from Bouchut, Klingenberg and Waagan (2012)

* proper handling of misdefined preprocessor symbols

* clear distinction between a analytic soloution and the initial and boundary values, which is needed for a new example

* Made some structural changes for the computation of the surface terms

* changed the USE_VR_* to USE_HLL_SPEED_*, #ifndef USE_PERIODIC to #ifdef USE_BOUNDARY_FLUX_HLL, max() to fmax and added 'const' in front of constant values, as Ranocha recommended.

* incooperated Ranochas comment on the preprocessor symbols and removed a compiler warning

* A new testcase: A dipole which stays next to the domain and oscillations of the magentic field on the opposite domain border.

* boundary conditions for MHD

* proper handling of misdefined preprocessor symbols

* clear distinction between a analytic soloution and the initial and boundary values, which is needed for a new example

* A new testcase: A dipole which stays next to the domain and oscillations of the magentic field on the opposite domain border.

* rebase on master

* ideal_MHD_Dipole.m

Author: simonius

examples/ideal_MHD_Dipole.m

@@ -0,0 +1,86 @@
+clc
+clear
+close all
+
+addpath('../matlab')
+
+BalanceLaws.prepare_vars();
+global I_Mesh I_TI I_BalanceLaws I_Tech I_RunOps I_Results
+
+N = uint32(100);
+I_Mesh('NODES_X') = N; I_Mesh('NODES_Y') = N; I_Mesh('NODES_Z') = N;
+I_Mesh('XMIN') = 0.0; I_Mesh('XMAX') = 1.0;
+I_Mesh('YMIN') = 0.0; I_Mesh('YMAX') = 1.0;
+I_Mesh('ZMIN') = 0.0; I_Mesh('ZMAX') = 1.0;
+
+I_TI('final_time') = 1;
+I_TI('cfl') = 0.65;
+
+dt = I_TI('cfl') * 1.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') = 'NONE'; % 'NONE', 'USE_PERIODIC'; must be set to 'USE_PERIODIC'
+                                       % if periodic boundary conditions should be used
+
+I_RunOps('order') = 4; I_RunOps('operator_form') = 'classical'; % order: 2, 4, 6; operator_form: classical, extended
+I_RunOps('conservation_laws') = 'ideal_MHD';
+I_RunOps('testcase') = 'far_dipole';
+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

include_defines/ideal_MHD_defines.h

@@ -57,13 +57,15 @@ TODO: JANHUNEN_CENTRAL, BRACKBILL_BARNES_CENTRAL, cf.
 */
 #define USE_SOURCE_GODUNOV_CENTRAL
 
+
 // Boundary Fluxes: USE_BOUNDARY_FLUX_HLL
 #define USE_BOUNDARY_FLUX_HLL
 
 // Speeds: USE_HLL_SPEED_KUSANO USE_HLL_SPEED_BKW
 #define USE_HLL_SPEED_BKW
 
 
+
 enum Fields {
              // Density
              Field_rho,

include_physics/ideal_MHD.h

@@ -362,9 +362,9 @@ inline void compute_ext_num_flux_z(REAL const* uk, REAL const* um, REAL* ext_num
 inline void init_fields(uint ix, uint iy, uint iz, global REAL* u) {
 
   REAL rho = rho_initial(ix, iy, iz, (REAL)(0));
-  REAL p = p_initial(ix, iy, iz, (REAL)(0)); 
-  REAL4 uinit = u_initial(ix, iy, iz, (REAL)(0)); 
-  REAL4 binit = b_initial(ix, iy, iz, (REAL)(0)); 
+  REAL p = p_initial(ix, iy, iz, (REAL)(0));
+  REAL4 uinit = u_initial(ix, iy, iz, (REAL)(0));
+  REAL4 binit = b_initial(ix, iy, iz, (REAL)(0));
 
   REAL um[NUM_TOTAL_VARS] = {0};
 
@@ -387,6 +387,7 @@ inline void init_fields(uint ix, uint iy, uint iz, global REAL* u) {
 // Surface terms
 //--------------------------------------------------------------------------------------------------
 
+
 #ifndef USE_PERIODIC
 
 #define sq(x) ((x)*(x))
@@ -432,6 +433,7 @@ inline void calc_flux_h(const REAL *u, REAL *flux){
 }
 
 
+
 #if defined USE_BOUNDARY_FLUX_HLL
 
 inline void calc_num_flux(REAL al, REAL ar, REAL *ul, REAL *ur, REAL *fluxl, REAL *fluxr, REAL *num_flux){
@@ -443,7 +445,7 @@ inline void calc_num_flux(REAL al, REAL ar, REAL *ul, REAL *ur, REAL *fluxl, REA
   else if(0 < ar) {
     for (i = 0; i < NUM_CONSERVED_VARS; i++)
       num_flux[i] = ((ar*fluxl[i]-al*fluxr[i]) + ar*al*(ur[i]-ul[i]))/(ar-al);
-  } 
+  }
   else {
     for (i = 0; i < NUM_CONSERVED_VARS; i++)
       num_flux[i] = fluxr[i];
@@ -453,7 +455,7 @@ inline void calc_num_flux(REAL al, REAL ar, REAL *ul, REAL *ur, REAL *fluxl, REA
 
 inline void calc_hll_speeds(REAL* ul, REAL* ur, REAL *cl, REAL* cr, int dir) {
 #ifdef USE_HLL_SPEED_BKW
-// Speeds from 
+// Speeds from
 // A multiwave approximate Riemann solver for ideal MHD based on relaxation I and II
 // Bouchut, Klingenberg and Waagan (2007 and 2012)
 
@@ -498,8 +500,7 @@ inline void calc_hll_speeds(REAL* ul, REAL* ur, REAL *cl, REAL* cr, int dir) {
       // Default?
   }
 
-
-  REAL dp_drho_cs_l = GAMMA*P_l/rho_l; 
+  REAL dp_drho_cs_l = GAMMA*P_l/rho_l;
   REAL dp_drho_cs_r = GAMMA*P_r/rho_r;
   REAL pi_l = P_l + 0.5 * Bsq_l - 0.5 * Bx_l * Bx_l;
   REAL pi_r = P_r + 0.5 * Bsq_r - 0.5 * Bx_r * Bx_r;
@@ -525,10 +526,10 @@ inline void calc_hll_speeds(REAL* ul, REAL* ur, REAL *cl, REAL* cr, int dir) {
   // eq. 3.13
   *cl = -rho_l * a_0l - alpha * rho_l * (PLUS(u_l - u_r) + PLUS(pi_r - pi_l) / (rho_l * a_ql + rho_r * a_qr));
   *cr = rho_r * a_0r + alpha * rho_r * (PLUS(u_l - u_r) + PLUS(pi_l - pi_r) / (rho_l * a_ql + rho_r * a_qr));
-  return; 
+  return;
 
 #elif defined USE_HLL_SPEED_KUSANO
-  // Speeds for the HLL solver are from equation (12) of 
+  // Speeds for the HLL solver are from equation (12) of
   // A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics
   // Miyoshi and Kusano (2005)
 
@@ -566,7 +567,7 @@ inline void calc_hll_speeds(REAL* ul, REAL* ur, REAL *cl, REAL* cr, int dir) {
 #endif //USE_HLL_SPEED_KUSANO
 }
 
- 
+
 #endif //USE_BOUNDARY_FLUX_HLL
 #endif //USE_PERIODIC
 
@@ -577,7 +578,7 @@ inline void add_surface_terms(REAL time, uint ix, uint iy, uint iz, const global
 #ifdef USE_BOUNDARY_FLUX_HLL
   int i;
   REAL um[NUM_TOTAL_VARS] = {0.0};
-  get_field(ix, iy, iz, 0, 0, 0, u, um); 
+  get_field(ix, iy, iz, 0, 0, 0, u, um);
   REAL4 b_bound = b_boundary(ix, iy, iz, time);
   REAL4 u_bound = u_boundary(ix, iy, iz, time);
   REAL rho_bound = rho_boundary(ix, iy, iz, time); //(REAL)1;
@@ -604,9 +605,9 @@ inline void add_surface_terms(REAL time, uint ix, uint iy, uint iz, const global
   REAL fluxb[NUM_CONSERVED_VARS] = {0.0};
 
   REAL alx, arx, aly, ary, alz, arz;
- 
+
   if (check_bound_xr(ix, 1)) {
-    calc_hll_speeds(um, ub, &alx, &arx, 0);    
+    calc_hll_speeds(um, ub, &alx, &arx, 0);
     calc_flux_f(um, fluxm);
     calc_flux_f(ub, fluxb);
     calc_num_flux(alx, arx, um, ub, fluxm, fluxb, flux);
@@ -622,7 +623,7 @@ inline void add_surface_terms(REAL time, uint ix, uint iy, uint iz, const global
       du_dt[i] += (REAL)(M_INV[0] / DX) * (flux[i] - fluxm[i]);
   }
 
-  if (check_bound_yr(iy, 1)) { 
+  if (check_bound_yr(iy, 1)) {
     calc_hll_speeds(um, ub, &aly, &ary, 1);
     calc_flux_g(um, fluxm);
     calc_flux_g(ub, fluxb);

include_testcases/ideal_MHD_alfven_periodic.h

@@ -90,7 +90,6 @@ inline REAL rho_boundary(uint ix, uint iy, uint iz, REAL time) {
 }
 
 
-
 /* initial conditions for the simulation */
 
 inline REAL4 b_initial(uint ix, uint iy, uint iz, REAL time){
@@ -104,11 +103,10 @@ inline REAL4 u_initial(uint ix, uint iy, uint iz, REAL 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);
-}	
+}

include_testcases/ideal_MHD_far_dipole.h

@@ -0,0 +1,106 @@
+// This project is licensed under the terms of the Creative Commons CC BY-NC-ND 4.0 license.
+
+//--------------------------------------------------------------------------------------------------
+// Functions for field initialisation
+//--------------------------------------------------------------------------------------------------
+
+
+// Example for inhomogeneous magnetic fields
+// A dipole outside of the domain
+// paramters
+REAL4 constant M = {0.0, 0.0, 1.0, 0.0};
+REAL4 constant X_0 = {-1.0, 0.5, 0.0, 0.0};
+REAL constant A0 = 0.1;
+REAL constant omega = 6.28 * 10;
+/*
+There is no analytical solution for this simulation
+*/
+inline REAL4 u_analytical(uint ix, uint iy, uint iz, REAL time) {
+
+  return (REAL4) {(REAL)(0), (REAL)(0), (REAL)(0), (REAL)(0)};
+}
+
+inline REAL4 b_analytical(uint ix, uint iy, uint iz, REAL time) {
+
+  return (REAL4) {(REAL)(0), (REAL)(0), (REAL)(0), (REAL)(0)};
+}
+
+// Magnetic Dipole definition out of
+// Course of Theoretical Physics: The Classical Theory of Fields
+// L.D. Landau and E. M. Lifschitz
+// (1951)
+
+inline REAL4 b_dipole(REAL4 X, REAL4 M) {
+  REAL inorm = rsqrt(X.x*X.x + X.y*X.y + X.z*X.z);
+  REAL4 n = inorm * X;
+  REAL4 H = (3 * n * (M * n) - M) * inorm * inorm * inorm;
+  return H;
+}
+
+/* initial values for p, rho, B and the velocity */
+
+inline REAL p_initial(uint ix, uint iy, uint iz, REAL time){
+
+  return (REAL)(1.0);
+}
+
+inline REAL rho_initial(uint ix, uint iy, uint iz, REAL time){
+
+  return (REAL)(1.0);
+}
+
+
+
+inline REAL4 b_initial(uint ix, uint iy, uint iz, REAL time){
+  REAL x = XMIN + ix * DX - X_0.x;
+  REAL y = YMIN + iy * DY - X_0.y;
+  REAL z = ZMIN + iz * DZ - X_0.z;
+
+  return b_dipole((REAL4){x, y, z, 0}, M);
+}
+
+inline REAL4 u_initial(uint ix, uint iy, uint iz, REAL time){
+
+  return (REAL4){0, 0, 0, 0} ;
+}
+
+
+/*
+Boundary condition of the magnetic field.
+The absolute magnetic field of the right boundary in x-direction oscillates. 
+We would like to see the waves which propagate to the left.
+*/
+inline REAL4 b_boundary(uint ix, uint iy, uint iz, REAL time) {
+  REAL x = XMIN + ix * DX;
+  REAL y = YMIN + iy * DY;
+  REAL z = ZMIN + iz * DZ;
+  REAL amplitude;
+  REAL4 init_b = b_initial(ix, iy, iz, (REAL)0.0);
+  if (x > 0.9) {
+    amplitude = 1 + A0 * sin(omega*time) * (x - 0.9) * 10;
+  }
+  else {
+    amplitude = 1;
+  }
+
+  return init_b * amplitude;
+}
+
+/*
+Boundary condition of the velocity field.
+*/
+inline REAL4 u_boundary(uint ix, uint iy, uint iz, REAL time) {
+
+  return u_initial(ix, iy, iz, time);
+}
+
+inline REAL rho_boundary(uint ix, uint iy, uint iz, REAL time) {
+
+  return rho_initial(ix, iy, iz, time);
+}
+
+inline REAL p_boundary(uint ix, uint iy, uint iz, REAL time){
+
+  return p_initial(ix, iy, iz, time);
+
+}