-
Notifications
You must be signed in to change notification settings - Fork 2
/
viscosity.cpp
104 lines (81 loc) · 3.61 KB
/
viscosity.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
/*
Crown Copyright 2012 AWE.
This file is part of CloverLeaf.
CloverLeaf is free software: you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your option)
any later version.
CloverLeaf is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
details.
You should have received a copy of the GNU General Public License along with
CloverLeaf. If not, see http://www.gnu.org/licenses/.
*/
#include "viscosity.h"
// @brief Fortran viscosity kernel.
// @author Wayne Gaudin
// @details Calculates an artificial viscosity using the Wilkin's method to
// smooth out shock front and prevent oscillations around discontinuities.
// Only cells in compression will have a non-zero value.
void viscosity_kernel(int x_min, int x_max, int y_min, int y_max,
Kokkos::View<double*>& celldx,
Kokkos::View<double*>& celldy,
Kokkos::View<double**>& density0,
Kokkos::View<double**>& pressure,
Kokkos::View<double**>& viscosity,
Kokkos::View<double**>& xvel0,
Kokkos::View<double**>& yvel0) {
// DO k=y_min,y_max
// DO j=x_min,x_max
Kokkos::MDRangePolicy<Kokkos::Rank<2>> policy({x_min+1, y_min+1}, {x_max+2, y_max+2});
Kokkos::parallel_for("viscosity", policy, KOKKOS_LAMBDA(const int j, const int k) {
double ugrad = (xvel0(j+1,k )+xvel0(j+1,k+1))-(xvel0(j ,k )+xvel0(j ,k+1));
double vgrad = (yvel0(j ,k+1)+yvel0(j+1,k+1))-(yvel0(j ,k )+yvel0(j+1,k ));
double div = (celldx(j)*(ugrad)+ celldy(k)*(vgrad));
double strain2 = 0.5*(xvel0(j, k+1) + xvel0(j+1,k+1)-xvel0(j ,k )-xvel0(j+1,k ))/celldy(k)
+ 0.5*(yvel0(j+1,k ) + yvel0(j+1,k+1)-yvel0(j ,k )-yvel0(j ,k+1))/celldx(j);
double pgradx=(pressure(j+1,k)-pressure(j-1,k))/(celldx(j)+celldx(j+1));
double pgrady=(pressure(j,k+1)-pressure(j,k-1))/(celldy(k)+celldy(k+1));
double pgradx2 = pgradx*pgradx;
double pgrady2 = pgrady*pgrady;
double limiter = ((0.5*(ugrad)/celldx(j))*pgradx2+(0.5*(vgrad)/celldy(k))*pgrady2+strain2*pgradx*pgrady)
/MAX(pgradx2+pgrady2,1.0e-16);
if ((limiter > 0.0) || (div >= 0.0)) {
viscosity(j,k) = 0.0;
} else {
double dirx=1.0;
if (pgradx < 0.0) dirx=-1.0;
pgradx = dirx*MAX(1.0e-16,fabs(pgradx));
double diry=1.0;
if (pgradx < 0.0) diry=-1.0;
pgrady = diry*MAX(1.0e-16,fabs(pgrady));
double pgrad = sqrt(pgradx*pgradx+pgrady*pgrady);
double xgrad = fabs(celldx(j)*pgrad/pgradx);
double ygrad = fabs(celldy(k)*pgrad/pgrady);
double grad = MIN(xgrad,ygrad);
double grad2 = grad*grad;
viscosity(j,k)=2.0*density0(j,k)*grad2*limiter*limiter;
}
});
}
// @brief Driver for the viscosity kernels
// @author Wayne Gaudin
// @details Selects the user specified kernel to caluclate the artificial
// viscosity.
void viscosity(global_variables& globals) {
for (int tile = 0; tile < globals.tiles_per_chunk; ++tile) {
viscosity_kernel(
globals.chunk.tiles[tile].t_xmin,
globals.chunk.tiles[tile].t_xmax,
globals.chunk.tiles[tile].t_ymin,
globals.chunk.tiles[tile].t_ymax,
globals.chunk.tiles[tile].field.celldx,
globals.chunk.tiles[tile].field.celldy,
globals.chunk.tiles[tile].field.density0,
globals.chunk.tiles[tile].field.pressure,
globals.chunk.tiles[tile].field.viscosity,
globals.chunk.tiles[tile].field.xvel0,
globals.chunk.tiles[tile].field.yvel0);
}
}