From 34554357f5f45baa2f6b90f6901b2693c3e4c0d2 Mon Sep 17 00:00:00 2001 From: Simon-Christian Klein Date: Wed, 5 Jun 2019 15:05:49 +0200 Subject: [PATCH 1/5] added a dissipation operator for the second order --- include/2ndOrder.h | 16 ++++++++++ include_defines/ideal_MHD_defines.h | 3 ++ kernel/kernel_time_integrator.cl | 48 +++++++++++++++++++++++++++++ 3 files changed, 67 insertions(+) diff --git a/include/2ndOrder.h b/include/2ndOrder.h index c1ce83b..63477d6 100644 --- a/include/2ndOrder.h +++ b/include/2ndOrder.h @@ -84,4 +84,20 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = { 2.0 }; +//------------------------------------------------------------------------------ +// Dissipation Operator +//------------------------------------------------------------------------------ + +/* A Dissipation Operator out of + * Stable and Accurate Artificial Dissipation (2003) + * Ken Mattsson, Magnus Svärd, and Jan Nordström */ + +#define NUM_BOUNDS_HO 1 +#define STENCIL_WIDTH_HOD 3 +REAL constant D_HO[2 * NUM_BOUNDS_HO+1][STENCIL_WIDTH_HOD] = { + {0.0, 2.0, -2.0}, + {1.0, -2.0, 1.0}, + {-2.0, 2.0, 0.0}, +}; + diff --git a/include_defines/ideal_MHD_defines.h b/include_defines/ideal_MHD_defines.h index 5365dc1..167d193 100644 --- a/include_defines/ideal_MHD_defines.h +++ b/include_defines/ideal_MHD_defines.h @@ -64,7 +64,10 @@ TODO: JANHUNEN_CENTRAL, BRACKBILL_BARNES_CENTRAL, cf. // Speeds: USE_HLL_SPEED_KUSANO USE_HLL_SPEED_BKW #define USE_HLL_SPEED_BKW +// Dissipation in the Volume +#define VOLUME_DISSIPATION_ACTIVE +#define VOL_DISS +0.01 enum Fields { // Density diff --git a/kernel/kernel_time_integrator.cl b/kernel/kernel_time_integrator.cl index 4240095..7290ba5 100644 --- a/kernel/kernel_time_integrator.cl +++ b/kernel/kernel_time_integrator.cl @@ -75,6 +75,51 @@ inline void add_volume_terms(REAL time, uint ix, uint iy, uint iz, global REAL c } } +#ifdef VOLUME_DISSIPATION_ACTIVE +inline void add_dissipative_terms(REAL time, uint ix, uint iy, uint iz, global REAL const* u, REAL *du_dt) { + int bound_x = get_bound_x(ix, NUM_BOUNDS_X); + int bound_y = get_bound_y(iy, NUM_BOUNDS_Y); + int bound_z = get_bound_z(iz, NUM_BOUNDS_Z); + + // the variables at remote positions used to compute the time derivative at (ix, iy, iz) + REAL um[NUM_TOTAL_VARS] = {0.0}; + + // the time derivative of the conserved variables, initalised as zero + REAL du[NUM_CONSERVED_VARS] = {0.0}; + + // x-direction + for (uint j = 0; j < STENCIL_WIDTH_HOD; j++) { + + get_field(ix, iy, iz, (j - (STENCIL_WIDTH_HOD - 1)/2), 0, 0, u, um); + + for (uint i = 0; i < NUM_CONSERVED_VARS; ++i) { + du_dt[i] = du_dt[i] + (REAL)(2.0/DX) * D_HO[NUM_BOUNDS_HO + bound_x][j] * um[i] * VOL_DISS; + } + } + + // y-direction + for (uint j = 0; j < STENCIL_WIDTH_HOD; j++) { + + get_field(ix, iy, iz, 0, (j - (STENCIL_WIDTH_HOD - 1) / 2), 0, u, um); + + for (uint i = 0; i < NUM_CONSERVED_VARS; ++i) { + du_dt[i] = du_dt[i] + (REAL)(2.0/DY) * D_HO[NUM_BOUNDS_HO + bound_x][j] * um[i] * VOL_DISS; + } + } + + // z-direction + for (uint j = 0; j < STENCIL_WIDTH_HOD; j++) { + + get_field(ix, iy, iz, 0, 0, (j - (STENCIL_WIDTH_HOD - 1) / 2), u, um); + + for (uint i = 0; i < NUM_CONSERVED_VARS; ++i) { + du_dt[i] = du_dt[i] + (REAL)(2.0/DZ) * D_HO[NUM_BOUNDS_HO + bound_x][j] * um[i] * VOL_DISS; + } + } + + +} +#endif // VOLUME_DISSIPATION_ACTIVE // Compute the time derivative du_dt using an SBP SAT discretisation with extended numerical fluxes. inline void compute_du_dt(REAL time, uint ix, uint iy, uint iz, global REAL const* u, REAL *du_dt) { @@ -85,6 +130,9 @@ inline void compute_du_dt(REAL time, uint ix, uint iy, uint iz, global REAL cons add_volume_terms(time, ix, iy, iz, u, du_dt); add_surface_terms(time, ix, iy, iz, u, du_dt); + #ifdef VOLUME_DISSIPATION_ACTIVE + add_dissipative_terms(time, ix, iy, iz, u, du_dt); + #endif // VOLUME_DISSIPATION_ACTIVE } //-------------------------------------------------------------------------------------------------- From 09cb4d6d212c5f657c67480cf256739e0b7090dd Mon Sep 17 00:00:00 2001 From: Simon-Christian Klein Date: Thu, 6 Jun 2019 23:37:54 +0200 Subject: [PATCH 2/5] added dissipation Operator from InductionEq --- include/4thOrder.h | 24 ++++++++++++++++++++++++ include/4thOrderExtended.h | 26 ++++++++++++++++++++++++++ include/6thOrder.h | 25 +++++++++++++++++++++++++ include/6thOrderExtended.h | 31 +++++++++++++++++++++++++++++++ 4 files changed, 106 insertions(+) diff --git a/include/4thOrder.h b/include/4thOrder.h index a4c402b..3e23e96 100644 --- a/include/4thOrder.h +++ b/include/4thOrder.h @@ -94,3 +94,27 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = { 48.0/17.0 }; + +//------------------------------------------------------------------------------ +// Coeffcients of the high-order dissipation operator +//------------------------------------------------------------------------------ + +#define STENCIL_WIDTH_HOD 3 +#define NUM_BOUNDS_HOD 4 + + +REAL constant D_HO[2*NUM_BOUNDS+1][STENCIL_WIDTH_HOD] = { + // left boundary coefficients + { 0.0, 0.0, -1.0}, + { 0.0, 2.0, -1.0}, + {-1.0, 2.0, -1.0}, + {-1.0, 2.0, -1.0}, + // central coefficients + {-1.0, 2.0, -1.0}, + // right boundary coefficients + {-1.0, 2.0, -1.0}, + {-1.0, 2.0, -1.0}, + {-1.0, 2.0, 0.0}, + {-1.0, 0.0, 0.0} +}; + diff --git a/include/4thOrderExtended.h b/include/4thOrderExtended.h index 45112f8..dfb55b8 100644 --- a/include/4thOrderExtended.h +++ b/include/4thOrderExtended.h @@ -100,4 +100,30 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = { 3.131115459882583 }; +//------------------------------------------------------------------------------ +// Coeffcients of the high-order dissipation operator +//------------------------------------------------------------------------------ + +#define STENCIL_WIDTH_HOD 3 +#define NUM_BOUNDS_HOD 6 + + +REAL constant D_HO[2*NUM_BOUNDS+1][STENCIL_WIDTH_HOD] = { + // left boundary coefficients + { 0.0, 0.0, -1.0}, + { 0.0, 2.0, -1.0}, + {-1.0, 2.0, -1.0}, + {-1.0, 2.0, -1.0}, + {-1.0, 2.0, -1.0}, + {-1.0, 2.0, -1.0}, + // central coefficients + {-1.0, 2.0, -1.0}, + // right boundary coefficients + {-1.0, 2.0, -1.0}, + {-1.0, 2.0, -1.0}, + {-1.0, 2.0, -1.0}, + {-1.0, 2.0, -1.0}, + {-1.0, 2.0, 0.0}, + {-1.0, 0.0, 0.0} +}; diff --git a/include/6thOrder.h b/include/6thOrder.h index 46d32fe..7746136 100644 --- a/include/6thOrder.h +++ b/include/6thOrder.h @@ -104,4 +104,29 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = { 3.1650670378782326 }; +//------------------------------------------------------------------------------ +// Coeffcients of the high-order dissipation operator +//------------------------------------------------------------------------------ + +#define STENCIL_WIDTH_HOD 4 +#define NUM_BOUNDS_HOD 6 + +REAL constant D_HO[2*NUM_BOUNDS_HOD+1][STENCIL_WIDTH_HOD] = { + // left boundary coefficients + {0.0, 0.0, 3.0, -1.0}, + {0.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + // central coefficients + {1.0, -3.0, 3.0, -1.0}, + // right boundary coefficients + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, 0.0}, + {1.0, -3.0, 0.0, 0.0}, +}; diff --git a/include/6thOrderExtended.h b/include/6thOrderExtended.h index 7ea0e33..561dcc8 100644 --- a/include/6thOrderExtended.h +++ b/include/6thOrderExtended.h @@ -111,3 +111,34 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = { }; +//------------------------------------------------------------------------------ +// Coeffcients of the high-order dissipation operator +//------------------------------------------------------------------------------ + +#define STENCIL_WIDTH_HOD 4 +#define NUM_BOUNDS_HOD 8 + + +REAL constant D_HO[2*NUM_BOUNDS_HOD+1][STENCIL_WIDTH_HOD] = { + // left boundary coefficients + {0.0, 0.0, 3.0, -1.0}, + {0.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + // central coefficients + {1.0, -3.0, 3.0, -1.0}, + // right boundary coefficients + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, -1.0}, + {1.0, -3.0, 3.0, 0.0}, + {1.0, -3.0, 0.0, 0.0}, +}; + From 23c5393c207d870562ce7a23b2f90658f9fad673 Mon Sep 17 00:00:00 2001 From: Simon-Christian Klein Date: Fri, 7 Jun 2019 14:40:33 +0200 Subject: [PATCH 3/5] ported the high order dissipation from InductionEq --- include/2ndOrder.h | 21 ++++--- include/2ndOrderExtended.h | 21 +++++++ include/4thOrder.h | 2 +- include/4thOrderExtended.h | 2 +- include_defines/ideal_MHD_defines.h | 2 +- kernel/artificial_dissipation.cl | 93 +++++++++++++++++++++++++++++ kernel/kernel_time_integrator.cl | 47 +-------------- 7 files changed, 128 insertions(+), 60 deletions(-) create mode 100644 kernel/artificial_dissipation.cl diff --git a/include/2ndOrder.h b/include/2ndOrder.h index 63477d6..0ba5ce7 100644 --- a/include/2ndOrder.h +++ b/include/2ndOrder.h @@ -88,16 +88,15 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = { // Dissipation Operator //------------------------------------------------------------------------------ -/* A Dissipation Operator out of - * Stable and Accurate Artificial Dissipation (2003) - * Ken Mattsson, Magnus Svärd, and Jan Nordström */ - -#define NUM_BOUNDS_HO 1 -#define STENCIL_WIDTH_HOD 3 -REAL constant D_HO[2 * NUM_BOUNDS_HO+1][STENCIL_WIDTH_HOD] = { - {0.0, 2.0, -2.0}, - {1.0, -2.0, 1.0}, - {-2.0, 2.0, 0.0}, -}; +#define NUM_BOUNDS_HOD 1 +#define STENCIL_WIDTH_HOD 2 +REAL constant D_HO[2*NUM_BOUNDS_HOD+1][STENCIL_WIDTH_HOD] = { + // left boundary coefficients + { 0.0, 1.0}, + // central coefficients + {-1.0, 1.0}, + // right boundary coefficients + {-1.0, 0.0} +}; diff --git a/include/2ndOrderExtended.h b/include/2ndOrderExtended.h index 013c4cd..37e8da4 100644 --- a/include/2ndOrderExtended.h +++ b/include/2ndOrderExtended.h @@ -91,3 +91,24 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = { 12.0/5.0 }; + +//------------------------------------------------------------------------------ +// Dissipation Operator +//------------------------------------------------------------------------------ + +#define NUM_BOUNDS_HOD 3 +#define STENCIL_WIDTH_HOD 2 + +REAL constant D_HO[2*NUM_BOUNDS_HOD+1][STENCIL_WIDTH_HOD] = { + // left boundary coefficients + { 0.0, 1.0}, + {-1.0, 1.0}, + {-1.0, 1.0}, + // central coefficients + {-1.0, 1.0}, + // right boundary coefficients + {-1.0, 1.0}, + {-1.0, 1.0}, + {-1.0, 0.0} +}; + diff --git a/include/4thOrder.h b/include/4thOrder.h index 3e23e96..61b5aa4 100644 --- a/include/4thOrder.h +++ b/include/4thOrder.h @@ -103,7 +103,7 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = { #define NUM_BOUNDS_HOD 4 -REAL constant D_HO[2*NUM_BOUNDS+1][STENCIL_WIDTH_HOD] = { +REAL constant D_HO[2*NUM_BOUNDS_HOD+1][STENCIL_WIDTH_HOD] = { // left boundary coefficients { 0.0, 0.0, -1.0}, { 0.0, 2.0, -1.0}, diff --git a/include/4thOrderExtended.h b/include/4thOrderExtended.h index dfb55b8..0d09257 100644 --- a/include/4thOrderExtended.h +++ b/include/4thOrderExtended.h @@ -108,7 +108,7 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = { #define NUM_BOUNDS_HOD 6 -REAL constant D_HO[2*NUM_BOUNDS+1][STENCIL_WIDTH_HOD] = { +REAL constant D_HO[2*NUM_BOUNDS_HOD+1][STENCIL_WIDTH_HOD] = { // left boundary coefficients { 0.0, 0.0, -1.0}, { 0.0, 2.0, -1.0}, diff --git a/include_defines/ideal_MHD_defines.h b/include_defines/ideal_MHD_defines.h index 167d193..22edded 100644 --- a/include_defines/ideal_MHD_defines.h +++ b/include_defines/ideal_MHD_defines.h @@ -67,7 +67,7 @@ TODO: JANHUNEN_CENTRAL, BRACKBILL_BARNES_CENTRAL, cf. // Dissipation in the Volume #define VOLUME_DISSIPATION_ACTIVE -#define VOL_DISS +0.01 +#define HO_DISSIPATION_FACTOR 1.0 enum Fields { // Density diff --git a/kernel/artificial_dissipation.cl b/kernel/artificial_dissipation.cl new file mode 100644 index 0000000..63cf257 --- /dev/null +++ b/kernel/artificial_dissipation.cl @@ -0,0 +1,93 @@ +// The high order artificial dissipation, ported here from InductionEq +// TODO: implement per field dissipation factors + +#ifndef HO_DISSIPATION_FACTOR + define HO_DISSIPATION_FACTOR 0.01 +#endif + +// Auxillary function for high order dissipation +// Used to compute the intermediate result D*a of equation (?) +void inline diff_HOD_x(uint ix, uint iy, uint iz, int bound_x, int start, global REAL const *d_field, REAL *HOD) { + + REAL um[NUM_TOTAL_VARS] = {0.0}; + int i, j; + + for (i = 0; i < STENCIL_WIDTH_HOD; i++) { + get_field(ix, iy, iz, i + start, 0, 0, d_field, um); + + for (j = 0; i < NUM_CONSERVED_VARS; i++) { + HOD[j] = HOD[j] + + M_INV_X[NUM_BOUNDS_HOD + bound_x] + * D_HO[NUM_BOUNDS_HOD + bound_x][start + STENCIL_WIDTH_HOD - 1] + * D_HO[NUM_BOUNDS_HOD][i] * pown(-1.0, (int)(ORDER / 2 + 1)) + * um[j]; + } + } + +} + +void inline diff_HOD_y(uint ix, uint iy, uint iz, int bound_y, int start, global REAL const *d_field, REAL *HOD) { + + REAL um[NUM_TOTAL_VARS] = {0.0}; + int i, j; + + for (i = 0; i < STENCIL_WIDTH_HOD; i++) { + get_field(ix, iy, iz, 0, i + start, 0, d_field, um); + + for (j = 0; j < NUM_CONSERVED_VARS; j++) { + HOD[j] = HOD[j] + + M_INV_Y[NUM_BOUNDS_HOD + bound_y] + * D_HO[NUM_BOUNDS_HOD + bound_y][start + STENCIL_WIDTH_HOD - 1] + * D_HO[NUM_BOUNDS_HOD][i] * pown(-1.0, (int)(ORDER / 2 + 1)) + * um[j]; + } + } + +} + +void inline diff_HOD_z(uint ix, uint iy, uint iz, int bound_z, int start, global REAL const *d_field, REAL *HOD) { + + REAL um[NUM_TOTAL_VARS] = {0.0}; + int i, j; + + for (i = 0; i < STENCIL_WIDTH_HOD; i++) { + get_field(ix, iy, iz, 0,0, i + start, d_field, um); + + for (j = 0; j < NUM_CONSERVED_VARS; j++) { + HOD[j] = HOD[j] + + M_INV_Y[NUM_BOUNDS_HOD + bound_z] + * D_HO[NUM_BOUNDS_HOD + bound_z][start + STENCIL_WIDTH_HOD - 1] + * D_HO[NUM_BOUNDS_HOD][i] * pown(-1.0, (int)(ORDER / 2 + 1)) + * um[j]; + } + } + +} + + + +// Computes the high order dissipation for fields in u at point (ix, iy, iz) +void inline high_order_dissipation(uint ix, uint iy, uint iz, global REAL const *u, REAL *HOD) { + int i, j; + uint idx = calc_idx(ix,iy,iz); + + int bound_x = get_bound_x(ix, NUM_BOUNDS_HOD); + int bound_y = get_bound_y(iy, NUM_BOUNDS_HOD); + int bound_z = get_bound_z(iz, NUM_BOUNDS_HOD); + REAL HOD_X[NUM_TOTAL_VARS] = {0.0}; + REAL HOD_Y[NUM_TOTAL_VARS] = {0.0}; + REAL HOD_Z[NUM_TOTAL_VARS] = {0.0}; + + for (i = -STENCIL_WIDTH_HOD + 1; i < 1; i++) { + diff_HOD_x(ix, iy, iz, bound_x, i, u, HOD_X); + diff_HOD_y(ix, iy, iz, bound_x, i, u, HOD_Y); + diff_HOD_z(ix, iy, iz, bound_x, i, u, HOD_Z); + for(j = 0; j < NUM_CONSERVED_VARS;j++) { + HOD[j] = HOD[j] + (REAL)(HO_DISSIPATION_FACTOR) * HOD_X[j] + + (REAL)(HO_DISSIPATION_FACTOR) * HOD_Y[j] + + (REAL)(HO_DISSIPATION_FACTOR) * HOD_Z[j]; + } + } + +} + diff --git a/kernel/kernel_time_integrator.cl b/kernel/kernel_time_integrator.cl index 7290ba5..ca8bb6a 100644 --- a/kernel/kernel_time_integrator.cl +++ b/kernel/kernel_time_integrator.cl @@ -75,51 +75,6 @@ inline void add_volume_terms(REAL time, uint ix, uint iy, uint iz, global REAL c } } -#ifdef VOLUME_DISSIPATION_ACTIVE -inline void add_dissipative_terms(REAL time, uint ix, uint iy, uint iz, global REAL const* u, REAL *du_dt) { - int bound_x = get_bound_x(ix, NUM_BOUNDS_X); - int bound_y = get_bound_y(iy, NUM_BOUNDS_Y); - int bound_z = get_bound_z(iz, NUM_BOUNDS_Z); - - // the variables at remote positions used to compute the time derivative at (ix, iy, iz) - REAL um[NUM_TOTAL_VARS] = {0.0}; - - // the time derivative of the conserved variables, initalised as zero - REAL du[NUM_CONSERVED_VARS] = {0.0}; - - // x-direction - for (uint j = 0; j < STENCIL_WIDTH_HOD; j++) { - - get_field(ix, iy, iz, (j - (STENCIL_WIDTH_HOD - 1)/2), 0, 0, u, um); - - for (uint i = 0; i < NUM_CONSERVED_VARS; ++i) { - du_dt[i] = du_dt[i] + (REAL)(2.0/DX) * D_HO[NUM_BOUNDS_HO + bound_x][j] * um[i] * VOL_DISS; - } - } - - // y-direction - for (uint j = 0; j < STENCIL_WIDTH_HOD; j++) { - - get_field(ix, iy, iz, 0, (j - (STENCIL_WIDTH_HOD - 1) / 2), 0, u, um); - - for (uint i = 0; i < NUM_CONSERVED_VARS; ++i) { - du_dt[i] = du_dt[i] + (REAL)(2.0/DY) * D_HO[NUM_BOUNDS_HO + bound_x][j] * um[i] * VOL_DISS; - } - } - - // z-direction - for (uint j = 0; j < STENCIL_WIDTH_HOD; j++) { - - get_field(ix, iy, iz, 0, 0, (j - (STENCIL_WIDTH_HOD - 1) / 2), u, um); - - for (uint i = 0; i < NUM_CONSERVED_VARS; ++i) { - du_dt[i] = du_dt[i] + (REAL)(2.0/DZ) * D_HO[NUM_BOUNDS_HO + bound_x][j] * um[i] * VOL_DISS; - } - } - - -} -#endif // VOLUME_DISSIPATION_ACTIVE // Compute the time derivative du_dt using an SBP SAT discretisation with extended numerical fluxes. inline void compute_du_dt(REAL time, uint ix, uint iy, uint iz, global REAL const* u, REAL *du_dt) { @@ -131,7 +86,7 @@ inline void compute_du_dt(REAL time, uint ix, uint iy, uint iz, global REAL cons add_volume_terms(time, ix, iy, iz, u, du_dt); add_surface_terms(time, ix, iy, iz, u, du_dt); #ifdef VOLUME_DISSIPATION_ACTIVE - add_dissipative_terms(time, ix, iy, iz, u, du_dt); + high_order_dissipation(ix, iy, iz, u, du_dt); #endif // VOLUME_DISSIPATION_ACTIVE } From ac59cdddb9ba869c38181ad900ae3a8a2e3adc07 Mon Sep 17 00:00:00 2001 From: Simon-Christian Klein Date: Wed, 19 Jun 2019 09:38:16 +0200 Subject: [PATCH 4/5] dedicated dissipation factor for every conserved variable --- include_defines/ideal_MHD_defines.h | 2 +- kernel/artificial_dissipation.cl | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/include_defines/ideal_MHD_defines.h b/include_defines/ideal_MHD_defines.h index 22edded..e2cf417 100644 --- a/include_defines/ideal_MHD_defines.h +++ b/include_defines/ideal_MHD_defines.h @@ -67,7 +67,7 @@ TODO: JANHUNEN_CENTRAL, BRACKBILL_BARNES_CENTRAL, cf. // Dissipation in the Volume #define VOLUME_DISSIPATION_ACTIVE -#define HO_DISSIPATION_FACTOR 1.0 +#define HO_DISSIPATION_FACTOR {0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01} enum Fields { // Density diff --git a/kernel/artificial_dissipation.cl b/kernel/artificial_dissipation.cl index 63cf257..a8f6226 100644 --- a/kernel/artificial_dissipation.cl +++ b/kernel/artificial_dissipation.cl @@ -70,7 +70,7 @@ void inline diff_HOD_z(uint ix, uint iy, uint iz, int bound_z, int start, global void inline high_order_dissipation(uint ix, uint iy, uint iz, global REAL const *u, REAL *HOD) { int i, j; uint idx = calc_idx(ix,iy,iz); - + REAL ho_dissipation_factor[NUM_CONSERVED_VARS] = HO_DISSIPATION_FACTOR; int bound_x = get_bound_x(ix, NUM_BOUNDS_HOD); int bound_y = get_bound_y(iy, NUM_BOUNDS_HOD); int bound_z = get_bound_z(iz, NUM_BOUNDS_HOD); @@ -83,9 +83,9 @@ void inline high_order_dissipation(uint ix, uint iy, uint iz, global REAL const diff_HOD_y(ix, iy, iz, bound_x, i, u, HOD_Y); diff_HOD_z(ix, iy, iz, bound_x, i, u, HOD_Z); for(j = 0; j < NUM_CONSERVED_VARS;j++) { - HOD[j] = HOD[j] + (REAL)(HO_DISSIPATION_FACTOR) * HOD_X[j] - + (REAL)(HO_DISSIPATION_FACTOR) * HOD_Y[j] - + (REAL)(HO_DISSIPATION_FACTOR) * HOD_Z[j]; + HOD[j] = HOD[j] + (REAL)(ho_dissipation_factor[j]) * HOD_X[j] + + (REAL)(ho_dissipation_factor[j]) * HOD_Y[j] + + (REAL)(ho_dissipation_factor[j]) * HOD_Z[j]; } } From d8619d92819c369aabb8906ff97aedd45bdc2aa3 Mon Sep 17 00:00:00 2001 From: Simon-Christian Klein Date: Fri, 28 Jun 2019 15:48:58 +0200 Subject: [PATCH 5/5] corrected a error in the new routines for artificail dissipation --- kernel/artificial_dissipation.cl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/kernel/artificial_dissipation.cl b/kernel/artificial_dissipation.cl index a8f6226..e82ded0 100644 --- a/kernel/artificial_dissipation.cl +++ b/kernel/artificial_dissipation.cl @@ -80,8 +80,8 @@ void inline high_order_dissipation(uint ix, uint iy, uint iz, global REAL const for (i = -STENCIL_WIDTH_HOD + 1; i < 1; i++) { diff_HOD_x(ix, iy, iz, bound_x, i, u, HOD_X); - diff_HOD_y(ix, iy, iz, bound_x, i, u, HOD_Y); - diff_HOD_z(ix, iy, iz, bound_x, i, u, HOD_Z); + diff_HOD_y(ix, iy, iz, bound_y, i, u, HOD_Y); + diff_HOD_z(ix, iy, iz, bound_z, i, u, HOD_Z); for(j = 0; j < NUM_CONSERVED_VARS;j++) { HOD[j] = HOD[j] + (REAL)(ho_dissipation_factor[j]) * HOD_X[j] + (REAL)(ho_dissipation_factor[j]) * HOD_Y[j]