Skip to content

High order artificial dissipation ported from InductionEq #39

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
simonius wants to merge 5 commits into
base: master
Choose a base branch
from
Open

High order artificial dissipation ported from InductionEq #39

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
3455435
added a dissipation operator for the second order
Jun 5, 2019
09cb4d6
added dissipation Operator from InductionEq
simonius Jun 6, 2019
23c5393
ported the high order dissipation from InductionEq
Jun 7, 2019
ac59cdd
dedicated dissipation factor for every conserved variable
Jun 19, 2019
d8619d9
corrected a error in the new routines for artificail dissipation
Jun 28, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 15 additions & 0 deletions include/2ndOrder.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,4 +84,19 @@ REAL constant M_INV_B[2*NUM_BOUNDS_B+1] = {
2.0
};

//------------------------------------------------------------------------------
// Dissipation Operator
//------------------------------------------------------------------------------

#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}
};

21 changes: 21 additions & 0 deletions include/2ndOrderExtended.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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}
};

24 changes: 24 additions & 0 deletions include/4thOrder.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_HOD+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}
};

26 changes: 26 additions & 0 deletions include/4thOrderExtended.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_HOD+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}
};

25 changes: 25 additions & 0 deletions include/6thOrder.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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},
};

31 changes: 31 additions & 0 deletions include/6thOrderExtended.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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},
};

3 changes: 3 additions & 0 deletions include_defines/ideal_MHD_defines.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 HO_DISSIPATION_FACTOR {0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01}

enum Fields {
// Density
Expand Down
93 changes: 93 additions & 0 deletions kernel/artificial_dissipation.cl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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);
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);
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_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]
+ (REAL)(ho_dissipation_factor[j]) * HOD_Z[j];
}
}

}

3 changes: 3 additions & 0 deletions kernel/kernel_time_integrator.cl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,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
high_order_dissipation(ix, iy, iz, u, du_dt);
#endif // VOLUME_DISSIPATION_ACTIVE
}

//--------------------------------------------------------------------------------------------------
Expand Down