-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathPhysics.cpp
278 lines (223 loc) · 13.2 KB
/
Physics.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
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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
// Copyright(C) 2021 Salvatore Cielo, LRZ
// Copyright(C) 2022 Alexander Pöppl, Intel Corp.
//
// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the
// License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an
// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific
// language governing permissions and limitations under the License.
#include "Physics.hpp"
#include "Metric.hpp"
using namespace sycl;
//-- Manual 3-vectors operation in GR
SYCL_EXTERNAL void matMul(field m[9], field *vIn, field *vOut, id<1> gid, unsigned offset){
const unsigned id0 = gid[0], id1 = id0+offset, id2 = id1+offset;
vOut[id0] = m[0]*vIn[id0] + m[1] *vIn[id1] + m[2] *vIn[id2];
vOut[id1] = m[3]*vIn[id0] + m[4] *vIn[id1] + m[5] *vIn[id2];
vOut[id2] = m[6]*vIn[id0] + m[7] *vIn[id1] + m[8] *vIn[id2];
} // Multiplication by STATIC 3x3 matrix. Sorry it can't be inlined
SYCL_EXTERNAL field dot(field* u, field* v, unsigned dims, id<1> gid, unsigned offset){
field sum = 0.0;
for(unsigned i=0; i<dims; i++){ sum+=u[gid[0]+i*offset]*v[gid[0]+i*offset]; }
return sum;
} // Sorry it can't be inlined. Otherwise it should go straigth to MKL
SYCL_EXTERNAL void cross(field *u, field *v, field *r , id<1> gid, unsigned offset){
const unsigned id0 = gid[0], id1 = id0+offset, id2 = id1+offset;
r[id0] = u[id1]*v[id2] - v[id1]*u[id2];
r[id1] = u[id2]*v[id0] - v[id2]*u[id0];
r[id2] = u[id0]*v[id1] - v[id0]*u[id1];
}
//-- Conversion of quantities
SYCL_EXTERNAL void prim2cons(id<1> myId, unsigned n, field_array v, field_array u, Metric &g){
const field gDet = g.gDet();
unsigned const gid = myId[0];
#if PHYSICS==MHD
field vCon[]={v[VX][gid], v[VY][gid], v[VZ][gid]}, vCov[3]; g.con2Cov(vCon, vCov); field v2=dot(vCon,vCov);
field bCon[]={v[BX][gid], v[BY][gid], v[BZ][gid]}, bCov[3]; g.con2Cov(bCon, bCov); field b2=dot(bCon,bCov);
field rh=v[RH][gid], pg=v[PG][gid], pt=pg+0.5*b2, wt=0.5*rh*v2+GAMMA1*pg+b2;
u[RH][gid]=gDet*rh; u[PG][gid]=gDet*(wt-pt);
u[VX][gid]=gDet*rh*vCov[0]; u[VY][gid]=gDet*rh*vCov[1]; u[VZ][gid]=gDet*rh*vCov[2];
u[BX][gid]=gDet *bCov[0]; u[BY][gid]=gDet *bCov[1]; u[BZ][gid]=gDet *bCov[2];
#elif PHYSICS==GRMHD
field rh=v[RH][gid], pg=v[PG][gid];
field vCon[]={v[VX][gid], v[VY][gid], v[VZ][gid]}, vCov[3]; g.con2Cov(vCon, vCov); field u2=dot(vCon,vCov); field glf=sycl::sqrt(1.0+u2);
field bCon[]={v[BX][gid], v[BY][gid], v[BZ][gid]}, bCov[3]; g.con2Cov(bCon, bCov); field b2=dot(bCon,bCov);
vCov[0] *= 1.0/glf; vCov[1] *= 1.0/glf; vCov[2] *= 1.0/glf;
vCon[0] *= 1.0/glf; vCon[1] *= 1.0/glf; vCon[2] *= 1.0/glf;
field eCov[3]; cross(vCon, bCon, eCov); eCov[0]*=-gDet; eCov[1]*=-gDet; eCov[2]*=-gDet;
field eCon[3]; g.cov2Con(eCov, eCon); field e2=dot(eCov,eCon), uem=0.5*(e2+b2);
field sCov[3]; cross(eCon, bCon, sCov); sCov[0]*= gDet; sCov[1]*= gDet; sCov[2]*= gDet;
field sCon[3]; g.cov2Con(sCov, sCon);
field d=rh*glf, h=rh+GAMMA1*pg, w=h*glf*glf, w1=d*u2/(1.+glf)+GAMMA1*pg*glf*glf; // W'=W-D
u[RH][gid] = d; u[PG][gid] = w1-pg+uem;
u[VX][gid] = w*vCov[0]+sCov[0]; u[VY][gid] = w*vCov[1]+sCov[1]; u[VZ][gid] = w*vCov[2]+sCov[2];
u[BX][gid] = bCon[0] ; u[BY][gid] = bCon[1] ; u[BZ][gid] = bCon[2] ;
for(unsigned short iFld=0; iFld<FLD_TOT; ++iFld){ u[iFld][gid] *= gDet; }
#endif // PHYSICS
}
SYCL_EXTERNAL void cons2prim(id<1> myId, unsigned n, field_array u, field_array v, Metric &g){
unsigned const gid = myId[0];
field const gDet1 = g.gDet1();
field sCov[3]={u[VX][gid]*gDet1,u[VY][gid]*gDet1,u[VZ][gid]*gDet1}, sCon[3]; g.cov2Con(sCov, sCon); field s2=dot(sCov, sCon);
field bCon[3]={u[BX][gid]*gDet1,u[BY][gid]*gDet1,u[BZ][gid]*gDet1}, bCov[3]; g.con2Cov(bCon, bCov); field b2=dot(bCov, bCon);
#if PHYSICS==MHD
field rh = u[RH][gid]*gDet1, et=u[PG][gid]*gDet1, rh1=1./rh, pg = (GAMMA-1.0)*(et-.5*(rh1*s2+b2));
v[RH][gid] = rh; v[PG][gid] = pg;
v[VX][gid] = sCon[0]*rh1; v[VY][gid] = sCon[1]*rh1; v[VZ][gid] = sCon[2]*rh1;
v[BX][gid] = bCon[0] ; v[BY][gid] = bCon[1] ; v[BZ][gid] = bCon[2] ;
#elif PHYSICS==GRMHD
field sb=dot(sCov,bCon), sb2=sb*sb, b2st2=sycl::max(s2*b2-sb2, 0.);
field d=u[RH][gid]*gDet1, et1=u[PG][gid]*gDet1, w1, u2, glf;
// The folowing two are alternative: Initial guess from old values ...
//field uCon[] = {v[VX][gid], v[VY][gid], v[VZ][gid]}, uCov[3]; g.con2Cov(uCon, uCov); field u2=dot(uCov,uCon); glf=sycl::sqrt(1.+u2);
//w1 = d*u2/(1.0+glf)+GAMMA1*v[PG][gid]*glf*glf;
// ... OR from quadratic equation (usually better)
field com = (et1+d-b2); w1 = 4.*com*com - 3.*(s2-b2*(2.*et1+d-b2));
w1 = sycl::max( (2.*com+sycl::sqrt(sycl::max(w1,0.)))/3.-d, 0.);
//-- Undetermined iteration. ALL LOCAL, luckily
field w, vv2, pg, fw, dv2, dpg, dfw, dw1; const field tol=1.e-9;
for(unsigned iter=0; iter<20; ++iter){
w = w1+d; vv2 = w*w*s2+(2.*w+b2)*sb2; com = w*(w+b2); u2 = vv2/(com*com-vv2); glf = sycl::sqrt(1.+u2);
pg = (w1-d*u2/(1.+glf))/(GAMMA1*glf*glf); com = w+b2; fw = w1-et1-pg + .5*(b2 + b2st2/(com*com));
dv2 = -2.*( s2+sb2* (3.*w*(com)+b2*b2)/(w*w*w) )/ (com*com*com);
dpg = 1./(GAMMA1*glf*glf) - glf*(.5*d/GAMMA1+glf*pg)*dv2; dfw = 1.-dpg-b2st2/(com*com*com); dw1 = -fw/dfw;
if (sycl::fabs(dw1) < tol*w1){ break; }else{ w1+= dw1;}
}
field rh = d/glf; pg = sycl::max(pg, (field)PGFLOOR);
field vCon[3]={ (sCon[0]+sb*bCon[0]/w)/(w+b2), (sCon[1]+sb*bCon[1]/w)/(w+b2), (sCon[2]+sb*bCon[2]/w)/(w+b2)};
v[RH][gid] = rh; v[PG][gid] = pg;
v[VX][gid] = vCon[0]*glf; v[VY][gid] = vCon[1]*glf; v[VZ][gid] = vCon[2]*glf;
v[BX][gid] = bCon[0] ; v[BY][gid] = bCon[1] ; v[BZ][gid] = bCon[2] ;
#endif // PHYSICS==GRMHD
}
//-- Fluxes and characteristic velocities.
// IMPORTANT: all local quantities, they have been sampled --> access simply by eg. f[VX]
SYCL_EXTERNAL void physicalFlux(int dir, Metric &g, field vD[FLD_TOT], field uD[FLD_TOT], field f[FLD_TOT], field vf[2], field vt[2] ){
field betai[3], gCov[9], gCon[9]; g.g3DCov(gCov); g.beta(betai);
field const alpha = g.alpha(), gDet = g.g3DCon(gCon);
const short k1 = (dir+1)-1, k2 = ((dir+1)%3)-1, k3 = ((dir+2)%3)-2;
#if PHYSICS==MHD
//-- Assignments
field vCov[3], vCon[3] = {vD[VX], vD[VY], vD[VZ]}; matMul(gCov, vCon, vCov);
field bCov[3], bCon[3] = {vD[BX], vD[BY], vD[BZ]}; matMul(gCov, bCon, bCov);
const field v2 = dot(vCon,vCov), b2 = dot(bCon,bCov), vb = dot(vCov,bCon);
const field rh = vD[RH], pg = vD[PG], pt = pg+.5*b2, wt =.5*rh*v2+GAMMA1*pg+b2;
//-- Conserved
uD[RH] = rh ; uD[VX] = rh *vCov[0]; uD[VY] = rh *vCov[1]; uD[VZ] = rh *vCov[2];
uD[PG] = wt - pt; uD[BX] = vD[BX]; uD[BY] = vD[BY]; uD[BZ] = vD[BZ];
for(unsigned iVar = 0; iVar< FLD_TOT; ++iVar){ uD[iVar] *= gDet;}
//-- Fluxes
f[VX]= rh*vCon[dir]*vCov[0] - bCon[dir]*bCov[0]; f[BX]= vCon[dir]*bCon[0]-bCon[dir]*vCon[0];
f[VY]= rh*vCon[dir]*vCov[1] - bCon[dir]*bCov[1]; f[BY]= vCon[dir]*bCon[1]-bCon[dir]*vCon[1];
f[VZ]= rh*vCon[dir]*vCov[2] - bCon[dir]*bCov[2]; f[BZ]= vCon[dir]*bCon[2]-bCon[dir]*vCon[2];
f[VX+k1] += pt;
f[PG] = wt*vD[VX+dir] - vb*vD[BX+dir]; f[RH] = rh * vD[BX+dir];
for(unsigned iVar = 0; iVar< FLD_TOT; ++iVar){ f[iVar] *= gDet;}
//-- Fast magnetosonic speeds (vCon along direction dir)
const field c2 = GAMMA *pg, a2 = c2+b2;
const field comfort = gCon[dir+3*dir]*a2*a2 - 4.0*c2*bCon[dir]*bCon[dir];
const field vfd = sycl::sqrt( 0.5*( a2+sycl::sqrt( sycl::max((field)0.0,(field)comfort) ) )/rh );
vf[0] = vCon[dir]+vfd; vf[1] = vCon[dir]-vfd;
//-- Transverse speeds
vt[0] = vCon[1+k2]; vt[1] = vCon[2+k3];
#elif PHYSICS==GRMHD
//-- Assignments
field rh = vD[RH], pg = vD[PG], vCon[] = {vD[VX], vD[VY], vD[VZ]}, bCon[] = {vD[BX], vD[BY], vD[BZ]};
field vCov[3], bCov[3]; matMul(gCov, vCon, vCov); matMul(gCov, bCon, bCov);
const field u2 = dot(vCov,vCon), b2 = dot(bCon,bCov), glfInv = sycl::rsqrt(1+u2), glf = 1.0/glfInv;
vCon[0] *= glfInv; vCon[1] *= glfInv; vCon[2] *= glfInv;
vCov[0] *= glfInv; vCov[1] *= glfInv; vCov[2] *= glfInv;
field eCov[3]; cross(vCon, bCon, eCov); eCov[0]*=-gDet; eCov[1]*=-gDet; eCov[2]*=-gDet;
field eCon[3]; matMul(gCon, eCov, eCon); const field e2 = dot(eCov,eCon);
field sCov[3]; cross(eCon, bCon, sCov); sCov[0]*= gDet; sCov[1]*= gDet; sCov[2]*= gDet;
field sCon[3]; matMul(gCon, sCov, sCon);
const field d = rh/glfInv, h = rh + GAMMA1*pg, w = h * glf * glf, uem = 0.5*(e2+b2);
const field w1= d*u2/(1.+glf) + GAMMA1*pg*glf*glf; // W' = W-D
//-- Conserved
uD[RH] = d; uD[PG] = w1-pg+uem;
uD[VX] = w*vCov[0]+sCov[0]; uD[VY] = w*vCov[1]+sCov[1]; uD[VZ] = w*vCov[2]+sCov[2];
uD[BX] = bCon[0] ; uD[BY] = bCon[1]; uD[BZ] = bCon[2];
for(unsigned iVar = 0; iVar< FLD_TOT; ++iVar){ uD[iVar]*= gDet;}
//-- Fluxes
f[VX] = w*vCon[dir]*vCov[0]-eCon[dir]*eCov[0]-bCon[dir]*bCov[0];
f[VY] = w*vCon[dir]*vCov[1]-eCon[dir]*eCov[1]-bCon[dir]*bCov[1];
f[VZ] = w*vCon[dir]*vCov[2]-eCon[dir]*eCov[2]-bCon[dir]*bCov[2];
f[VX+k1]+= pg+uem;
f[RH] = d*vCon[dir]; f[PG] = w1*vCon[dir]+sCon[dir]; // Seems odd
f[RH] = gDet*alpha*f[RH] - betai[dir]*uD[RH];
f[VX] = gDet*alpha*f[VX] - betai[dir]*uD[VX];
f[VY] = gDet*alpha*f[VY] - betai[dir]*uD[VY];
f[VZ] = gDet*alpha*f[VZ] - betai[dir]*uD[VZ];
f[PG] = gDet*alpha*f[PG] - betai[dir]*uD[PG];
field tmp[3]; cross(betai, bCon, tmp);
eCov[0] = alpha*eCov[0] + gDet * tmp[0];
eCov[1] = alpha*eCov[1] + gDet * tmp[1];
eCov[2] = alpha*eCov[2] + gDet * tmp[2];
f[BX+k1] = 0.0; f[BY+k2] =-eCov[2+k3]; f[BZ+k3] = eCov[1+k2];
//-- Fast magnetosonic speeds (vCon along direction dir)
const field cs2=GAMMA*pg/h, ca2=1.0-h/(h+sycl::max(b2-e2,0.)), a2=cs2+ca2-cs2*ca2, v2 = u2/(1.0+u2);
const field vf1 = vCon[dir]*(1.0-a2)/(1.0-v2*a2);
const field vf2 = sycl::sqrt( a2*glfInv*glfInv* ( (1.-v2*a2)*gCon[4*dir]-(1.-a2)*vCon[dir]*vCon[dir]))/(1.-v2*a2);
vf[0] = alpha*(vf1+vf2)-betai[dir];
vf[1] = alpha*(vf1-vf2)-betai[dir];
//-- Transverse speeds (vCon)
vt[0] = alpha*vCon[1+k2]-betai[1+k2];
vt[1] = alpha*vCon[2+k3]-betai[2+k3];
#endif // PHYSICS
}
// IMPORTANT: all cell-centered quantities here, take from the v!
void physicalSource(id<1> myId, field_array v, Metric &g, field src[4]){
#if METRIC > CARTESIAN
for(short unsigned i=0; i<4; i++){ src[i] = 0; }; return;
#else
field ssCon[3][3];
//-- Metric components
field betai[3], gCov[9], gCon[9]; g.g3DCov(gCov); g.beta(betai);
field const alpha=g.alpha(), gDet=g.g3DCon(gCon), gDet1=g.gDet1();
//-- Readouts (all cell-centered values, unlike in physicalFlux)
const field rh=v[RH][myId], pg=v[PG][myId];
field bCov[3], bCon[]={v[BX][myId],v[BY][myId],v[BZ][myId]};
matMul(gCov, bCon, bCov); const field b2 = dot(bCon,bCov);
#if PHYSICS==MHD
const field pt=pg+.5*b2, et=rh, vCon[]={v[VX][myId],v[VY][myId],v[VZ][myId]};
field sCon[3] = {rh*vCon[0], rh*vCon[1], rh*vCon[2]};
for(short unsigned i=0; i<3; i++)
for(short unsigned j=0; j<3; j++)
ssCon[i][j] = rh*vCon[i]*vCon[j]-bCon[i]*bCon[j]+pt*gCon[i*3+j];
#elif PHYSICS==GRMHD
field vCov[3], vCon[]={v[VX][myId],v[VY][myId],v[VZ][myId]}; matMul(gCov, vCon, vCov);
const field u2 = dot(vCov,vCon), glfInv = sycl::rsqrt(1+u2), glf = 1.0/glfInv;
for(short unsigned i=0; i<3; i++){ vCon[i]*= glfInv; vCov[i]*= glfInv;}
field eCov[3]; cross (vCon, bCon, eCov); eCov[0]*=-gDet; eCov[1]*=-gDet; eCov[2]*=-gDet;
field eCon[3]; matMul(gCon, eCov, eCon); const field e2 = dot(eCov,eCon);
const field d=rh/glfInv, h=rh+GAMMA1*pg, w=h*glf*glf, uem=0.5*(e2+b2), pt=pg+uem, et=w-pg+uem;
field sCon[3]; cross (eCov, bCov, sCon); //was: s_con=w*v_con+g%det1*cross_product(e_cov,b_cov)
for(short unsigned i=0; i<3; i++){ sCon[i] = gDet1*sCon[i]+ w*vCon[i];}
for(short unsigned i=0; i<3; i++)
for(short unsigned j=0; j<3; j++)
ssCon[i][j] = w*vCon[i]*vCon[j]+pt*g.gCon(i,j)-bCon[i]*bCon[j]-eCon[i]*eCon[j];
#endif // PHYSICS
// WARNING: source tested only for CARTESIAN case (trivial!)
// This was in mod_metric::source. Executed only once, so let's just have it here. -SC
field sum1[3] = {0.0, 0.0, 0.0}, sum2 = 0.0;
for(short unsigned i=0; i<3; i++)
for(short unsigned j=0; j<3; j++)
for(short unsigned iS=0; iS<3; iS++){
sum1[iS]+= 0.5*ssCon[i ][j] *g.dgCov(i,j,iS); // sum1(1:3)=sum1(1:3)+0.5*ss_con(i,j)*dg%cov(i,j,1:3)
sum2 += gCov[i*3+iS]*ssCon[iS][j] *g.dgBeta(i,j); // sum2=sum2+dot_product(g%cov(i,:),ss_con(:,j))*dg%beta(i,j)
}
field sCov[3]; matMul(gCov,sCon,sCov);
// source(1:3)=g%alpha*sum1+matmul(s_cov,dg%beta)-et*dg%alpha
for(short unsigned iS=0; iS<3; iS++){
src[iS] = g.alpha()*sum1[iS] - et*g.dgAlpha(iS);
for(short unsigned k=0; k<3; k++){ src[iS]+= sCov[k] * g.dgBeta(k,iS); }
}
// source(4) =dot_product(g%beta,sum1)+sum2-dot_product(s_con,dg%alpha)
src[3] = sum2;
for(short unsigned k=0; k<3; k++){ src[3]+= betai[k]*sum1[k]-sCon[k]*g.dgAlpha(k); }
// This was located back in physics_source.
for(short unsigned iS=0; iS<4; iS++){ src[iS]*= gDet; }
return;
#endif // METRIC type
}