-
Notifications
You must be signed in to change notification settings - Fork 1
Expand file tree
/
Copy pathsetup.cpp
More file actions
141 lines (119 loc) · 3.76 KB
/
setup.cpp
File metadata and controls
141 lines (119 loc) · 3.76 KB
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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
#include <algorithm>
#include "idefix.hpp"
#include "setup.hpp"
#include "soundspeed.hpp"
real sigma0Glob;
real sigmaSlopeGlob;
real h0Glob;
real betaGlob;
real HidealGlob;
real AmMidGlob;
real gammaGlob;
real densityFloorGlob;
SoundSpeed *ss;
void MySoundSpeed(DataBlock &data, const real t, IdefixArray3D<real> &cs) {
// Call our class
ss->Compute(cs);
}
// User-defined boundaries
void UserdefBoundary(Hydro *hydro, int dir, BoundarySide side, real t) {
IdefixArray4D<real> Vc = hydro->Vc;
auto *data = hydro->data;
IdefixArray1D<real> x1 = data->x[IDIR];
IdefixArray1D<real> x3 = data->x[KDIR];
if(dir==IDIR) {
int ighost,ibeg,iend;
if(side == left) {
ighost = data->beg[IDIR];
ibeg = 0;
iend = data->beg[IDIR];
idefix_for("UserDefBoundary",
0, data->np_tot[KDIR],
0, data->np_tot[JDIR],
ibeg, iend,
KOKKOS_LAMBDA (int k, int j, int i) {
real R=x1(i);
real z=x3(k);
real Vk = 1.0/sqrt(R);
Vc(RHO,k,j,i) = Vc(RHO,k,j,2*ighost - i +1);
Vc(VX1,k,j,i) = - Vc(VX1,k,j,2*ighost - i +1);
Vc(VX2,k,j,i) = Vk;
Vc(VX3,k,j,i) = Vc(VX3,k,j,2*ighost - i +1);
});
}
else if(side==right) {
ighost = data->end[IDIR]-1;
ibeg=data->end[IDIR];
iend=data->np_tot[IDIR];
idefix_for("UserDefBoundary",
0, data->np_tot[KDIR],
0, data->np_tot[JDIR],
ibeg, iend,
KOKKOS_LAMBDA (int k, int j, int i) {
real R=x1(i);
real z=x3(k);
real Vk = 1.0/sqrt(R);
Vc(RHO,k,j,i) = Vc(RHO,k,j,ighost);
Vc(VX1,k,j,i) = Vc(VX1,k,j,ighost);
Vc(VX2,k,j,i) = Vk;
Vc(VX3,k,j,i) = Vc(VX3,k,j,ighost);
});
}
}
}
void FargoVelocity(DataBlock &data, IdefixArray2D<real> &Vphi) {
IdefixArray1D<real> x1 = data.x[IDIR];
idefix_for("FargoVphi",0,data.np_tot[KDIR], 0, data.np_tot[IDIR],
KOKKOS_LAMBDA (int k, int i) {
Vphi(k,i) = 1.0/sqrt(x1(i));
});
}
// Default constructor
// Initialisation routine. Can be used to allocate
// Arrays or variables which are used later on
Setup::Setup(Input &input, Grid &grid, DataBlock &data, Output &output)// : m_planet(0)//, Planet &planet)
{
// Set the function for userdefboundary
data.hydro->EnrollUserDefBoundary(&UserdefBoundary);
data.hydro->EnrollIsoSoundSpeed(&MySoundSpeed);
if(data.haveFargo)
data.fargo->EnrollVelocity(&FargoVelocity);
sigma0Glob = input.Get<real>("Setup","sigma0",0);
sigmaSlopeGlob = input.Get<real>("Setup","sigmaSlope",0);
h0Glob = input.Get<real>("Setup","h0",0);
ss = new SoundSpeed(input, data);
}
// Destructor
Setup::~Setup() {
delete ss;
}
// This routine initialize the flow
// Note that data is on the device.
// One can therefore define locally
// a datahost and sync it, if needed
void Setup::InitFlow(DataBlock &data) {
// Create a host copy
DataBlockHost d(data);
real h0=h0Glob;
real sigma0=sigma0Glob;
real sigmaSlope = sigmaSlopeGlob;
for(int k = 0; k < d.np_tot[KDIR] ; k++) {
for(int j = 0; j < d.np_tot[JDIR] ; j++) {
for(int i = 0; i < d.np_tot[IDIR] ; i++) {
real R=d.x[IDIR](i);
real z=d.x[KDIR](k);
real Vk=1.0/sqrt(R);
real cs2=(h0*Vk)*(h0*Vk);
d.Vc(RHO,k,j,i) = sigma0/sqrt(2.0*M_PI)/(h0*R)*pow(R,-sigmaSlope) * exp(1.0/ (cs2) * (1.0/sqrt(R*R+z*z)-1.0/R)) ;
d.Vc(VX1,k,j,i) = 0.0;//+sin(8*d.x[JDIR](j));
d.Vc(VX3,k,j,i) = 0.0;
d.Vc(VX2,k,j,i) = Vk*sqrt( R / sqrt(R*R + z*z) -(2.0+sigmaSlope)*h0*h0);
}
}
}
// Send it all, if needed
d.SyncToDevice();
}
// Analyse data to produce an output
void MakeAnalysis(DataBlock & data) {
}