Peano 4
Loading...
Searching...
No Matches
SSInfall.cpp
Go to the documentation of this file.
2
3#include "Constants.h"
5
8
9
10tarch::logging::Log applications::exahype2::euler::sphericalaccretion::SSInfall::_log( "applications::exahype2::euler::sphericalaccretion::SSInfall" );
11
12
14
15
16#if defined(GPUOffloadingOMP)
17#pragma omp declare target
18#endif
22#if defined(GPUOffloadingOMP)
23#pragma omp end declare target
24#endif
25
26
28 double globalMinTimeStamp,
29 double globalMaxTimeStamp,
30 double globalMinTimeStepSize,
31 double globalMaxTimeStepSize
32){
33 AbstractSSInfall::startTimeStep(globalMinTimeStamp, globalMaxTimeStamp, globalMinTimeStepSize, globalMaxTimeStepSize);
34 constexpr double pi = M_PI;
35 if (AbstractSSInfall::isFirstGridSweepOfTimeStep()){
36 for (int i=0;i<sample_number;i++) {
37 m_tot_copy[i] = global_m_tot[i];
38 global_m_tot[i] = 0;
39 global_cell_tot[i] = 0;
40 m_tot[i] = 0;
41 cell_tot[i] = 0;
42 }
43 double inter_time=_minTimeStamp+getAdmissibleTimeStepSize()/2.0;
44 #ifdef useTable
45 int time_position=0;
46 //logInfo( "startTImeStep()", "current min time=" << globalMinTimeStamp<<" current min timestep=" << globalMinTimeStepSize<<" inter time=" << inter_time);
47 //logInfo( "startTImeStep()", "current max time=" << globalMaxTimeStamp<<" current max timestep=" << globalMaxTimeStepSize);
48 //logInfo( "startTImeStep()", "current time=" << _minTimeStamp<<" current timestep=" << getAdmissibleTimeStepSize());
49
50 if (tarch::la::equals(inter_time,0.0) ){
51 interpolated_a=scale_factor_points[0];
52 //logInfo( "startTImeStep()", "reset scale a to " << interpolated_a << " because t=0.0");
53 }else{
54 for(int i=0; i<table_point_number;i++){
55 if (tarch::la::smallerEquals(inter_time,code_time_points[i]) ){time_position=i; break;}
56 time_position=i;
57 }
58 //logInfo( "startTImeStep()", inter_time << " is between time point " << code_time_points[time_position-1] << " and " << code_time_points[time_position]);
59 interpolated_a=scale_factor_points[time_position-1]*(code_time_points[time_position]-inter_time)/(code_time_points[time_position]-code_time_points[time_position-1])+scale_factor_points[time_position]*(inter_time-code_time_points[time_position-1])/(code_time_points[time_position]-code_time_points[time_position-1]);
60 //logInfo( "startTImeStep()", "interpolated to " << interpolated_a << " from " << scale_factor_points[time_position-1] << " and " << scale_factor_points[time_position]);
61 }
62 #endif
63
64 if (MassCal==1){
65 for(int i=0;i<sample_number;i++){
66 m_tot_copy[i]=std::max(0.0,(4/3)*pi*pow(r_s[0],3)*((rho_0+rho_x[0])/2-1));
67 }
68 for(int i=1;i<sample_number;i++){
69 double m_layer=(4/3)*pi*(pow(r_s[i],3)-pow(r_s[i-1],3))*((rho_x[i]+rho_x[i-1])/2-1);
70 for(int j=i;j<sample_number;j++){
71 m_tot_copy[j]+=std::max(0.0,m_layer);
72 }
73 }
74 }
75 }
76}
77
78
80 AbstractSSInfall::finishTimeStep();
81
82 if (AbstractSSInfall::isLastGridSweepOfTimeStep()){
83 #ifdef Parallel
85 m_tot,
86 global_m_tot,
87 sample_number, MPI_DOUBLE,
88 MPI_SUM,
90 );
92 cell_tot,
93 global_cell_tot,
94 sample_number, MPI_DOUBLE,
95 MPI_SUM,
97 );
98 #else
99 for (int i=0;i<sample_number;i++) {
100 global_m_tot[i]=m_tot[i];
101 global_cell_tot[i]=cell_tot[i];
102 }
103
104 // @Han is this ok?
105 _rMax = r_s[sample_number-1];
106 _rhoMax = rho_x[sample_number-1];
107 _mTotMax = m_tot_copy[sample_number-1];
108 #endif
109 plotGlobalMTot(false);
110 }
111}
112
113
115 std::ostringstream msg;
116 msg << "r_s=(";
117 for (int i=0;i<sample_number;i++) {
118 if (i!=0) msg << ",";
119 msg << r_s[i];
120 }
121 msg << ")";
122 return msg.str();
123}
124
125
127 std::ostringstream msg;
128 if (plotCopy) {
129 msg << "m_tot_copy=(";
130 for (int i=0;i<sample_number;i++) {
131 if (i!=0) msg << ",";
132 msg << m_tot_copy[i];
133 }
134 msg << ")";
135 }
136 else {
137 msg << "global_m_tot=(";
138 for (int i=0;i<sample_number;i++) {
139 if (i!=0) msg << ",";
140 msg << global_m_tot[i];
141 }
142 msg << ")";
143 }
144 return msg.str();
145}
146
147
149 const double * __restrict__ Q, // Q[5+0]
152 double t
153) {
154
156 double radius=volumeX(0)*volumeX(0)+volumeX(1)*volumeX(1)+volumeX(2)*volumeX(2); radius=pow(radius,0.5);
157 if (ReSwi==1){ //radius based
158 constexpr int NumberOfRefinementLayers = 2;
159 double Radius[NumberOfRefinementLayers] = {0.7,0.5};
160 double MaxH[NumberOfRefinementLayers] = {0.02,0.006};
161 //constexpr int NumberOfRefinementLayers = 1;
162 //double Radius[NumberOfRefinementLayers] = {0.5};
163 //double MaxH[NumberOfRefinementLayers] = {0.02};
164 for (int i=0; i<NumberOfRefinementLayers; i++) {
165 if (radius<Radius[i] and tarch::la::max(volumeH)>MaxH[i]) {
167 }
168 }
169 }
170 if (ReSwi==2){ //radius based
171 constexpr int NumberOfRefinementLayers = 5;
172 double Radius[NumberOfRefinementLayers] = {1.5,0.5,0.3,0.1,0.05};
173 double MaxH[NumberOfRefinementLayers] = {0.05,0.02,0.006,0.002,0.0007};
174 for (int i=0; i<NumberOfRefinementLayers; i++) {
175 if (radius<Radius[i] and tarch::la::max(volumeH)>MaxH[i]) {
177 }
178 }
179 }
180 if (ReSwi==3){ //radius based, for test
181 if (radius<0.1) {result=::exahype2::RefinementCommand::Refine;}
182 }
183 if (ReSwi==4){ //dynamic one
184 double t_i=pow(a_i,1.5)*2.057*1e17;
185 double t_used=t;
186 double Mpc=3.086*1e19; //in Km;
187 double t_real=pow( (pow((Mpc/150),(-1.0/3.0))*pow(a_i,(-0.5))-t_used*pow(150,(4.0/3.0))/300/pow(Mpc,(1.0/3.0))),-3);
188 double a_scale=pow((150/Mpc),(2.0/3.0))*pow( (pow((Mpc/150),(-1.0/3.0))*pow(a_i,(-0.5))-t_used*pow(150,(4.0/3.0))/300/pow(Mpc,(1.0/3.0))),-2);
189 double R_ta_code=(a_i/a_scale)*pow((3.0*delta_m)/(4.0*M_PI*tilde_rho_ini),(1.0/3.0))*pow((4.0/3.0/M_PI),(8.0/9.0))*pow((t_real/t_i),(8.0/9.0));
190 double delta_R=std::abs(radius-R_ta_code*0.55);
191 //double delta_R=std::abs(radius-0.0875);
192
193 constexpr int NumberOfRefinementLayers = 3;
194 double distance[NumberOfRefinementLayers+1] = {1.5,0.6,0.3,0};
195 double MaxH[NumberOfRefinementLayers+1] = {0.05,0.018,0.006,0};
196
197 /*if (delta_R>distance[0]) {result=::exahype2::RefinementCommand::Erase;}
198 else if (delta_R<distance[NumberOfRefinementLayers-1]) {result=::exahype2::RefinementCommand::Refine;}
199 else{
200 for (int i=0; i<(NumberOfRefinementLayers-1); i++) {
201 if (delta_R<distance[i] and delta_R>distance[i+1]) {
202 if (tarch::la::max(volumeH)>MaxH[i]) {result=::exahype2::RefinementCommand::Refine;}
203 if (tarch::la::max(volumeH)<MaxH[i+1]) {result=::exahype2::RefinementCommand::Erase;}
204 }
205 }
206 }*/
207 for (int i=0; i<(NumberOfRefinementLayers); i++) {
208 if (delta_R<distance[i] and delta_R>distance[i+1]) {
209 if (tarch::la::max(volumeH)>MaxH[i]) {result=::exahype2::RefinementCommand::Refine;}
210 //if (tarch::la::max(volumeH)<MaxH[i+1]) {result=::exahype2::RefinementCommand::Erase;}
211 }
212 }
213 }
214 if (ReSwi==5){ //erase test
215 double delta_R=radius;
216
217 constexpr int NumberOfRefinementLayers = 1;
218 double distance[NumberOfRefinementLayers] = {0.1};
219 double MaxH[NumberOfRefinementLayers] = {0.01};
220
221 if (delta_R>distance[0]) {result=::exahype2::RefinementCommand::Erase;}
222 else if (delta_R<distance[NumberOfRefinementLayers-1]) {result=::exahype2::RefinementCommand::Refine;}
223 else{
224 for (int i=0; i<(NumberOfRefinementLayers-1); i++) {
225 if (delta_R<distance[i] and delta_R>distance[i+1]) {
226 if (tarch::la::max(volumeH)>MaxH[i]) {result=::exahype2::RefinementCommand::Refine;}
227 if (tarch::la::max(volumeH)<MaxH[i+1]) {result=::exahype2::RefinementCommand::Erase;}
228 }
229 }
230 }
231 }
232 return result;
233}
234
235
236
238 double * __restrict__ Q,
239 const tarch::la::Vector<Dimensions,double>& volumeCentre,
241 bool gridIsConstructed
242) {
243 logDebug( "initialCondition(...)", "init volume at " << volumeCentre << "x" << volumeH << " (grid constructed=" << gridIsConstructed << ")" );
244
245 constexpr double pi = M_PI;
246 double x=volumeCentre(0)-OverDensityCentre(0);
247 double y=volumeCentre(1)-OverDensityCentre(1);
248 double z=volumeCentre(2)-OverDensityCentre(2);
249
250 bool isInTheSphere = ( (x*x+y*y+z*z) < r_ini*r_ini );//overdensity region
251 double r_coor=x*x+y*y+z*z;
252 r_coor=pow(r_coor,0.5);
253
254 //double H_i=2/(3*t_ini); //Hubble constant
255 //double rho_ini=1/(6*pi*G*t_ini*t_ini);
256
257 double rho_ini=tilde_rho_ini;
258 //constexpr double gamma = 5.0/3.0;
259 if (iseed==0) {
260 Q[0] = isInTheSphere ? rho_ini*(1+delta_rho) : rho_ini;
261 } // rho
262 if (iseed==1)
263 {Q[0] = rho_ini;}
264 if (v_scale==0)
265 {Q[1] = 0; Q[2] = 0; Q[3] = 0;} // velocities
266 else
267 {if (r_coor<r_point)
268 {Q[1]=-v_scale*x*1.5*Omega_m*a_i*delta_m/4/pi/pow(r_point,3)*Q[0];
269 Q[2]=-v_scale*y*1.5*Omega_m*a_i*delta_m/4/pi/pow(r_point,3)*Q[0];
270 Q[3]=-v_scale*z*1.5*Omega_m*a_i*delta_m/4/pi/pow(r_point,3)*Q[0];}
271 else
272 {Q[1]=-v_scale*x*1.5*Omega_m*a_i*delta_m/4/pi/pow(r_coor,3)*Q[0];
273 Q[2]=-v_scale*y*1.5*Omega_m*a_i*delta_m/4/pi/pow(r_coor,3)*Q[0];
274 Q[3]=-v_scale*z*1.5*Omega_m*a_i*delta_m/4/pi/pow(r_coor,3)*Q[0];}
275 }
276 Q[4] = 0.5*(Q[1]*Q[1]+Q[2]*Q[2]+Q[3]*Q[3])/Q[0]+tilde_P_ini/(gamma-1); // inner energy
277
278 for (int i=NumberOfUnknowns; i<NumberOfUnknowns+NumberOfAuxiliaryVariables; i++) {
279 Q[i] = 0.0;
280 }
281
282 const double irho = 1./Q[0];
283 #if Dimensions==3
284 const double p = (gamma-1) * (Q[4] - 0.5*irho*(Q[1]*Q[1]+Q[2]*Q[2]+Q[3]*Q[3]));
285 #else
286 const double p = (gamma-1) * (Q[4] - 0.5*irho*(Q[1]*Q[1]+Q[2]*Q[2]));
287 #endif
288
289 Q[0] = rho_ini;
290 nonCriticalAssertion6( p>=0.0, Q[0], Q[1], Q[2], Q[3], Q[4], volumeH );
291 assertion5( Q[0]>1e-12, Q[1], Q[2], Q[3], Q[4], volumeH );
292 // initial conditions
293/*
294 Q[0] = isInTheSphere ? rho_ini*(1+delta_rho) : rho_ini; // rho
295 Q[1] = Q[0]*H_i*x; // velocities
296 Q[2] = Q[0]*H_i*y;
297 Q[3] = Q[0]*H_i*z;
298 Q[4] = 0.5*(Q[1]*Q[1]+Q[2]*Q[2]+Q[3]*Q[3])/Q[0]+initial_internal_energy; // inner energy
299*/
300 tarch::la::Vector<Dimensions,double> circleCentre = {0,0,0};
301/*
302 // initial conditions
303 bool isInTheCentre = ( tarch::la::norm2( volumeCentre-circleCentre ) < 0.2 );
304 Q[0] = 0.1; // rho
305 Q[1] = 0; // velocities
306 Q[2] = 0;
307 Q[3] = 0;
308 Q[4] = isInTheCentre ? 1.0 : 0.0; // inner energy
309*/
310 if ( Q[0]<1e-12 ){
311 ::tarch::triggerNonCriticalAssertion( __FILE__, __LINE__, "Q[0]>0", "density become 0 at volume: ["+std::to_string(volumeCentre(0))+", "+std::to_string(volumeCentre(1))+", "+std::to_string(volumeCentre(2))+"] Q[] array: "+std::to_string(Q[0])+" "+std::to_string(Q[1])+" "+std::to_string(Q[2])+" "+std::to_string(Q[3])+" "+std::to_string(Q[4])+".");
312 }
313
314 for (int i=5; i<NumberOfUnknowns+NumberOfAuxiliaryVariables; i++) {
315 Q[i] = 0.0;
316 }
317}
318
319
321 const double rho_ini=tilde_rho_ini;
322 if (iseed==0) {//not using
323 if (r_coor<r_ini) {
324 return rho_ini*delta_rho*pow(r_coor,3)/3;
325 }
326 else {
327 return rho_ini*delta_rho*pow(r_ini,3)/3;
328 }
329 }
330 if (iseed==1) {//tophat
331 if (r_coor<r_point){
332 return delta_m*pow(r_coor/r_point,3)/4/tarch::la::PI;
333 }
334 else {
335 return delta_m/4/tarch::la::PI;
336 }
337 }
338 return 0.0;
339}
340
341
342double applications::exahype2::euler::sphericalaccretion::SSInfall::getForceDensityNorm(double r_coor,double m_in, double t, double density, double a_input) {
343 double a=0.0287*pow((-t/11.8+0.1694*pow(a_i,-0.5)),-2);//when code time ~ 2*(a_i^(-0.5)-1), a~1
344 #ifdef useTable
345 a=a_input; //if we use a table to interpolate a (for LCDM), we replace it here.
346 #endif
347 double force_density_norm=density*G*m_in/pow(r_coor,3)*Omega_m*a*1.5;
348 #ifdef DGP
349 //double y3_over_Delta=(4.0/3.0)*tarch::la::PI*pow(r_coor,3)/delta_m;
350 double l3_over_tildeM=(4.0/3.0)*tarch::la::PI*pow(r_coor,3)/m_in;
351 double modifier=0.0;
352 if (tarch::la::equals(DGP_zeta,0.0)){modifier=1.0/3.0/DGP_beta;}
353 else{modifier=(27*DGP_beta/16/DGP_zeta/DGP_zeta)*l3_over_tildeM*(pow(1+(32*DGP_zeta*DGP_zeta/81/DGP_beta/DGP_beta/l3_over_tildeM),0.5)-1);}
354 //logInfo( "source()", "modifier" << modifier );
355 force_density_norm=force_density_norm*(1.0+modifier);
356 #endif
357 return force_density_norm;
358}
359
360
362 const double * __restrict__ Q, // Q[5+0]
365 double t,
366 double dt,
367 double * __restrict__ S // S[5
368) {
369 logTraceInWith4Arguments( "sourceTerm(...)", volumeX, volumeH, t, dt );
370
371 nonCriticalAssertion4( not std::isnan(Q[0]), volumeX, volumeH, t, dt );
372
373 assertion4( not std::isnan(Q[0]), volumeX, volumeH, t, dt );
374 assertion4( not std::isnan(Q[1]), volumeX, volumeH, t, dt );
375 assertion4( not std::isnan(Q[2]), volumeX, volumeH, t, dt );
376 assertion4( not std::isnan(Q[3]), volumeX, volumeH, t, dt );
377 assertion4( not std::isnan(Q[4]), volumeX, volumeH, t, dt );
378 assertion4( not std::isinf(Q[0]), volumeX, volumeH, t, dt );
379 assertion4( not std::isinf(Q[1]), volumeX, volumeH, t, dt );
380 assertion4( not std::isinf(Q[2]), volumeX, volumeH, t, dt );
381 assertion4( not std::isinf(Q[3]), volumeX, volumeH, t, dt );
382 assertion4( not std::isinf(Q[4]), volumeX, volumeH, t, dt );
383
385 not std::isnan(Q[0]) or t<0.5,
386 volumeX, volumeH, t, dt
387 );
388
389 double x=volumeX(0)-OverDensityCentre(0);
390 double y=volumeX(1)-OverDensityCentre(1);
391 double z=volumeX(2)-OverDensityCentre(2);
392
393 double r_coor=x*x+y*y+z*z;
394 r_coor=pow(r_coor,0.5);
395 double m_in=0;
396
397 if (tarch::la::equals(t,0)){//we know the mass distri at the beginning
398 m_in = getInitialMIn( r_coor );
399 }
400 else {
401 if (iseed==0){
402 m_in=mass_interpolate(r_coor,MassCal)/4.0/tarch::la::PI; //remove the overall 4\pi coefficient.
403 nonCriticalAssertion4( not std::isnan(m_in), volumeX, volumeH, t, dt );
404 nonCriticalAssertion4( m_in<1e100, volumeX, volumeH, t, dt );
405 }
406 if (iseed==1){
407 if (r_coor<r_point){
408 m_in=(mass_interpolate(r_coor, MassCal)+delta_m*pow(r_coor/r_point,3))/4.0/tarch::la::PI;
409 nonCriticalAssertion4( not std::isnan(m_in), volumeX, volumeH, t, dt );
410 nonCriticalAssertion4( m_in<1e100, volumeX, volumeH, t, dt );
411 }
412 else {
413 assertion3( isOutsideOfLargestRadius(volumeX, volumeH), volumeX, volumeH, _rMax );
414 m_in=(mass_interpolate(r_coor, MassCal)+delta_m)/4.0/tarch::la::PI;
415 nonCriticalAssertion4( not std::isnan(m_in), volumeX, volumeH, t, dt );
416 nonCriticalAssertion4( m_in<1e100, volumeX, volumeH, t, dt );
417 }
418 }
419 }
420
421 double force_density_norm = getForceDensityNorm(
422 r_coor,
423 m_in,
424 t, Q[0],interpolated_a
425 );
426
427 S[0] = 0; // rho
428 S[1] = -force_density_norm*x; // velocities
429 S[2] = -force_density_norm*y;
430 S[3] = -force_density_norm*z;
431 S[4] = -force_density_norm*(Q[1]*x+Q[2]*y+Q[3]*z)/Q[0];
432
433 for (int i=5; i<NumberOfUnknowns+NumberOfAuxiliaryVariables; i++) {
434 S[i] = 0.0;
435 }
436
437 logTraceOut( "sourceTerm(...)" );
438}
439
440
441#if defined(GPUOffloadingOMP)
442#pragma omp declare target
443#endif
445 const double * __restrict__ Q, // Q[5+0]
448 double t,
449 double dt,
450 double * __restrict__ S, // S[5
451 Offloadable
452) {
453 double r_coor = tarch::la::norm2(volumeX - OverDensityCentre);
454
455 double m_in = (4.0/3.0)*tarch::la::PI*(pow(r_coor,3)-pow(_rMax,3))*(_rMax-1.0) + _mTotMax;
456
457 double force_density_norm = getForceDensityNorm(
458 r_coor,
459 m_in,
460 t, Q[0],1.0
461 );
462
463 double x = volumeX(0)-OverDensityCentre(0);
464 double y = volumeX(1)-OverDensityCentre(1);
465 double z = volumeX(2)-OverDensityCentre(2);
466
467 S[0] = 0.0; // rho
468 S[1] = -force_density_norm*x; // velocities
469 S[2] = -force_density_norm*y;
470 S[3] = -force_density_norm*z;
471 S[4] = -force_density_norm*(Q[1]*x+Q[2]*y+Q[3]*z)/Q[0];
472
473 for (int i=5; i<NumberOfUnknowns+NumberOfAuxiliaryVariables; i++) {
474 S[i] = 0.0;
475 }
476}
477#if defined(GPUOffloadingOMP)
478#pragma omp end declare target
479#endif
480
481
483 const tarch::la::Vector<Dimensions,double>& patchCentre,
485) const {
486 double r_coor = tarch::la::norm2(patchCentre - OverDensityCentre);
487 double safetyDistance = std::sqrt(3) * tarch::la::max(patchH) / 2.0;
488 return r_coor + safetyDistance > _rMax;
489}
490
491
493 const tarch::la::Vector<Dimensions,double>& patchCentre,
495 double t,
496 double dt
497) const {
498 return tarch::la::greater(t,0.0) and isOutsideOfLargestRadius(patchCentre,patchH);
499}
500
501
503 const double * __restrict__ Qinside, // Qinside[5+0]
504 double * __restrict__ Qoutside, // Qoutside[5+0]
505 const tarch::la::Vector<Dimensions,double>& faceCentre,
507 double t,
508 int normal
509) {
510 logTraceInWith4Arguments( "boundaryConditions(...)", faceCentre, volumeH, t, normal );
511
512 nonCriticalAssertion4( Qinside[0]==Qinside[0], faceCentre, volumeH, t, normal );
513 nonCriticalAssertion4( Qinside[1]==Qinside[1], faceCentre, volumeH, t, normal );
514 nonCriticalAssertion4( Qinside[2]==Qinside[2], faceCentre, volumeH, t, normal );
515 nonCriticalAssertion4( Qinside[3]==Qinside[3], faceCentre, volumeH, t, normal );
516 nonCriticalAssertion4( Qinside[4]==Qinside[4], faceCentre, volumeH, t, normal );
517
518 nonCriticalAssertion5( Qinside[0]>1e-12, faceCentre, volumeH, t, normal, Qinside[0] );
519/*
520 assertion5( Qinside[0]>1e-12, faceCentre, volumeH, t, normal, Qinside[0] );
521 if ( Qinside[0]<1e-12 ){
522 ::tarch::triggerNonCriticalAssertion( __FILE__, __LINE__, "Qinside[0]>0", "density become 0 at face: ["+std::to_string(faceCentre(0))+", "+std::to_string(faceCentre(1))+", "+std::to_string(faceCentre(2))+"] Q[] array: "+std::to_string(Qinside[0])+" "+std::to_string(Qinside[1])+" "+std::to_string(Qinside[2])+" "+std::to_string(Qinside[3])+" "+std::to_string(Qinside[4])+".");
523 }
524*/
525
526 if (extrapolate_bc==0){
527 Qoutside[0] = Qinside[0];
528 Qoutside[1] = Qinside[1];
529 Qoutside[2] = Qinside[2];
530 Qoutside[3] = Qinside[3];
531 Qoutside[4] = Qinside[4];
532 }
533 else if (extrapolate_bc==1)
534 {
535 Qoutside[0] = Qinside[0];
536 Qoutside[4] = Qinside[4];
537 for (int i=0; i<5; i++){
538 if (normal<3) {
539 Qoutside[i]=Qinside[i]+Qinside[5+i*3+normal]*(-volumeH(normal));
540 }
541 else if (normal>=3) {
542 Qoutside[i]=Qinside[i]+Qinside[5+i*3+normal-3]*(volumeH(normal-3));
543 }
544 }
545 for (int i=NumberOfUnknowns; i<NumberOfUnknowns+NumberOfAuxiliaryVariables; i++) {
546 Qoutside[i]=0;
547 }
548 }
549 else if (extrapolate_bc==2)
550 {
552 not std::isnan(Qinside[0]),
553 normal, faceCentre, Qinside[0], Qinside[1], Qinside[2], Qinside[3], Qinside[4]
554 );
555 for (int i=NumberOfUnknowns; i<NumberOfUnknowns+NumberOfAuxiliaryVariables; i++) {
556 Qoutside[i]=0;
557 }
558 for (int i=0; i<5; i++){
559 if (normal<3) {
560 Qoutside[i]=Qinside[i]+Qinside[5+i*3+normal]*(-volumeH(normal));
561 }
562 else if (normal>=3) {
563 Qoutside[i]=Qinside[i]+Qinside[5+i*3+normal-3]*(volumeH(normal-3));
564 }
565 }
566 if (Qinside[0]<tilde_rho_ini){
567 for (int i=0; i<5; i++){Qoutside[i] = Qinside[i];}
568 }
569 }
570 //add more constraints here
571 //for (int j=1; j<=3; j++){
572 // if (Qoutside[j]*Qinside[j]<0) {Qoutside[j]=0;}
573 //}
574 const double p = (gamma-1) * (Qoutside[4] - 0.5*(Qoutside[1]*Qoutside[1]+Qoutside[2]*Qoutside[2]+Qoutside[3]*Qoutside[3])/Qoutside[0]);
575 if (p<0){
576 Qoutside[4]=0.5*(Qoutside[1]*Qoutside[1]+Qoutside[2]*Qoutside[2]+Qoutside[3]*Qoutside[3])/Qoutside[0]+1e-10;
577 }
579 not std::isnan(Qinside[0]),
580 normal, faceCentre, Qinside[0], Qinside[1], Qinside[2], Qinside[3], Qinside[4]
581 );
582
583 for (int i=5; i<NumberOfUnknowns+NumberOfAuxiliaryVariables; i++) {
584 Qoutside[i]=0;
585 }
586
587 logTraceOut( "boundaryConditions(...)" );
588}
589
590
592 const double * __restrict__ Q, // Q[5+0],
593 const tarch::la::Vector<Dimensions,double>& faceCentre,
595 double t,
596 double dt,
597 int normal
598) {
599 assertion(normal>=0);
600 assertion(normal<Dimensions);
601 assertion( Q[0]>0.0 );
602
603 assertion4( not std::isnan(Q[0]), faceCentre, volumeH, t, normal );
604 assertion4( not std::isnan(Q[1]), faceCentre, volumeH, t, normal );
605 assertion4( not std::isnan(Q[2]), faceCentre, volumeH, t, normal );
606 assertion4( not std::isnan(Q[3]), faceCentre, volumeH, t, normal );
607 assertion4( not std::isnan(Q[4]), faceCentre, volumeH, t, normal );
608 assertion4( not std::isinf(Q[0]), faceCentre, volumeH, t, normal );
609 assertion4( not std::isinf(Q[1]), faceCentre, volumeH, t, normal );
610 assertion4( not std::isinf(Q[2]), faceCentre, volumeH, t, normal );
611 assertion4( not std::isinf(Q[3]), faceCentre, volumeH, t, normal );
612 assertion4( not std::isinf(Q[4]), faceCentre, volumeH, t, normal );
613
614 double result = maxEigenvalue(Q,faceCentre,volumeH,t,dt,normal,Offloadable::Yes);
615
617 result<10000,
618 result, Q[0], Q[1], Q[2], Q[3], Q[4], faceCentre, volumeH, t, normal
619 );
620 //if (t>1.0 and result<1e-8) {
621 // logInfo( "maxEigenvalue()", "max eigenvalue too small: " << result );
622 // std::abort();
623 //}
624 return result;
625}
626
627
628#if defined(GPUOffloadingOMP)
629#pragma omp declare target
630#endif
632 const double * __restrict__ Q, // Q[5+0],
633 const tarch::la::Vector<Dimensions,double>& faceCentre,
635 double t,
636 double dt,
637 int normal,
638 Offloadable
639) {
640 //constexpr double gamma = 5.0/3.0;
641 const double irho = 1./Q[0];
642 #if Dimensions==3
643 double p = (gamma-1) * (Q[4] - 0.5*irho*(Q[1]*Q[1]+Q[2]*Q[2]+Q[3]*Q[3]));
644 #else
645 double p = (gamma-1) * (Q[4] - 0.5*irho*(Q[1]*Q[1]+Q[2]*Q[2]));
646 #endif
647
648 if ( p<0 ) {
649 p=0;
650 //::tarch::triggerNonCriticalAssertion( __FILE__, __LINE__, "p>=0", "negative pressure "+std::to_string(p)+" detected at t=" + std::to_string(t) + " at face position ["+std::to_string(faceCentre(0))+", "+std::to_string(faceCentre(1))+", "+std::to_string(faceCentre(2))+"] Q[] array: "+std::to_string(Q[0])+" "+std::to_string(Q[1])+" "+std::to_string(Q[2])+" "+std::to_string(Q[3])+" "+std::to_string(Q[4])+".");
651 }
652
653 // nonCriticalAssertion9( p>=0.0, Q[0], Q[1], Q[2], Q[3], Q[4], faceCentre, volumeH, t, normal );
654 const double c = std::sqrt(gamma * p * irho);
655
656 const double u_n = Q[normal + 1] * irho;
657 double result = std::max( std::abs(u_n - c), std::abs(u_n + c)); //result=1;
658 //nonCriticalAssertion14( result>0.0, result, p, u_n, irho, c, Q[0], Q[1], Q[2], Q[3], Q[4], faceCentre, volumeH, t, normal );
659 result=result*(1+C_1*exp(-C_2*t));
660
661 return result;
662}
663#if defined(GPUOffloadingOMP)
664#pragma omp end declare target
665#endif
666
667
669 const double * __restrict__ Q, // Q[5+0],
670 const tarch::la::Vector<Dimensions,double>& faceCentre,
672 double t,
673 double dt,
674 int normal,
675 double * __restrict__ F // F[5]
676) {
677 logTraceInWith4Arguments( "flux(...)", faceCentre, volumeH, t, normal );
678
679 assertion4( Q[0]>0.0, faceCentre, volumeH, t, normal );
680 assertion4( not std::isnan(Q[0]), faceCentre, volumeH, t, normal );
681 assertion4( not std::isnan(Q[1]), faceCentre, volumeH, t, normal );
682 assertion4( not std::isnan(Q[2]), faceCentre, volumeH, t, normal );
683 assertion4( not std::isnan(Q[3]), faceCentre, volumeH, t, normal );
684 assertion4( not std::isnan(Q[4]), faceCentre, volumeH, t, normal );
685 assertion4( not std::isinf(Q[0]), faceCentre, volumeH, t, normal );
686 assertion4( not std::isinf(Q[1]), faceCentre, volumeH, t, normal );
687 assertion4( not std::isinf(Q[2]), faceCentre, volumeH, t, normal );
688 assertion4( not std::isinf(Q[3]), faceCentre, volumeH, t, normal );
689 assertion4( not std::isinf(Q[4]), faceCentre, volumeH, t, normal );
690
691 flux(Q,faceCentre,volumeH,t,dt,normal,F,Offloadable::Yes);
692
693 assertion4( not std::isnan(F[0]), faceCentre, volumeH, t, normal );
694 assertion4( not std::isnan(F[1]), faceCentre, volumeH, t, normal );
695 assertion4( not std::isnan(F[2]), faceCentre, volumeH, t, normal );
696 assertion4( not std::isnan(F[3]), faceCentre, volumeH, t, normal );
697 assertion4( not std::isnan(F[4]), faceCentre, volumeH, t, normal );
698 assertion4( not std::isinf(F[0]), faceCentre, volumeH, t, normal );
699 assertion4( not std::isinf(F[1]), faceCentre, volumeH, t, normal );
700 assertion4( not std::isinf(F[2]), faceCentre, volumeH, t, normal );
701 assertion4( not std::isinf(F[3]), faceCentre, volumeH, t, normal );
702 assertion4( not std::isinf(F[4]), faceCentre, volumeH, t, normal );
703
704 logTraceOutWith4Arguments( "flux(...)", faceCentre, volumeH, t, normal );
705}
706
707
708#if defined(GPUOffloadingOMP)
709#pragma omp declare target
710#endif
712 const double * __restrict__ Q, // Q[5+0],
713 const tarch::la::Vector<Dimensions,double>& faceCentre,
715 double t,
716 double dt,
717 int normal,
718 double * __restrict__ F, // F[5]
719 Offloadable
720) {
721 //constexpr double gamma = 5.0/3.0;
722 const double irho = 1./Q[0];
723 #if Dimensions==3
724 const double p = (gamma-1) * (Q[4] - 0.5*irho*(Q[1]*Q[1]+Q[2]*Q[2]+Q[3]*Q[3]));
725 #else
726 const double p = (gamma-1) * (Q[4] - 0.5*irho*(Q[1]*Q[1]+Q[2]*Q[2]));
727 #endif
728
729 const double coeff = irho*Q[normal+1];
730 F[0] = coeff*Q[0];
731 F[1] = coeff*Q[1];
732 F[2] = coeff*Q[2];
733 F[3] = coeff*Q[3];
734 F[4] = coeff*Q[4];
735 F[normal+1] += p;
736 F[4] += coeff*p;
737}
738#if defined(GPUOffloadingOMP)
739#pragma omp end declare target
740#endif
741
742
744 const double r_coor,
745 const double rho,
746 const double size
747) {
748 static tarch::multicore::BooleanSemaphore _mySemaphore;
749
750 //double m=(rho-1)*pow(size,3); //notice here we use overdensity
751 double m=(rho-1)*pow(size,3);
752 if (m<0){m=0.0;}
753 //m=1;
754 tarch::multicore::Lock myLock( _mySemaphore );
755 for (int i=0;i<sample_number;i++){
756 if ((r_coor+size/2)<r_s[i]) {
757 m_tot[i]+=m;
758 cell_tot[i]+=1;
759 global_m_tot[i]+=m;
760 global_cell_tot[i]+=1;
761 }
762 else if ((r_coor-size/2)>r_s[i]) {
763 m_tot[i]+=0;
764 }
765 else {
766 m_tot[i]+=m*std::max(0.0,pow((r_s[i]-r_coor+size/2),3))/pow(size,3);cell_tot[i]+=1;
767 global_m_tot[i]+=m*std::max(0.0,pow((r_s[i]-r_coor+size/2),3))/pow(size,3); global_cell_tot[i]+=1;
768 }
769 }
770}
771
772
774 const double r_coor,
775 const int MassCal
776) {
777 double a,b;
778 double m_a,m_b;
779 double m_result;
780
781 if (MassCal==0){ //which means we use cell counting
782 bool IsCenter=false; // means that we are right in the centre
783 bool IsOutSkirt=false; // means that we are outside of the accreditation area
784
785 if ( tarch::la::smallerEquals( r_coor,r_s[0] ) ) {
786 a=0; b=r_s[0];
787 m_a=0;
788 m_b=m_tot_copy[0];
789 IsCenter=true;
790 }
791 else if ( tarch::la::greater( r_coor, r_s[sample_number-1] ) ) {
792 a=r_s[sample_number-2]; b=r_s[sample_number-1];
793 m_a=m_tot_copy[sample_number-2];
794 m_b=m_tot_copy[sample_number-1];
795 IsOutSkirt=true;
796 }
797 else{
798 for (int i=1;i<sample_number;i++){
799 if ( tarch::la::greater(r_coor, r_s[i-1]) and tarch::la::smallerEquals(r_coor,r_s[i])) {
800 a=r_s[i-1]; b=r_s[i];
801 m_a=m_tot_copy[i-1];
802 m_b=m_tot_copy[i];
803 }
804 }
805 }
806
807 if (IsCenter){
808 m_result=m_b*pow((r_coor),3)/pow(b,3);
809 assertion10( m_result>=0.0, m_result, r_coor, MassCal, m_a, m_b, a, b, IsOutSkirt, plotGlobalMTot(false), plotGlobalMTot(true) );
810 assertion10( m_result<1e100, m_result, r_coor, MassCal, m_a, m_b, a, b, IsOutSkirt, plotGlobalMTot(false), plotGlobalMTot(true) );
811 }
812 else if (IsOutSkirt){
813 double vol_tem=(4.0/3.0)*tarch::la::PI*(pow(b,3)-pow(a,3));
814 double rho_tem=(m_b-m_a)/vol_tem;
815 double vol_out=(4.0/3.0)*tarch::la::PI*(pow(r_coor,3)-pow(b,3));
816 m_result=m_b+rho_tem*vol_out;
817 assertion14( m_result>=0.0, m_result, r_coor, MassCal, m_a, m_b, a, b, IsOutSkirt, IsCenter, plotRs(), plotGlobalMTot(false), plotGlobalMTot(true), rho_tem, vol_out );
818 assertion14( m_result<1e100, m_result, r_coor, MassCal, m_a, m_b, a, b, IsOutSkirt, IsCenter, plotRs(), plotGlobalMTot(false), plotGlobalMTot(true), rho_tem, vol_out );
819 }
820 else { //linear interpolation
821 m_result=m_a*(b-r_coor)/(b-a)+m_b*(r_coor-a)/(b-a);
822 //try to use a more precise mass calculation scheme here
823 /*double vol_tem=(4/3)*pi*(pow(b,3)-pow(a,3));
824 double rho_tem=(m_b-m_a)/vol_tem;
825 double vol_in=(4/3)*pi*(pow(r_coor,3)-pow(a,3));
826 m_result=m_a+rho_tem*vol_in;*/
827 //if (not m_b==0){
828 assertion11( m_result>=0.0, m_result, r_coor, MassCal, m_a, m_b, a, b, IsOutSkirt, plotGlobalMTot(false), plotGlobalMTot(true), plotRs() );
829 assertion11( m_result<1e100, m_result, r_coor, MassCal, m_a, m_b, a, b, IsOutSkirt, plotGlobalMTot(false), plotGlobalMTot(true), plotRs() );
830 }
831 }
832 else if (MassCal==1) { //which means we use rho interpolation
833 if (r_coor<r_s[0]) {
834 double rho_currentpos=rho_0+(rho_x[0]-rho_0)*r_coor/r_s[0];
835 m_result=(4.0/3.0)*tarch::la::PI*pow(r_coor,3)*((rho_0+rho_currentpos)/2-1);
836 assertion9( m_result>=0.0, m_result, r_coor, MassCal, m_a, m_b, a, b, plotGlobalMTot(false), plotGlobalMTot(true) );
837 assertion9( m_result<1e100, m_result, r_coor, MassCal, m_a, m_b, a, b, plotGlobalMTot(false), plotGlobalMTot(true) );
838 }
839 else{
840 for (int i=1;i<sample_number;i++) {
841 if ((r_coor>r_s[i-1]) and (r_coor<r_s[i])) {
842 double rho_currentpos=rho_x[i-1]*(r_s[i]-r_coor)/(r_s[i]-r_s[i-1])+rho_x[i]*(r_coor-r_s[i-1])/(r_s[i]-r_s[i-1]);
843 m_result=(4.0/3.0)*tarch::la::PI*(pow(r_coor,3)-pow(r_s[i-1],3))*((rho_x[i-1]+rho_currentpos)/2-1);
844 m_result+=m_tot_copy[i-1];
845 assertion10( m_result>=0.0, m_result, r_coor, MassCal, m_a, m_b, a, b, i, plotGlobalMTot(false), plotGlobalMTot(true) );
846 assertion10( m_result<1e100, m_result, r_coor, MassCal, m_a, m_b, a, b, i, plotGlobalMTot(false), plotGlobalMTot(true) );
847 }
848 }
849 }
850 if (r_coor>r_s[sample_number-1]){
851 m_result=(4.0/3.0)*tarch::la::PI*(pow(r_coor,3)-pow(r_s[sample_number-1],3))*(rho_x[sample_number-1]-1);
852 m_result+=m_tot_copy[sample_number-1];
853 assertion10( m_result>=0.0, m_result, r_coor, MassCal, m_a, m_b, a, b, r_s[sample_number-1], plotGlobalMTot(false), plotGlobalMTot(true) );
854 assertion10( m_result<1e100, m_result, r_coor, MassCal, m_a, m_b, a, b, r_s[sample_number-1], plotGlobalMTot(false), plotGlobalMTot(true) );
855 }
856 }
857
858 return m_result;
859}
860
#define assertion10(expr, param0, param1, param2, param3, param4, param5, param6, param7, param8, param9)
#define assertion4(expr, param0, param1, param2, param3)
#define assertion3(expr, param0, param1, param2)
#define assertion14(expr, param0, param1, param2, param3, param4, param5, param6, param7, param8, param9, param10, param11, param12, param13)
#define assertion9(expr, param0, param1, param2, param3, param4, param5, param6, param7, param8)
#define assertion(expr)
#define assertion11(expr, param0, param1, param2, param3, param4, param5, param6, param7, param8, param9, param10)
#define assertion5(expr, param0, param1, param2, param3, param4)
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
#define logDebug(methodName, logMacroMessageStream)
Definition Log.h:50
#define logTraceOut(methodName)
Definition Log.h:379
#define logTraceInWith4Arguments(methodName, argument0, argument1, argument2, argument3)
Definition Log.h:373
#define logTraceOutWith4Arguments(methodName, argument0, argument1, argument2, argument3)
Definition Log.h:383
#define nonCriticalAssertion10(expr, param0, param1, param2, param3, param4, param5, param6, param7, param8, param9)
#define nonCriticalAssertion5(expr, param0, param1, param2, param3, param4)
#define nonCriticalAssertion7(expr, param0, param1, param2, param3, param4, param5, param6)
#define nonCriticalAssertion4(expr, param0, param1, param2, param3)
#define nonCriticalAssertion6(expr, param0, param1, param2, param3, param4, param5)
_minTimeStamp(0.0)
static double _rMax
This is a straight 1:1 copy of r_s[sample_number-1].
Definition SSInfall.h:27
virtual double maxEigenvalue(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, double dt, int normal) override
See flux() for a description why there are two maxEigenvalue() versions.
Definition SSInfall.cpp:591
static double getForceDensityNorm(double r_coord, double m_in, double t, double density, double a_input)
Definition SSInfall.cpp:342
bool isOutsideOfLargestRadius(const tarch::la::Vector< Dimensions, double > &patchCentre, const tarch::la::Vector< Dimensions, double > &patchH) const
This is a very conservative routine.
Definition SSInfall.cpp:482
virtual void boundaryConditions(const double *__restrict__ Qinside, double *__restrict__ Qoutside, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, int normal) override
Definition SSInfall.cpp:502
static double _rhoMax
This is a straight 1:1 copy of rho_x[sample_number-1].
Definition SSInfall.h:33
void add_mass(const double r_coor, const double rho, const double size)
Definition SSInfall.cpp:743
virtual void flux(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, double dt, int normal, double *__restrict__ F) override
The flux in the SSInfall does not depend on any variable of the solver.
Definition SSInfall.cpp:668
virtual bool patchCanUseStatelessPDETerms(const tarch::la::Vector< Dimensions, double > &patchCentre, const tarch::la::Vector< Dimensions, double > &patchH, double t, double dt) const
Within the mass accumulation radius, I return false.
Definition SSInfall.cpp:492
double getInitialMIn(double r_coord) const
Returns m_in for time t=0.
Definition SSInfall.cpp:320
static double _mTotMax
From m_tot_copy[sample_number-1].
Definition SSInfall.h:38
double mass_interpolate(const double r_coor, const int MassCal)
Definition SSInfall.cpp:773
void initialCondition(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, bool gridIsConstructed) override
Definition SSInfall.cpp:237
static const tarch::la::Vector< Dimensions, double > OverDensityCentre
Definition SSInfall.h:13
::exahype2::RefinementCommand refinementCriterion(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t) override
Definition SSInfall.cpp:148
void startTimeStep(double globalMinTimeStamp, double globalMaxTimeStamp, double globalMinTimeStepSize, double globalMaxTimeStepSize) override
Definition SSInfall.cpp:27
void sourceTerm(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, double dt, double *__restrict__ S) override
There are two variants of the source term.
Definition SSInfall.cpp:361
Log Device.
Definition Log.h:516
static Rank & getInstance()
This operation returns the singleton instance.
Definition Rank.cpp:538
void allReduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, std::function< void()> waitor=[]() -> void {})
Definition Rank.cpp:279
Create a lock around a boolean semaphore region.
Definition Lock.h:19
virtual void receiveDanglingMessages() override
Answer to MPI Messages.
static ServiceRepository & getInstance()
static void flux(const double *__restrict__ Q, int normal, double *__restrict__ F) InlineMethod
static double maxEigenvalue(const double *__restrict__ Q, int normal) InlineMethod
bool greater(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
double max(double a, double b, double c)
I need the maximum of three values all the time, to I decided to write a function for this.
Definition Scalar.cpp:8
bool smallerEquals(double lhs, double rhs, double tolerance=NUMERICAL_ZERO_DIFFERENCE)
constexpr double PI
Definition Scalar.h:12
bool equals(const Matrix< Rows, Cols, Scalar > &lhs, const Matrix< Rows, Cols, Scalar > &rhs, const Scalar &tolerance=NUMERICAL_ZERO_DIFFERENCE)
Compares to matrices on equality by means of a numerical accuracy.
Scalar norm2(const Vector< Size, Scalar > &vector)
Computes the 2-norm of the vector, i.e.
void triggerNonCriticalAssertion(std::string file, int line, std::string expression, std::string parameterValuePairs)
Trigger a non-critical assertion.
Simple vector class.
Definition Vector.h:134