HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
bounds.c
Go to the documentation of this file.
1 
2 
10 #include "decs.h"
11 
12 int bound_x1dn_radbeamflatinflow(
13  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
14  int *inboundloop,
15  int *outboundloop,
16  int *innormalloop,
17  int *outnormalloop,
18  int (*inoutlohi)[NUMUPDOWN][NDIM],
19  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
20  int *dosetbc,
21  int enerregion,
22  int *localenerpos
23  );
24 
25 int bound_radbeam2dbeaminflow(int dir,
26  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
27  int *inboundloop,
28  int *outboundloop,
29  int *innormalloop,
30  int *outnormalloop,
31  int (*inoutlohi)[NUMUPDOWN][NDIM],
32  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
33  int *dosetbc,
34  int enerregion,
35  int *localenerpos
36  );
37 
38 int bound_radbeam2dflowinflow(int dir,
39  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
40  int *inboundloop,
41  int *outboundloop,
42  int *innormalloop,
43  int *outnormalloop,
44  int (*inoutlohi)[NUMUPDOWN][NDIM],
45  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
46  int *dosetbc,
47  int enerregion,
48  int *localenerpos
49  );
50 
51 
52 int bound_radshadowinflow(int dir,
53  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
54  int *inboundloop,
55  int *outboundloop,
56  int *innormalloop,
57  int *outnormalloop,
58  int (*inoutlohi)[NUMUPDOWN][NDIM],
59  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
60  int *dosetbc,
61  int enerregion,
62  int *localenerpos
63  );
64 
65 int bound_radatmbeaminflow(int dir,
66  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
67  int *inboundloop,
68  int *outboundloop,
69  int *innormalloop,
70  int *outnormalloop,
71  int (*inoutlohi)[NUMUPDOWN][NDIM],
72  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
73  int *dosetbc,
74  int enerregion,
75  int *localenerpos
76  );
77 
78 int bound_radwallinflow(int dir,
79  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
80  int *inboundloop,
81  int *outboundloop,
82  int *innormalloop,
83  int *outnormalloop,
84  int (*inoutlohi)[NUMUPDOWN][NDIM],
85  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
86  int *dosetbc,
87  int enerregion,
88  int *localenerpos
89  );
90 
91 
92 int bound_radbondiinflow(int dir,
93  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
94  int *inboundloop,
95  int *outboundloop,
96  int *innormalloop,
97  int *outnormalloop,
98  int (*inoutlohi)[NUMUPDOWN][NDIM],
99  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
100  int *dosetbc,
101  int enerregion,
102  int *localenerpos
103  );
104 
105 int bound_raddot(
106  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
107  int *inboundloop,
108  int *outboundloop,
109  int *innormalloop,
110  int *outnormalloop,
111  int (*inoutlohi)[NUMUPDOWN][NDIM],
112  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
113  int *dosetbc,
114  int enerregion,
115  int *localenerpos
116  );
117 
118 int bound_radnt(int dir,
119  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
120  int *inboundloop,
121  int *outboundloop,
122  int *innormalloop,
123  int *outnormalloop,
124  int (*inoutlohi)[NUMUPDOWN][NDIM],
125  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
126  int *dosetbc,
127  int enerregion,
128  int *localenerpos
129  );
130 
131 
132 int bound_x1dn_cylaxis(
133  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
134  int *inboundloop,
135  int *outboundloop,
136  int *innormalloop,
137  int *outnormalloop,
138  int (*inoutlohi)[NUMUPDOWN][NDIM],
139  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
140  int *dosetbc,
141  int enerregion,
142  int *localenerpos
143  );
144 
145 
146 int bound_x1up_radcylbeam(
147  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
148  int *inboundloop,
149  int *outboundloop,
150  int *innormalloop,
151  int *outnormalloop,
152  int (*inoutlohi)[NUMUPDOWN][NDIM],
153  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
154  int *dosetbc,
155  int enerregion,
156  int *localenerpos
157  );
158 
159 
160 int bound_radbeam2dksvertbeaminflow(int dir,
161  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
162  int *inboundloop,
163  int *outboundloop,
164  int *innormalloop,
165  int *outnormalloop,
166  int (*inoutlohi)[NUMUPDOWN][NDIM],
167  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
168  int *dosetbc,
169  int enerregion,
170  int *localenerpos
171  );
172 
173 int bound_radcylbeamcart(int dir,
174  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
175  int *inboundloop,
176  int *outboundloop,
177  int *innormalloop,
178  int *outnormalloop,
179  int (*inoutlohi)[NUMUPDOWN][NDIM],
180  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
181  int *dosetbc,
182  int enerregion,
183  int *localenerpos
184  );
185 
186 int bound_staticset(int dir,
187  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
188  int *inboundloop,
189  int *outboundloop,
190  int *innormalloop,
191  int *outnormalloop,
192  int (*inoutlohi)[NUMUPDOWN][NDIM],
193  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
194  int *dosetbc,
195  int enerregion,
196  int *localenerpos
197  );
198 
199 
200 
201 
202 int bound_waldmono(int dir,
203  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
204  int *inboundloop,
205  int *outboundloop,
206  int *innormalloop,
207  int *outnormalloop,
208  int (*inoutlohi)[NUMUPDOWN][NDIM],
209  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
210  int *dosetbc,
211  int enerregion,
212  int *localenerpos
213  );
214 
215 
216 
217 
218 
219 
220 /* bound array containing entire set of primitive variables */
221 
222 
223 
224 static int bound_prim_user_general(int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int ispstag, int* dirprim, FTYPE (*prim)[NSTORE2][NSTORE3][NPR]);
225 
226 
227 
228 int bound_prim_user_dir(int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
229 {
230  int dirprim[NPR];
231  int pl,pliter;
232 
233 
234  // specify location of primitives
235  PALLLOOP(pl) dirprim[pl]=CENT;
236  // dualfprintf(fail_file,"start bound_prim\n"); // CHANGINGMARK
237  bound_prim_user_general(boundstage, finalstep, boundtime, whichdir, boundvartype, BOUNDPRIMLOC, dirprim, prim);
238  // dualfprintf(fail_file,"end bound_prim\n"); // CHANGINGMARK
239 
240  if(WHICHPROBLEM==RADDONUT && OUTERDEATH==1){
241  if(whichdir==1 && boundvartype!=BOUNDPSTAGTYPE && boundvartype!=BOUNDPSTAGSIMPLETYPE){// assumes always calls whichdir==1, which is true if N1>1
242  extern void debugfixupaltdeath_bc(FTYPE (*prim)[NSTORE2][NSTORE3][NPR]);
243  debugfixupaltdeath_bc(prim);
244  }
245  }
246 
247  return(0);
248 }
249 
250 
252 int bound_pstag_user_dir(int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
253 {
254 
255  int dirprim[NPR];
256  int pl,pliter;
257 
258 
259 
260 
261  if(FLUXB!=FLUXCTSTAG) return(0); // nothing to do
262 
263 
264  // specify location of primitives
265  PALLLOOP(pl) dirprim[pl]=CENT;
266  dirprim[B1]=FACE1;
267  dirprim[B2]=FACE2;
268  dirprim[B3]=FACE3;
269 
270 #if(0)
271  // Assume JCM setup bound_prim_user_general() correctly, so turned this off
272  // GODMARK: assume non-field velocity not set and have to have something reasonable
273  // use global values of non-field parts at present time
274  // note that bound_pstag() setup loops to be over only B1..B3, but user may violate this and just stick in something so no failures even if not using data
275  FULLLOOP PLOOPNOB1(pl) MACP0A1(prim,i,j,k,pl)=MACP0A1(p,i,j,k,pl);
276  FULLLOOP PLOOPNOB2(pl) MACP0A1(prim,i,j,k,pl)=MACP0A1(p,i,j,k,pl);
277 #endif
278 
279 
280  // assume before calling this that bound_pstag() setup PLOOPINTERP so only doing B1,B2,B3 (even though user may not respect this in bound_prim_user_general() -- which is ok since non-field quantities in pstag aren't needed -- may be problem if user_general() assumes primitive is reasonable)
281  // dualfprintf(fail_file,"start bound_pstag\n"); // CHANGINGMARK
282  bound_prim_user_general(boundstage, finalstep, boundtime, whichdir, boundvartype, BOUNDPSTAGLOC, dirprim, prim);
283  // dualfprintf(fail_file,"end bound_pstag\n"); // CHANGINGMARK
284 
285 
286 
287  return (0);
288 }
289 
290 
291 
293 int bound_prim_user_general(int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int ispstag, int* dirprim, FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
294 {
295  int inboundloop[NDIM];
296  int outboundloop[NDIM];
297  int innormalloop[NDIM];
298  int outnormalloop[NDIM];
299  int inoutlohi[NUMUPDOWN][NUMUPDOWN][NDIM];
300  int riin,riout,rjin,rjout,rkin,rkout;
301  int dosetbc[COMPDIM*2];
302  int enerregion;
303  int *localenerpos;
304  int dir;
305  int donebc[COMPDIM*2];
306 
307 
308 
309  DIRLOOP(dir) donebc[dir]=0; // assume BC not done yet
310 
311 
313  //
314  // set bound loop
315  //
317  set_boundloop(boundvartype, inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi, &riin, &riout, &rjin, &rjout, &rkin, &rkout, dosetbc);
318 
319 
320  // for CZLOOP:
321  // avoid looping over region outside active+ghost grid
322  // good because somewhat general and avoid bad inversions, etc.
323  // enerregion=TRUEGLOBALENERREGION;
324  enerregion=ACTIVEREGION; // now replaces TRUEGLOBALENERREGION
325  localenerpos=enerposreg[enerregion];
326  // if(WITHINENERREGION(localenerpos,i,j,k)){
327  // note that localenerpos[X1DN] sets first evolved cell
328 
329 
330 
331 
332  // first do non-directional internal grid assignments overwritting any evolution (ok to do for each whichdir)
333  if(WHICHPROBLEM==RADDOT){
334  bound_raddot(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
335  }
336 
337 
338 
339  if(whichdir==1){
340 
341  if((dosetbc[X1DN] || dosetbc[X1UP]) && (donebc[X1DN]==0 && donebc[X1UP]==0)){
342  if( (BCtype[X1DN]==PERIODIC)&&(BCtype[X1UP]==PERIODIC) ){
343  bound_x1_periodic(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
344  donebc[X1DN]=donebc[X1UP]=1;
345  }
346  }
347 
348  dir=X1DN;
349  if(dosetbc[dir] && donebc[dir]==0){
350  if(BCtype[dir]==HORIZONOUTFLOW){
351  bound_x1dn_outflow(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
352  donebc[dir]=1;
353  }
354  else if(BCtype[dir]==HORIZONOUTFLOWSTATIC){
355  BCtype[dir]=HORIZONOUTFLOW;
356  bound_x1dn_outflow(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
357  BCtype[dir]=HORIZONOUTFLOWSTATIC;
358  bound_staticset(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
359  donebc[dir]=1;
360  }
361  else if((BCtype[dir]==OUTFLOW)||(BCtype[dir]==FIXEDOUTFLOW)||(BCtype[dir]==FREEOUTFLOW)){
362  bound_x1dn_outflow_simple(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
363  donebc[dir]=1;
364  }
365  else if((BCtype[dir]==R0SING)||(BCtype[dir]==SYMM)||(BCtype[dir]==ASYMM) ){
366  bound_x1dn_sym(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
367  donebc[dir]=1;
368  }
369  else if(BCtype[dir]==FIXEDUSEPANALYTIC){
370  bound_x1dn_analytic(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim);
371  donebc[dir]=1;
372  }
373  else if(BCtype[dir]==RADBEAMFLATINFLOW){
374  bound_x1dn_radbeamflatinflow(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
375  donebc[dir]=1;
376  }
377  else if(BCtype[dir]==RADSHADOWINFLOW){
378  bound_radshadowinflow(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
379  donebc[dir]=1;
380  }
381  else if(BCtype[dir]==RADATMBEAMINFLOW){
382  bound_radatmbeaminflow(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
383  donebc[dir]=1;
384  }
385  else if(BCtype[dir]==RADWALLINFLOW){
386  bound_radwallinflow(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
387  donebc[dir]=1;
388  }
389  else if(BCtype[dir]==CYLAXIS){
390  bound_x1dn_cylaxis(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
391  donebc[dir]=1;
392  }
393  else if(BCtype[dir]==RADCYLBEAMCARTBC){
394  bound_radcylbeamcart(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
395  donebc[dir]=1;
396  }
397  else{
398  dualfprintf(fail_file,"No x1dn boundary condition specified: %d\n",BCtype[dir]);
399  myexit(7598730);
400  }
401  }
402 
403 
404  dir=X1UP;
405  if(dosetbc[dir] && donebc[dir]==0){
406  if((BCtype[dir]==OUTFLOW)||(BCtype[dir]==FIXEDOUTFLOW)||(BCtype[dir]==FREEOUTFLOW)){
407  bound_x1up_outflow_simple(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
408  donebc[dir]=1;
409  }
410  else if(BCtype[dir]==HORIZONOUTFLOW){
411  bound_x1up_outflow(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
412  donebc[dir]=1;
413  }
414  else if(BCtype[dir]==HORIZONOUTFLOWSTATIC){
415  BCtype[dir]=HORIZONOUTFLOW;
416  bound_x1up_outflow(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
417  BCtype[dir]=HORIZONOUTFLOWSTATIC;
418  bound_staticset(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
419  donebc[dir]=1;
420  }
421  else if(BCtype[dir]==RADBEAM2DFLOWINFLOW){
422  bound_radbeam2dflowinflow(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
423  donebc[dir]=1;
424  }
425  else if(BCtype[dir]==RADATMBEAMINFLOW){
426  bound_radatmbeaminflow(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
427  donebc[dir]=1;
428  }
429  else if(BCtype[dir]==FIXEDUSEPANALYTIC){
430  bound_x1up_analytic(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim);
431  donebc[dir]=1;
432  }
433  else if(BCtype[dir]==RADBONDIINFLOW){
434  bound_radbondiinflow(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
435  donebc[dir]=1;
436  }
437  else if(BCtype[dir]==RADNTBC){
438  bound_radnt(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
439  donebc[dir]=1;
440  }
441  else if(BCtype[dir]==RADCYLBEAMBC){
442  bound_x1up_radcylbeam(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
443  donebc[dir]=1;
444  }
445  else if(BCtype[dir]==RADCYLBEAMCARTBC){
446  bound_radcylbeamcart(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
447  donebc[dir]=1;
448  }
449  else if(BCtype[dir]==WALDMONOBC){
450  bound_waldmono(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
451  donebc[dir]=1;
452  }
453  else if(BCtype[dir]==RADCYLJETBC){
454  bound_x1up_radcyljet(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
455  donebc[dir]=1;
456  }
457  else{
458  dualfprintf(fail_file,"No x1up boundary condition specified: %d\n",BCtype[dir]);
459  myexit(7598731);
460  }
461  }
462 
463 
464  }
465  else if(whichdir==2){
466 
467 
468  if((dosetbc[X2DN] || dosetbc[X2UP]) && (donebc[X2DN]==0 && donebc[X2UP]==0)){
469  if( (BCtype[X2DN]==PERIODIC)&&(BCtype[X2UP]==PERIODIC) ){
470  bound_x2_periodic(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
471  donebc[X2DN]=donebc[X2UP]=1;
472  }
473  }
474 
475 
476  dir=X2DN;
477  if(dosetbc[dir] && donebc[dir]==0){
478  if((BCtype[dir]==OUTFLOW)||(BCtype[dir]==FIXEDOUTFLOW)||(BCtype[dir]==FREEOUTFLOW)){
479  bound_x2dn_outflow_simple(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
480  donebc[dir]=1;
481  }
482  else if(BCtype[dir]==POLARAXIS && special3dspc){
483  int whichcall=1;
484  bound_x2dn_polaraxis_full3d(whichcall,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
485  donebc[dir]=1;
486  }
487  else if((BCtype[dir]==POLARAXIS)||(BCtype[dir]==SYMM)||(BCtype[dir]==ASYMM) ){
488  bound_x2dn_polaraxis(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
489  donebc[dir]=1;
490  }
491  else if(BCtype[dir]==FIXEDUSEPANALYTIC){
492  bound_x2dn_analytic(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim);
493  donebc[dir]=1;
494  }
495  else if(BCtype[dir]==RADCYLBEAMCARTBC){
496  bound_radcylbeamcart(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
497  donebc[dir]=1;
498  }
499  else if(BCtype[dir]==RADCYLJETBC){
500  bound_x2dn_radcyljet(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
501  donebc[dir]=1;
502  }
503  else{
504  dualfprintf(fail_file,"No x2dn boundary condition specified: %d\n",BCtype[dir]);
505  myexit(7598732);
506  }
507  }
508 
509 
510  dir=X2UP;
511  if(dosetbc[dir] && donebc[dir]==0){
512  if((BCtype[dir]==OUTFLOW)||(BCtype[dir]==FIXEDOUTFLOW)||(BCtype[dir]==FREEOUTFLOW)){
513  bound_x2up_outflow_simple(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
514  donebc[dir]=1;
515  }
516  else if(BCtype[dir]==POLARAXIS && special3dspc){
517  int whichcall=1;
518  bound_x2up_polaraxis_full3d(whichcall,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
519  donebc[dir]=1;
520  }
521  else if((BCtype[dir]==POLARAXIS)||(BCtype[dir]==SYMM)||(BCtype[dir]==ASYMM) ){
522  bound_x2up_polaraxis(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
523  donebc[dir]=1;
524  }
525  else if(BCtype[dir]==FIXEDUSEPANALYTIC){
526  bound_x2up_analytic(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim);
527  donebc[dir]=1;
528  }
529  else if(BCtype[dir]==RADSHADOWINFLOWX2UP){
530  bound_radshadowinflow(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
531  donebc[dir]=1;
532  }
533  else if(BCtype[dir]==RADWALLINFLOW){
534  bound_radwallinflow(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
535  donebc[dir]=1;
536  }
537  else if(BCtype[dir]==RADNTBC){
538  // do ASYMM first
539  BCtype[dir]=ASYMM;
540  bound_x2up_polaraxis(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
541  // now do anything else that needs to be done
542  BCtype[dir]=RADNTBC;
543  bound_radnt(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
544  donebc[dir]=1;
545  }
546  else if(BCtype[dir]==RADBEAM2DKSVERTBEAMINFLOW){
547  bound_radbeam2dksvertbeaminflow(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
548  donebc[dir]=1;
549  }
550  else if(BCtype[dir]==RADCYLBEAMCARTBC){
551  bound_radcylbeamcart(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
552  donebc[dir]=1;
553  }
554  else if(BCtype[dir]==WALDMONOBC){
555  // do ASYMM first
556  BCtype[dir]=ASYMM;
557  bound_x2up_polaraxis(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
558  // now do anything else that needs to be done
559  BCtype[dir]=WALDMONOBC;
560  bound_waldmono(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
561  donebc[dir]=1;
562  }
563  else{
564  dualfprintf(fail_file,"No x2dn boundary condition specified: %d\n",BCtype[dir]);
565  myexit(7598733);
566  }
567  }
568 
569 
570 
571 
572  }
573  else if(whichdir==3){
574 
575 
576 
577  if((dosetbc[X3DN] || dosetbc[X3UP]) && (donebc[X3DN]==0 && donebc[X3UP]==0)){
578  if( (BCtype[X3DN]==PERIODIC)&&(BCtype[X3UP]==PERIODIC) ){
579  bound_x3_periodic(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
580  donebc[X3DN]=donebc[X3UP]=1;
581  }
582  }
583 
584 
585  dir=X3DN;
586  if(dosetbc[dir] && donebc[dir]==0){
587  if((BCtype[dir]==OUTFLOW)||(BCtype[dir]==FIXEDOUTFLOW)||(BCtype[dir]==FREEOUTFLOW)){
588  bound_x3dn_outflow_simple(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
589  donebc[dir]=1;
590  }
591  else if(BCtype[dir]==FIXEDUSEPANALYTIC){
592  bound_x3dn_analytic(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim);
593  donebc[dir]=1;
594  }
595  else if(BCtype[dir]==RADBEAM2DBEAMINFLOW){
596  bound_radbeam2dbeaminflow(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
597  donebc[dir]=1;
598  }
599  else if(BCtype[dir]==RADCYLBEAMCARTBC){
600  bound_radcylbeamcart(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
601  donebc[dir]=1;
602  }
603  else{
604  dualfprintf(fail_file,"No x3dn boundary condition specified: %d\n",BCtype[dir]);
605  myexit(34672546);
606  }
607  }
608 
609 
610  dir=X3UP;
611  if(dosetbc[dir] && donebc[dir]==0){
612  if((BCtype[dir]==OUTFLOW)||(BCtype[dir]==FIXEDOUTFLOW)||(BCtype[dir]==FREEOUTFLOW)){
613  bound_x3up_outflow_simple(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
614  donebc[dir]=1;
615  }
616  else if(BCtype[dir]==FIXEDUSEPANALYTIC){
617  bound_x3up_analytic(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim);
618  donebc[dir]=1;
619  }
620  else if(BCtype[dir]==RADCYLBEAMCARTBC){
621  bound_radcylbeamcart(dir,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
622  donebc[dir]=1;
623  }
624  else{
625  dualfprintf(fail_file,"No x3up boundary condition specified: %d\n",BCtype[dir]);
626  myexit(34672547);
627  }
628  }
629 
630 
631 
632 
633 
634  }
635  else{
636  dualfprintf(fail_file,"No such whichdir=%d\n",whichdir);
637  myexit(2436262);
638  }
639 
640 
641 
642 
643 
644 
645 
646  return(0);
647 }
648 
649 
650 
651 
652 
654 int apply_bc_line(int nprlocalstart, int nprlocalend, int*nprlocallist, int doinverse, int iterdir, int recontype, int bs, int be, FTYPE (*yin)[2][NBIGM], FTYPE (*yout)[2][NBIGM], FTYPE (*youtpolycoef)[MAXSPACEORDER][NBIGM])
655 {
656  int flip_y(int nprlocalstart, int nprlocalend, int*nprlocallist, int iterdir, int recontype, int bs, int be, FTYPE (*y)[2][NBIGM]);
657 
658  if(doinverse==0){
659  flip_y(nprlocalstart,nprlocalend,nprlocallist,iterdir, recontype, bs, be, yin);
660  }
661  else{
662  flip_y(nprlocalstart,nprlocalend,nprlocallist,iterdir, recontype, bs, be, yin);
663  flip_y(nprlocalstart,nprlocalend,nprlocallist,iterdir, recontype, bs, be, yout);
664  }
665 
666  return(0);
667 
668 }
669 
670 
671 #include "reconstructeno.h"
672 
673 int flip_y(int nprlocalstart, int nprlocalend, int*nprlocallist, int iterdir, int recontype, int bs, int be, FTYPE (*y)[2][NBIGM])
674 {
675  int pllocal,pl,myi;
676 
677 
678 #if( WENO_DIR_FLIP_CONS_SIGN_DN ) //flip the sign of the consrved quantities at the cylindrical axis so that they do not have a kink due to multiplication by gdet = |R|
679  if( iterdir == WENO_DIR_FLIP_CONS_SIGN_DN && (recontype == CVT_C2A || recontype == CVT_A2C) && mycpupos[iterdir] == 0 ) {
680  NUMPRIMLOOP(pllocal,pl)
681  for( myi = bs; myi < 0; myi++ ) {
682  y[pl][0][myi] = - y[pl][0][myi];
683  }
684  }
685 #endif
686 
687 #if( WENO_DIR_FLIP_CONS_SIGN_UP ) //flip the sign of the consrved quantities at the cylindrical axis so that they do not have a kink due to multiplication by gdet = |R|
688  if( iterdir == WENO_DIR_FLIP_CONS_SIGN_UP && (recontype == CVT_C2A || recontype == CVT_A2C) && mycpupos[iterdir] == numbercpu[iterdir] - 1 ) {
689  NUMPRIMLOOP(pllocal,pl)
690  for( myi = N1*(iterdir==1) + N2*(iterdir==2) + N3*(iterdir==3); myi <= be; myi++ ) {
691  y[pl][0][myi] = - y[pl][0][myi];
692  }
693  }
694 #endif
695 
696 
697  return(0);
698 
699 }
700 
701 void remapdq( int dir, int idel, int jdel, int kdel, int i, int j, int k, FTYPE (*p2interp)[NSTORE2][NSTORE3][NPR2INTERP],
702  FTYPE (*dq)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pleft)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pright)[NSTORE2][NSTORE3][NPR2INTERP],
703  FTYPE *p2interp_l, FTYPE *p2interp_r )
704 {
705 }
706 
707 void remapplpr( int dir, int idel, int jdel, int kdel, int i, int j, int k,
708  FTYPE (*p2interp)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*dq)[NSTORE2][NSTORE3][NPR2INTERP],
709  FTYPE (*pleft)[NSTORE2][NSTORE3][NPR2INTERP], FTYPE (*pright)[NSTORE2][NSTORE3][NPR2INTERP],
710  FTYPE *p2interp_l, FTYPE *p2interp_r )
711 {
712 }
713 
714 
717 int bound_prim_user_after_mpi_dir(int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR])
718 {
719  int dirprim[NPR];
720  int pliter,pl;
721 
722  int inboundloop[NDIM];
723  int outboundloop[NDIM];
724  int innormalloop[NDIM];
725  int outnormalloop[NDIM];
726  int inoutlohi[NUMUPDOWN][NUMUPDOWN][NDIM];
727  int riin,riout,rjin,rjout,rkin,rkout;
728  int dosetbc[COMPDIM*2];
729  int enerregion;
730  int *localenerpos;
731  int dir;
732 
733 
734 
735 
737  //
738  // set bound loop
739  //
741  set_boundloop(boundvartype, inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi, &riin, &riout, &rjin, &rjout, &rkin, &rkout, dosetbc);
742 
743 
744  // for CZLOOP:
745  // avoid looping over region outside active+ghost grid
746  // good because somewhat general and avoid bad inversions, etc.
747  // enerregion=TRUEGLOBALENERREGION;
748  enerregion=ACTIVEREGION; // now replaces TRUEGLOBALENERREGION
749  localenerpos=enerposreg[enerregion];
750  // if(WITHINENERREGION(localenerpos,i,j,k)){
751  // note that localenerpos[X1DN] sets first evolved cell
752 
753 
754  if(ispstag==BOUNDPRIMLOC){
755  // specify location of primitives
756  PALLLOOP(pl) dirprim[pl]=CENT;
757  }
758 
759  if(ispstag==BOUNDPSTAGLOC){
760  // specify location of primitives
761  PALLLOOP(pl) dirprim[pl]=CENT;
762  dirprim[B1]=FACE1;
763  dirprim[B2]=FACE2;
764  dirprim[B3]=FACE3;
765  }
766 
767 
768  // use post-MPI values to operate poledeath() and/or polesmooth() since otherwise boundary values poledeath uses are not set yet
769  if(whichdir==2){
770 
771  dir=X2DN;
772  if(dosetbc[dir]){
773  if(BCtype[dir]==POLARAXIS && special3dspc){
774  bound_x2dn_polaraxis_full3d(2,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
775  }
776  }
777 
778 
779  dir=X2UP;
780  if(dosetbc[dir]){
781  if(BCtype[dir]==POLARAXIS && special3dspc){
782  bound_x2up_polaraxis_full3d(2,boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
783  }
784  }
785 
786  }
787 
788 
789  // can only really check boundaries after MPI is done
790  // e.g. periodicx3 and ncpux3==2 and then MPi required to set k<0 and k>=ncpux2*N3 cells
791  if(whichdir==1 && N2==1 && N3==1 || N3==1 && whichdir==2 || N3>1 && whichdir==3){ // not completely general conditional
792  bound_checks1(boundstage,finalstep,boundtime,whichdir,boundvartype,dirprim,ispstag,prim,inboundloop,outboundloop,innormalloop,outnormalloop,inoutlohi,riin,riout,rjin,rjout,rkin,rkout,dosetbc,enerregion,localenerpos);
793  }
794 
795 
796  return(0);
797 }
798 
799 
800 
801 
802 
803 
804 
805 
806 
807 
808 
809 
810 
811 
812 
818 
819 
820 
822 int bound_x1dn_radbeamflatinflow(
823  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
824  int *inboundloop,
825  int *outboundloop,
826  int *innormalloop,
827  int *outnormalloop,
828  int (*inoutlohi)[NUMUPDOWN][NDIM],
829  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
830  int *dosetbc,
831  int enerregion,
832  int *localenerpos
833  )
834 
835 {
836 
837 
838 #pragma omp parallel // assume don't require EOS
839  {
840 
841  int i,j,k,pl,pliter;
842  FTYPE vcon[NDIM],X[NDIM],V[NDIM];
843 #if(WHICHVEL==VEL3)
844  int failreturn;
845 #endif
846  int ri, rj, rk; // reference i,j,k
847  FTYPE prescale[NPR];
848  int jj,kk;
849  struct of_geom geomdontuse[NPR];
850  struct of_geom *ptrgeom[NPR];
851  struct of_geom rgeomdontuse[NPR];
852  struct of_geom *ptrrgeom[NPR];
853 
854  // assign memory
855  PALLLOOP(pl){
856  ptrgeom[pl]=&(geomdontuse[pl]);
857  ptrrgeom[pl]=&(rgeomdontuse[pl]);
858  }
859 
860 
861 
862  if(BCtype[X1DN]==RADBEAMFLATINFLOW && totalsize[1]>1 && mycpupos[1] == 0 ) {
863 
864 
865 
866  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
868 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
869  OPENMPBCLOOPBLOCK{
870  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
871 
872 
873 
874  ri=riin;
875  rj=j;
876  rk=k;
877 
878 
879  // ptrrgeom : i.e. ref geom
880  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
881 
882 
883 
884  LOOPBOUND1IN{
885  FTYPE *pr = &MACP0A1(prim,i,j,k,0);
886 
887  //initially copying everything
888  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
889 
890 
891  if(ispstag==0){ // only do something special with non-field primitives
892 
893  // local geom
894  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
895 
896  //coordinates of the ghost cell
897  bl_coord_ijk_2(i,j,k,CENT,X, V);
898 
899  // set radiation quantities as R^t_\nu in orthonormal fluid frame using whichvel velocity and whichcoord coordinates
900 
901  extern FTYPE RADBEAMFLAT_FRATIO,RADBEAMFLAT_ERAD, RADBEAMFLAT_RHO, RADBEAMFLAT_UU;
902 
903  // FTYPE ERADAMB=RADBEAMFLAT_ERAD;
904  // FTYPE ERADINJ=1000.0*ERADAMB;
905 
906  FTYPE ERADAMB=RADBEAMFLAT_ERAD;
907  // FTYPE ERADINJ=1000.0*ERADAMB; // old harm setup
908  FTYPE ERADINJ=100.0*ERADAMB; // koral's current setup
909 
910 
911 
912  if(0){
913  // correct version in general
914 
915  pr[RHO] = RADBEAMFLAT_RHO ;
916  pr[UU] = RADBEAMFLAT_UU;
917  SLOOPA(jj) pr[U1+jj-1] = 0.0;
918 
919  //E, F^i in orthonormal fluid frame
920  FTYPE pradffortho[NPR];
921  FTYPE Fx=0,Fy=0,Fz=0;
922  //primitives in whichvel,whichcoord
923  if(V[2]>.4 && V[2]<.6){//beam to be imposed
924  Fx=RADBEAMFLAT_FRATIO*ERADINJ;
925  Fy=Fz=0.0;
926 
927  pradffortho[PRAD0] = ERADINJ;
928  pradffortho[PRAD1] = Fx;
929  pradffortho[PRAD2] = Fy;
930  pradffortho[PRAD3] = Fz;
931  }
932  else{ //no beam
933  Fx=Fy=Fz=0.0;
934  pradffortho[PRAD0] = ERADAMB;
935  pradffortho[PRAD1] = Fx;
936  pradffortho[PRAD2] = Fy;
937  pradffortho[PRAD3] = Fz;
938  }
939 
940  int whichvel=VEL4;
941  int whichcoordfluid=MCOORD;
942  int whichcoordrad=whichcoordfluid;
943  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
944  }
945  else if(0){
946  // set vradx
947 
948  pr[RHO] = RADBEAMFLAT_RHO ;
949  pr[UU] = RADBEAMFLAT_UU;
950  SLOOPA(jj) pr[U1+jj-1] = 0.0;
951 
952  // assume RADBEAMFLAT_FRATIO is just vradx
953  FTYPE uradx=1.0/sqrt(1.0 - RADBEAMFLAT_FRATIO*RADBEAMFLAT_FRATIO); // radiation 4-velocity
954 
955  //primitives in whichvel,whichcoord
956  if(V[2]>.4 && V[2]<.6){//beam to be imposed
957  pr[URAD0] = ERADINJ;
958  pr[URAD1] = uradx;
959  pr[URAD2] = 0.;
960  pr[URAD3] = 0.;
961  }
962  else{ //no beam
963  pr[URAD0] = ERADAMB;
964  pr[URAD1] = 0.;
965  pr[URAD2] = 0.;
966  pr[URAD3] = 0.;
967  }
968 
969  // get all primitives in WHICHVEL/PRIMECOORDS value
970  int whichvel=VEL4;
971  int whichcoord=MCOORD;
972  primefluid_EVrad_to_primeall(&whichvel, &whichcoord, ptrgeom[RHO],MAC(prim,i,j,k),MAC(prim,i,j,k)); // assumes ptrgeom[RHO] is same location as all other primitives (as is currently true).
973  }
974  else if(1){
975  // koral mixed way (must use ff ortho choice in init.koral.c (i.e. first if(1))
976 
977  pr[RHO] = RADBEAMFLAT_RHO ;
978  pr[UU] = RADBEAMFLAT_UU;
979  SLOOPA(jj) pr[U1+jj-1] = 0.0;
980 
981 
982  //primitives in whichvel,whichcoord
983  if(V[2]>.4 && V[2]<.6){//beam to be imposed
984  FTYPE dxdxp[NDIM][NDIM];
985  dxdxprim_ijk(i, j, k, CENT, dxdxp);
986 
987  // like koral
988  // FTYPE uradx=ERADINJ*RADBEAMFLAT_FRATIO;
989 
990  // fix:
991  FTYPE uradx=1.0/sqrt(1.0 - RADBEAMFLAT_FRATIO*RADBEAMFLAT_FRATIO); // radiation 4-velocity
992 
993 
994  pr[URAD0] = ERADINJ;
995  pr[URAD1] = uradx/dxdxp[1][1];
996  pr[URAD2] = 0.;
997  pr[URAD3] = 0.;
998  }
999  else{ //no beam and ERADAMB is assumed as radiation frame -- as in koral.
1000  pr[URAD0] = ERADAMB;
1001  pr[URAD1] = 0.;
1002  pr[URAD2] = 0.;
1003  pr[URAD3] = 0.;
1004  }
1005 
1006  }
1007 
1008 
1009  }// end if not staggered field
1010 
1011 
1012  }// end loop over inner i's
1013  }
1014  }
1015 
1016 
1017  }// end parallel region
1018 
1019  return(0);
1020 }
1021 
1022 
1023 
1024 
1026 int bound_radshadowinflow(int dir,
1027  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
1028  int *inboundloop,
1029  int *outboundloop,
1030  int *innormalloop,
1031  int *outnormalloop,
1032  int (*inoutlohi)[NUMUPDOWN][NDIM],
1033  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
1034  int *dosetbc,
1035  int enerregion,
1036  int *localenerpos
1037  )
1038 
1039 {
1040 
1041  if(dir==X2DN){
1042  dualfprintf(fail_file,"Shouldn't be in bound_radshadowinflow() with dir=%d -- should be using ASYMM for RADDBLSHADOW\n");
1043  myexit(87243534);
1044  }
1045 
1046 
1047  FTYPE RHOAMB=1.e-4; // matches init.koral.c
1048  FTYPE TAMB=1.e7/TEMPBAR; // matches init.koral.c
1049  FTYPE BLOBW=0.22;
1050  FTYPE RHOBLOB=1.e3;
1051  //
1052  extern FTYPE RADSHADOW_NLEFT,RADSHADOW_ANGLE;
1053  extern FTYPE RADSHADOW_TLEFTOTAMB;
1054  extern FTYPE RADSHADOW_BEAMY;
1055 
1056  extern FTYPE RADDBLSHADOW_NLEFT,RADDBLSHADOW_ANGLE;
1057  extern FTYPE RADDBLSHADOW_TLEFTOTAMB;
1058  extern FTYPE RADDBLSHADOW_BEAMY;
1059 
1060  FTYPE NLEFT,angle,TLEFT,BEAMY;
1061 
1062  if(WHICHPROBLEM==RADSHADOW){
1063  NLEFT=RADSHADOW_NLEFT;
1064  angle=RADSHADOW_ANGLE;
1065  TLEFT=TAMB*RADSHADOW_TLEFTOTAMB;
1066  BEAMY=RADSHADOW_BEAMY;
1067  }
1068  else if(WHICHPROBLEM==RADDBLSHADOW){
1069  NLEFT=RADDBLSHADOW_NLEFT;
1070  angle=RADDBLSHADOW_ANGLE;
1071  TLEFT=TAMB*RADDBLSHADOW_TLEFTOTAMB;
1072  BEAMY=RADDBLSHADOW_BEAMY;
1073  }
1074 
1075 
1076 #pragma omp parallel // assume don't require EOS
1077  {
1078  // make rho,u consistent with on-domain values
1079  FTYPE rho;
1080  FTYPE uint;
1081  FTYPE Trad;
1082 
1083  int i,j,k,pl,pliter;
1084  FTYPE vcon[NDIM],X[NDIM],V[NDIM];
1085 #if(WHICHVEL==VEL3)
1086  int failreturn;
1087 #endif
1088  int ri, rj, rk; // reference i,j,k
1089  FTYPE prescale[NPR];
1090  int jj,kk;
1091  struct of_geom geomdontuse[NPR];
1092  struct of_geom *ptrgeom[NPR];
1093  struct of_geom rgeomdontuse[NPR];
1094  struct of_geom *ptrrgeom[NPR];
1095 
1096  // assign memory
1097  PALLLOOP(pl){
1098  ptrgeom[pl]=&(geomdontuse[pl]);
1099  ptrrgeom[pl]=&(rgeomdontuse[pl]);
1100  }
1101 
1102 
1103 
1105  // dir==X1DN
1107  if(dir==X1DN && BCtype[X1DN]==RADSHADOWINFLOW && (totalsize[1]>1) && (mycpupos[1] == 0) ){
1108 
1109 
1110  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
1112 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1113  OPENMPBCLOOPBLOCK{
1114  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
1115 
1116 
1117  ri=riin;
1118  rj=j;
1119  rk=k;
1120 
1121 
1122  // ptrrgeom : i.e. ref geom
1123  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
1124 
1125 
1126  LOOPBOUND1IN{
1127  FTYPE *pr=&MACP0A1(prim,i,j,k,0);
1128 
1129  //initially copying everything
1130  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
1131 
1132  if(ispstag==0){
1133  // local geom
1134  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
1135 
1136  //coordinates of the ghost cell
1137  bl_coord_ijk_2(i,j,k,CENT,X, V);
1138  FTYPE xx,yy,zz,rsq;
1139  xx=V[1];
1140  yy=V[2];
1141  zz=V[3];
1142 
1143  // default for beam or not
1144  rsq=xx*xx+yy*yy+zz*zz;
1145  rho=(RHOBLOB-RHOAMB)*exp(-sqrt(rsq)/(BLOBW*BLOBW))+RHOAMB;
1146  pr[RHO] = rho;
1147  Trad=TAMB*RHOAMB/rho;
1148  uint=calc_PEQ_ufromTrho(Trad,rho);
1149  pr[UU] = uint;
1150 
1151 
1152  if(yy>BEAMY && WHICHPROBLEM==RADDBLSHADOW || WHICHPROBLEM==RADSHADOW ){
1153 
1154  FTYPE ERAD;
1155  ERAD=calc_LTE_EfromT(TLEFT);
1156 
1157  FTYPE ux=0.0; // orthonormal 4-velocity. Matches init.koral.c
1158  pr[U1] = ux/sqrt(ptrgeom[U1]->gcov[GIND(1,1)]); // assumed no spatial mixing
1159 
1160  //E, F^i
1161  if(1){
1162  // correct way of using ERAD(TLEFT) and NLEFT
1163  FTYPE Fx=NLEFT*ERAD/sqrt(1+angle*angle);
1164  FTYPE Fy=-NLEFT*ERAD*angle/sqrt(1+angle*angle);
1165  FTYPE Fz=0.0;
1166 
1167  //E, F^i in orthonormal fluid frame
1168  FTYPE pradffortho[NPR];
1169  pradffortho[PRAD0] = ERAD;
1170  pradffortho[PRAD1] = Fx;
1171  pradffortho[PRAD2] = Fy;
1172  pradffortho[PRAD3] = Fz;
1173 
1174 
1175  int whichvel=VEL4; // in which vel U1-U3 set
1176  int whichcoordfluid=MCOORD; // in which coordinates U1-U3 set
1177  int whichcoordrad=whichcoordfluid; // in which coordinates E,F are orthonormal
1178  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
1179  }
1180  else{
1181  // old harm way
1182  FTYPE gammax=1.0/sqrt(1.0-NLEFT*NLEFT);
1183  FTYPE uradx=NLEFT*gammax/sqrt(1.0+angle*angle);
1184  FTYPE urady=-NLEFT*gammax*angle/sqrt(1.0+angle*angle);
1185 
1186  pr[URAD0] = ERAD;
1187  pr[URAD1] = uradx;
1188  pr[URAD2] = urady;
1189  pr[URAD3] = 0.;
1190 
1191  // get all primitives in WHICHVEL/PRIMECOORDS value
1192  int whichvel;
1193  whichvel=VEL4;
1194  int whichcoord;
1195  whichcoord=MCOORD;
1196  primefluid_EVrad_to_primeall(&whichvel, &whichcoord, ptrgeom[RHO],MAC(prim,i,j,k),MAC(prim,i,j,k)); // assumes ptrgeom[RHO] is same location as all other primitives (as is currently true).
1197  }
1198 
1199  } // over spatially relevant region
1200 
1201  }// end if not staggered fields
1202 
1203  }// end loop over inner i's
1204 
1205  }// end over loop
1206  }// end if correct boundary condition and core
1207 
1208 
1210  // dir==X2DN handled via AYMM BC
1212 
1213 
1215  // dir==X2UP
1217  if(dir==X2UP && BCtype[X2UP]==RADSHADOWINFLOWX2UP && (totalsize[2]>1) && (mycpupos[2] == ncpux2-1) ){
1218 
1219 
1220  OPENMPBCLOOPVARSDEFINELOOPX2DIR; OPENMPBCLOOPSETUPLOOPX2DIR;
1222 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1223  OPENMPBCLOOPBLOCK{
1224  OPENMPBCLOOPBLOCK2IJKLOOPX2DIR(i,k);
1225 
1226 
1227  ri=i;
1228  rj=rjout;
1229  rk=k;
1230 
1231 
1232  // ptrrgeom : i.e. ref geom
1233  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
1234 
1235 
1236  LOOPBOUND2OUT{
1237  FTYPE *pr=&MACP0A1(prim,i,j,k,0);
1238 
1239  //initially copying everything
1240  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
1241 
1242  if(ispstag==0){
1243  // local geom
1244  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
1245 
1246  //coordinates of the ghost cell
1247  bl_coord_ijk_2(i,j,k,CENT,X, V);
1248  FTYPE xx,yy,zz,rsq;
1249  xx=V[1];
1250  yy=V[2];
1251  zz=V[3];
1252 
1253 
1254 
1255  // default for beam or not
1256  rsq=xx*xx+yy*yy+zz*zz;
1257  rho=(RHOBLOB-RHOAMB)*exp(-sqrt(rsq)/(BLOBW*BLOBW))+RHOAMB;
1258  pr[RHO] = rho;
1259  Trad=TAMB*RHOAMB/rho;
1260  uint= calc_PEQ_ufromTrho(Trad,rho);
1261  pr[UU] = uint;
1262 
1263 
1264  if(WHICHPROBLEM==RADDBLSHADOW || WHICHPROBLEM==RADSHADOW){
1265 
1266  FTYPE ERAD;
1267  ERAD=calc_LTE_EfromT(TLEFT);
1268 
1269  FTYPE ux=0.0; // orthonormal 4-velocity. Matches init.koral.c
1270  pr[U1] = ux/sqrt(ptrgeom[U1]->gcov[GIND(1,1)]); // assumed no spatial mixing
1271 
1272  //E, F^i
1273  if(1){
1274  // correct way of using ERAD(TLEFT) and NLEFT
1275  FTYPE Fx=NLEFT*ERAD/sqrt(1+angle*angle);
1276  FTYPE Fy=-NLEFT*ERAD*angle/sqrt(1+angle*angle);
1277  FTYPE Fz=0.0;
1278 
1279  //E, F^i in orthonormal fluid frame
1280  FTYPE pradffortho[NPR];
1281  pradffortho[PRAD0] = ERAD;
1282  pradffortho[PRAD1] = Fx;
1283  pradffortho[PRAD2] = Fy;
1284  pradffortho[PRAD3] = Fz;
1285 
1286 
1287  int whichvel=VEL4; // in which vel U1-U3 set
1288  int whichcoordfluid=MCOORD; // in which coordinates U1-U3 set
1289  int whichcoordrad=whichcoordfluid; // in which coordinates E,F are orthonormal
1290  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
1291  }
1292  else{
1293  // old harm way
1294  FTYPE gammax=1.0/sqrt(1.0-NLEFT*NLEFT);
1295  FTYPE uradx=NLEFT*gammax/sqrt(1.0+angle*angle);
1296  FTYPE urady=-NLEFT*gammax*angle/sqrt(1.0+angle*angle);
1297 
1298  pr[URAD0] = ERAD;
1299  pr[URAD1] = uradx;
1300  pr[URAD2] = urady;
1301  pr[URAD3] = 0.;
1302 
1303  // get all primitives in WHICHVEL/PRIMECOORDS value
1304  int whichvel;
1305  whichvel=VEL4;
1306  int whichcoord;
1307  whichcoord=MCOORD;
1308  primefluid_EVrad_to_primeall(&whichvel, &whichcoord, ptrgeom[RHO],MAC(prim,i,j,k),MAC(prim,i,j,k)); // assumes ptrgeom[RHO] is same location as all other primitives (as is currently true).
1309  }
1310 
1311 
1312 
1313  }// if spatially relevant region
1314 
1315  }// end if not staggered fields
1316 
1317  }// end loop over inner j's
1318  }// end over loop
1319  }// end if correct boundary condition and core
1320 
1321 
1322  }// end parallel region
1323 
1324 
1325 
1326 
1327 
1328 
1329 
1330  return(0);
1331 }
1332 
1333 
1334 
1335 
1336 
1338 int bound_radbeam2dbeaminflow(int dir,
1339  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
1340  int *inboundloop,
1341  int *outboundloop,
1342  int *innormalloop,
1343  int *outnormalloop,
1344  int (*inoutlohi)[NUMUPDOWN][NDIM],
1345  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
1346  int *dosetbc,
1347  int enerregion,
1348  int *localenerpos
1349  )
1350 
1351 {
1352 
1353 
1354 #pragma omp parallel // assume don't require EOS
1355  {
1356 
1357  int i,j,k,pl,pliter;
1358  FTYPE X[NDIM],V[NDIM];
1359  int ri, rj, rk; // reference i,j,k
1360  int jj,kk;
1361  struct of_geom geomdontuse[NPR];
1362  struct of_geom *ptrgeom[NPR];
1363  struct of_geom rgeomdontuse[NPR];
1364  struct of_geom *ptrrgeom[NPR];
1365 
1366  // assign memory
1367  PALLLOOP(pl){
1368  ptrgeom[pl]=&(geomdontuse[pl]);
1369  ptrrgeom[pl]=&(rgeomdontuse[pl]);
1370  }
1371 
1372  extern int RADBEAM2D_BEAMNO;
1373  extern int RADBEAM2D_FLATBACKGROUND;
1374  extern FTYPE RADBEAM2D_RHOAMB;
1375  extern FTYPE RADBEAM2D_TAMB;
1376  extern int RADBEAM2D_BLOB;
1377  extern FTYPE RADBEAM2D_BLOBW;
1378  extern FTYPE RADBEAM2D_BLOBP;
1379  extern FTYPE RADBEAM2D_BLOBX;
1380  extern FTYPE RADBEAM2D_BLOBZ;
1381  extern FTYPE RADBEAM2D_PAR_D;
1382  extern FTYPE RADBEAM2D_PAR_E;
1383  extern int RADBEAM2D_IFBEAM;
1384  extern FTYPE RADBEAM2D_TLEFT;
1385  extern FTYPE RADBEAM2D_NLEFT;
1386  extern FTYPE RADBEAM2D_BEAML;
1387  extern FTYPE RADBEAM2D_BEAMR;
1388 
1389 
1390  if(dir==X3DN && BCtype[X3DN]==RADBEAM2DBEAMINFLOW && totalsize[3]>1 && mycpupos[3] == 0 ){
1391 
1392 
1393 
1394 
1395  OPENMPBCLOOPVARSDEFINELOOPX3DIR; OPENMPBCLOOPSETUPLOOPX3DIR;
1397 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1398  OPENMPBCLOOPBLOCK{
1399  OPENMPBCLOOPBLOCK2IJKLOOPX3DIR(i,j);
1400 
1401  ri=i;
1402  rj=j;
1403  rk=rkin;
1404 
1405 
1406  // ptrrgeom : i.e. ref geom
1407  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
1408 
1409 
1410  LOOPBOUND3IN{
1411  FTYPE *pr=&MACP0A1(prim,i,j,k,0);
1412 
1413  //initially copying everything
1414  PBOUNDLOOP(pliter,pl) pr[pl] = MACP0A1(prim,ri,rj,rk,pl);
1415 
1416 
1417  if(0){
1418  // NOTEMARK: only really makes sense near the hole if in KSCOORDS
1419  PBOUNDLOOP(pliter,pl) if(pl==U3) if(pr[U3]>0.0) pr[U3]=0.0; // limit so no arbitrary fluid inflow
1420  PBOUNDLOOP(pliter,pl) if(pl==URAD3) if(pr[URAD3]>0.0) pr[URAD3]=0.0; // limit so no arbitrary radiative inflow
1421  }
1422 
1423  // store this outflow result
1424  FTYPE pr0[NPR];
1425  PBOUNDLOOP(pliter,pl) pr0[pl]=pr[pl];
1426 
1427  // only overwrite copy if not inside i<0 since want to keep consistent with outflow BCs used for other \phi at those i
1428  if(1 && startpos[1]+i>=0 && ispstag==0){ // only do something special with non-field primitives
1429 
1430  // local geom
1431  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
1432 
1433  int whichcoord;
1434  whichcoord=MCOORD;
1435 
1436  // get metric grid geometry for these ICs
1437  int getprim=0;
1438  struct of_geom geomrealdontuse;
1439  struct of_geom *ptrgeomreal=&geomrealdontuse;
1440  gset(getprim,whichcoord,i,j,k,ptrgeomreal);
1441 
1442 
1443  //coordinates of the ghost cell
1444  bl_coord_ijk_2(i,j,k,CENT,X, V);
1445 
1446 
1447  FTYPE ERADAMB;
1448  FTYPE rho,uint,Vr;
1449  if(RADBEAM2D_FLATBACKGROUND){
1450  Vr=0.0;
1451  rho=RADBEAM2D_RHOAMB;
1452  uint=calc_PEQ_ufromTrho(RADBEAM2D_TAMB,rho);
1453  ERADAMB=calc_LTE_EfromT(RADBEAM2D_TAMB);
1454 
1455  // override so like outflow conditions to avoid shear at boundary
1456  if(0) ERADAMB=pr[URAD0]; // only makes sense with urad method
1457 
1458  }
1459  else{
1460  //zaczynam jednak od profilu analitycznego:
1461  FTYPE r=V[1];
1462  FTYPE mD=RADBEAM2D_PAR_D/(r*r*sqrt(2./r*(1.-2./r)));
1463  FTYPE mE=RADBEAM2D_PAR_E/(pow(r*r*sqrt(2./r),gamideal)*pow(1.-2./r,(gamideal+1.)/4.));
1464  Vr=sqrt(2./r)*(1.-2./r);
1465 
1466 
1467  FTYPE W=1./sqrt(1.-Vr*Vr*ptrgeomreal->gcov[GIND(1,1)]); // assumes RHO location is good for all these quantities
1468  rho=RADBEAM2D_PAR_D/(r*r*sqrt(2./r));
1469  FTYPE T=RADBEAM2D_TAMB;
1470  // FTYPE ERAD=calc_LTE_EfromT(T);
1471  uint=mE/W;
1472  ERADAMB=calc_LTE_Efromurho(uint,rho);
1473 
1474  // override so like outflow conditions to avoid shear at boundary
1475  if(0) ERADAMB=pr[URAD0];
1476  }
1477 
1478 
1479  // beam parameters
1480  FTYPE ERADINJ=calc_LTE_EfromT(RADBEAM2D_TLEFT);
1481  // original top-hat beam
1482  //FTYPE beamshape=(FTYPE)(V[1]>RADBEAM2D_BEAML && V[1]<RADBEAM2D_BEAMR && RADBEAM2D_IFBEAM);
1483  // gaussian-like beam:
1484  FTYPE powbeam=6.0;
1485  FTYPE beamhalfwidth=0.5*(RADBEAM2D_BEAMR-RADBEAM2D_BEAML);
1486  FTYPE beamcenter=(RADBEAM2D_BEAML+RADBEAM2D_BEAMR)*0.5;
1487  FTYPE beamshape=exp(-pow(V[1]-beamcenter,powbeam)/(2.0*pow(beamhalfwidth,powbeam)))*RADBEAM2D_IFBEAM;
1488 
1489 
1490 #define INJECTINFLUIDFRAME 0 // need beam to go out as designed, not with fluid, so usually 0
1491 
1492  if(INJECTINFLUIDFRAME){
1493  //E, F^i in orthonormal fluid frame
1494  FTYPE Fx,Fy,Fz;
1495  // default flux
1496  Fx=Fy=Fz=0;
1497  // beam flux
1498  Fz=RADBEAM2D_NLEFT*ERADINJ;
1499 
1500  FTYPE pradffortho[NPR];
1501  pradffortho[PRAD0] = ERADAMB + ERADINJ*beamshape;
1502  pradffortho[PRAD1] = Fx*beamshape;
1503  pradffortho[PRAD2] = Fy*beamshape;
1504  pradffortho[PRAD3] = Fz*beamshape;
1505 
1506  int whichvel=WHICHVEL; // in which vel U1-U3 set
1507  int whichcoordfluid=PRIMECOORDS; // in which coordinates U1-U3 set
1508  int whichcoordrad=MCOORD; // in which coordinates E,F are orthonormal
1509  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
1510 
1511 
1512 #if(0)
1513  // try using outflow outside of beam
1514  if(){
1515  }
1516  else{
1517  pr[PRAD0]=pr0[PRAD0];
1518  pr[PRAD1]=pr0[PRAD1];
1519  pr[PRAD2]=pr0[PRAD2];
1520  pr[PRAD3]=pr0[PRAD3];
1521  if(pr[PRAD3]>0.0) pr[PRAD3]=0.0; // but don't let radiative inflow
1522  }
1523 #endif
1524 
1525  }
1526  else{
1527 
1528  int whichvel;
1529  whichvel=VEL4;
1530 
1531  // get coordinate basis in VEL4 format
1532  FTYPE uradcon[NDIM],othersrad[NUMOTHERSTATERESULTS];
1533  ucon_calc(&pr[URAD1-U1],ptrgeom[URAD1],uradcon,othersrad);
1534  // get coordinate basis in MCOORD basis
1535  FTYPE uradx,urady,uradz;
1536  uradx=urady=0.0;
1537 
1538  uradx=uradcon[1]*sqrt(fabs(ptrgeom[URAD1]->gcov[GIND(1,1)]))/sqrt(fabs(ptrgeomreal->gcov[GIND(1,1)]));
1539  urady=uradcon[2]*sqrt(fabs(ptrgeom[URAD2]->gcov[GIND(2,2)]))/sqrt(fabs(ptrgeomreal->gcov[GIND(2,2)]));
1540  uradz=uradcon[3]*sqrt(fabs(ptrgeom[URAD3]->gcov[GIND(3,3)]))/sqrt(fabs(ptrgeomreal->gcov[GIND(3,3)]));
1541  if(uradz>0.0) uradz=0.0; // limit so no arbitrary radiative inflow
1542 
1543  // override uradz
1544  if(V[1]>RADBEAM2D_BEAML && V[1]<RADBEAM2D_BEAMR && RADBEAM2D_IFBEAM) uradz=1.0/sqrt(1.0 - RADBEAM2D_NLEFT*RADBEAM2D_NLEFT);
1545  // else uradz=0.0;
1546 
1547  PBOUNDLOOP(pliter,pl) if(pl==URAD0) pr[URAD0] = ERADINJ;
1548  PBOUNDLOOP(pliter,pl) if(pl==URAD1) pr[URAD1] = uradx;
1549  PBOUNDLOOP(pliter,pl) if(pl==URAD2) pr[URAD2] = urady;
1550  PBOUNDLOOP(pliter,pl) if(pl==URAD3) pr[URAD3] = uradz;
1551 
1552  // get all primitives in WHICHVEL/PRIMECOORDS value
1553  primefluid_EVrad_to_primeall(&whichvel, &whichcoord, ptrgeom[RHO],MAC(prim,i,j,k),MAC(prim,i,j,k)); // assumes ptrgeom[RHO] is same location as all other primitives (as is currently true).
1554  }
1555 
1556 
1557  }// end if not staggered field
1558 
1559 
1560  }// end loop over inner i's
1561  } // over block
1562  }// end if correct BC and should be doind BC for this core
1563  }// end parallel region
1564 
1565  return(0);
1566 }
1567 
1568 
1569 
1570 
1571 
1573 int bound_radbeam2dflowinflow(int dir,
1574  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
1575  int *inboundloop,
1576  int *outboundloop,
1577  int *innormalloop,
1578  int *outnormalloop,
1579  int (*inoutlohi)[NUMUPDOWN][NDIM],
1580  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
1581  int *dosetbc,
1582  int enerregion,
1583  int *localenerpos
1584  )
1585 
1586 {
1587 
1588 
1589 #pragma omp parallel // assume don't require EOS
1590  {
1591 
1592  int i,j,k,pl,pliter;
1593  FTYPE vcon[NDIM],X[NDIM],V[NDIM];
1594 #if(WHICHVEL==VEL3)
1595  int failreturn;
1596 #endif
1597  int ri, rj, rk; // reference i,j,k
1598  FTYPE prescale[NPR];
1599  int jj,kk;
1600  struct of_geom geomdontuse[NPR];
1601  struct of_geom *ptrgeom[NPR];
1602  struct of_geom rgeomdontuse[NPR];
1603  struct of_geom *ptrrgeom[NPR];
1604 
1605  // assign memory
1606  PALLLOOP(pl){
1607  ptrgeom[pl]=&(geomdontuse[pl]);
1608  ptrrgeom[pl]=&(rgeomdontuse[pl]);
1609  }
1610 
1611 
1612  extern int RADBEAM2D_BEAMNO;
1613  extern int RADBEAM2D_FLATBACKGROUND;
1614  extern FTYPE RADBEAM2D_RHOAMB;
1615  extern FTYPE RADBEAM2D_TAMB;
1616  extern int RADBEAM2D_BLOB;
1617  extern FTYPE RADBEAM2D_BLOBW;
1618  extern FTYPE RADBEAM2D_BLOBP;
1619  extern FTYPE RADBEAM2D_BLOBX;
1620  extern FTYPE RADBEAM2D_BLOBZ;
1621  extern FTYPE RADBEAM2D_PAR_D;
1622  extern FTYPE RADBEAM2D_PAR_E;
1623  extern int RADBEAM2D_IFBEAM;
1624  extern FTYPE RADBEAM2D_TLEFT;
1625  extern FTYPE RADBEAM2D_NLEFT;
1626  extern FTYPE RADBEAM2D_BEAML;
1627  extern FTYPE RADBEAM2D_BEAMR;
1628 
1629  if(dir==X1UP && BCtype[X1UP]==RADBEAM2DFLOWINFLOW && totalsize[1]>1 && mycpupos[1] == ncpux1-1 ){
1630 
1631  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
1633 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1634  OPENMPBCLOOPBLOCK{
1635  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
1636 
1637 
1638 
1639  ri=riout;
1640  rj=j;
1641  rk=k;
1642 
1643 
1644  // ptrrgeom : i.e. ref geom
1645  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
1646 
1647 
1648  LOOPBOUND1OUT{
1649  FTYPE *pr=&MACP0A1(prim,i,j,k,0);
1650 
1651  //initially copying everything
1652  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
1653 
1654 
1655  if(ispstag==0){ // only do something special with non-field primitives
1656 
1657  // local geom
1658  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
1659 
1660  //coordinates of the ghost cell
1661  bl_coord_ijk_2(i,j,k,CENT,X, V);
1662 
1663  int whichcoord;
1664  whichcoord=MCOORD;
1665 
1666 
1667  FTYPE ERADAMB;
1668  FTYPE rho,uint,Vr;
1669  if(RADBEAM2D_FLATBACKGROUND){
1670  Vr=0.0;
1671  rho=RADBEAM2D_RHOAMB;
1672  uint=calc_PEQ_ufromTrho(RADBEAM2D_TAMB,rho);
1673  ERADAMB=calc_LTE_EfromT(RADBEAM2D_TAMB);
1674  }
1675  else{
1676  //zaczynam jednak od profilu analitycznego:
1677  FTYPE r=V[1];
1678  FTYPE mD=RADBEAM2D_PAR_D/(r*r*sqrt(2./r*(1.-2./r)));
1679  FTYPE mE=RADBEAM2D_PAR_E/(pow(r*r*sqrt(2./r),gamideal)*pow(1.-2./r,(gamideal+1.)/4.));
1680  Vr=sqrt(2./r)*(1.-2./r);
1681 
1682  // get metric grid geometry for these ICs
1683  int getprim=0;
1684  struct of_geom geomrealdontuse;
1685  struct of_geom *ptrgeomreal=&geomrealdontuse;
1686  gset(getprim,whichcoord,i,j,k,ptrgeomreal);
1687 
1688  FTYPE W=1./sqrt(1.-Vr*Vr*ptrgeomreal->gcov[GIND(1,1)]); // assumes RHO location is good for all these quantities
1689  rho=RADBEAM2D_PAR_D/(r*r*sqrt(2./r));
1690  FTYPE T=RADBEAM2D_TAMB;
1691  // FTYPE ERAD=calc_LTE_EfromT(T);
1692  uint=mE/W;
1693  ERADAMB=calc_LTE_Efromurho(uint,rho);
1694  }
1695 
1696  // set quantities at outer radial edge
1697  pr[RHO] = rho;
1698  pr[UU] = uint;
1699  pr[U1] = -Vr;
1700  pr[U2] = 0.;
1701  pr[U3] = 0.;
1702 
1703 
1704  FTYPE ERAD=ERADAMB;
1705 
1706 
1707  if(1){
1708  //E, F^i in orthonormal fluid frame
1709  FTYPE Fx,Fy,Fz;
1710  // default flux
1711  Fx=Fy=Fz=0;
1712 
1713  FTYPE pradffortho[NPR];
1714  pradffortho[PRAD0] = ERAD;
1715  pradffortho[PRAD1] = Fx;
1716  pradffortho[PRAD2] = Fy;
1717  pradffortho[PRAD3] = Fz;
1718 
1719  int whichvel=VEL4;
1720  int whichcoordfluid=whichcoord; // in which coordinates U1-U3 set
1721  int whichcoordrad=whichcoord; // in which coordinates E,F are orthonormal
1722  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
1723 
1724  }
1725  else if(0){
1726 
1727  FTYPE uradx,urady,uradz;
1728  uradx=urady=uradz=0.0;
1729 
1730  pr[PRAD0] = ERAD;
1731  pr[PRAD1] = uradx;
1732  pr[PRAD2] = urady;
1733  pr[PRAD3] = uradz;
1734 
1735  // get all primitives in WHICHVEL/PRIMECOORDS value
1736  int whichvel;
1737  whichvel=VEL4;
1738  if (bl2met2metp2v(whichvel, whichcoord,pr, i,j,k) >= 1){
1739  FAILSTATEMENT("bounds.koral.c:bound_radbeam2dflowinflow()", "bl2ks2ksp2v()", 1);
1740  }
1741  }
1742 
1743 
1744  }// end if not staggered field
1745 
1746 
1747  }// end loop over inner i's
1748  }// over block
1749  }// if correct BC and core
1750  }// end parallel region
1751 
1752  return(0);
1753 }
1754 
1755 
1756 
1757 
1758 
1760 int bound_radatmbeaminflow(int dir,
1761  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
1762  int *inboundloop,
1763  int *outboundloop,
1764  int *innormalloop,
1765  int *outnormalloop,
1766  int (*inoutlohi)[NUMUPDOWN][NDIM],
1767  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
1768  int *dosetbc,
1769  int enerregion,
1770  int *localenerpos
1771  )
1772 
1773 {
1774 
1775 
1776 
1777 
1778 
1779 #pragma omp parallel // assume don't require EOS
1780  {
1781 
1782  int i,j,k,pl,pliter;
1783  FTYPE vcon[NDIM],X[NDIM],V[NDIM];
1784 #if(WHICHVEL==VEL3)
1785  int failreturn;
1786 #endif
1787  int ri, rj, rk; // reference i,j,k
1788  FTYPE prescale[NPR];
1789  int jj,kk;
1790  struct of_geom geomdontuse[NPR];
1791  struct of_geom *ptrgeom[NPR];
1792  struct of_geom rgeomdontuse[NPR];
1793  struct of_geom *ptrrgeom[NPR];
1794 
1795  // assign memory
1796  PALLLOOP(pl){
1797  ptrgeom[pl]=&(geomdontuse[pl]);
1798  ptrrgeom[pl]=&(rgeomdontuse[pl]);
1799  }
1800 
1801 
1802  extern FTYPE RADATM_MDOTEDD;
1803  extern FTYPE RADATM_LUMEDD;
1804  extern int RADATM_THINRADATM;
1805  extern FTYPE RADATM_FERATIO;
1806  extern FTYPE RADATM_FRATIO;
1807  extern FTYPE RADATM_RHOAMB;
1808  extern FTYPE RADATM_TAMB;
1809 
1810 
1811  FTYPE MINX=Rin_array[1];
1812  FTYPE kappaesperrho=calc_kappaes_user(1,0, 0,0,0);
1813  FTYPE FLUXLEFT=RADATM_FRATIO/kappaesperrho/pow(MINX,2.0);
1814 
1815  //at boundary
1816  FTYPE f = (FTYPE)kappaesperrho*FLUXLEFT*MINX*MINX;
1817 
1818  FTYPE p0=RADATM_RHOAMB*RADATM_TAMB;
1819  FTYPE KKK=p0/pow(RADATM_RHOAMB,gamideal);
1820  FTYPE C3=gamideal*KKK/(gamideal-1.)*pow(RADATM_RHOAMB,gamideal-1.)-(1.-f)*(1./MINX+0.*1./MINX/MINX+0.*4./3./MINX/MINX/MINX);
1821 
1822  // dualfprintf(fail_file,"IT: %g %g %g : %g : %g %g %g\n",MINX,kappaesperrho,FLUXLEFT,f,p0,KKK,C3);
1823 
1824  if(dir==X1DN && BCtype[X1DN]==RADATMBEAMINFLOW && (totalsize[1]>1) && (mycpupos[1] == 0) ){
1825 
1826 
1827  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
1829 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1830  OPENMPBCLOOPBLOCK{
1831  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
1832 
1833 
1834  ri=riin;
1835  rj=j;
1836  rk=k;
1837 
1838 
1839  // ptrrgeom : i.e. ref geom
1840  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
1841 
1842  FTYPE *pr;
1843  LOOPBOUND1IN{
1844 
1845  pr = &MACP0A1(prim,i,j,k,0);
1846 
1847  //initially copying everything
1848  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
1849 
1850  if(ispstag==0){
1851  // local PRIMECOORDS geom
1852  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
1853 
1854  //coordinates of the ghost cell
1855  bl_coord_ijk_2(i,j,k,CENT,X, V);
1856 
1857  FTYPE xx,yy,zz,rsq;
1858  coord(i, j, k, CENT, X);
1859  bl_coord(X, V);
1860  xx=V[1];
1861  yy=V[2];
1862  zz=V[3];
1863 
1864 
1865  FTYPE rho=pow((gamideal-1.0)/gamideal/KKK*(C3+(1.-f)*(1./xx+0.*1./xx/xx+0.*4./3./xx/xx/xx)),1./(gamideal-1.0));
1866 
1867  FTYPE pre=KKK*pow(rho,gamideal);
1868 
1869  FTYPE uint=pre/(gamideal-1.0);
1870 
1871  FTYPE Fz=0;
1872  FTYPE Fy=0.;
1873  FTYPE Fx=FLUXLEFT*(MINX/xx)*(MINX/xx);
1874 
1875  FTYPE ERAD;
1876  if(RADATM_THINRADATM){
1877  ERAD=Fx/RADATM_FERATIO;
1878  }
1879  else{
1880  ERAD=calc_LTE_EfromT(calc_PEQ_Tfromurho(uint,rho));
1881  }
1882 
1883  // dualfprintf(fail_file,"BC: i=%d j=%d rho=%g Trad=%g uint=%g ERAD=%g\n",i,j,rho,uint,ERAD);
1884 
1885 
1886  pr[RHO] = rho;
1887  pr[UU] = uint;
1888  pr[U1] = 0.0; // static
1889  pr[U2] = 0.0;
1890  pr[U3] = 0.0;
1891 
1892  //E, F^i in orthonormal fluid frame
1893  FTYPE pradffortho[NPR];
1894  pradffortho[PRAD0] = ERAD;
1895  pradffortho[PRAD1] = Fx;
1896  pradffortho[PRAD2] = Fy;
1897  pradffortho[PRAD3] = Fz;
1898 
1899 
1900  int whichvel=VEL4; // in which vel U1-U3 set
1901  int whichcoordfluid=MCOORD; // in which coordinates U1-U3 set
1902  int whichcoordrad=whichcoordfluid; // in which coordinates E,F are orthonormal
1903  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
1904 
1905  // PLOOP(pliter,pl) dualfprintf(fail_file,"ijk=%d %d %d pl=%d pr=%g\n",i,j,k,pl,pr[pl]);
1906 
1907  }// end if not staggered fields
1908 
1909  }// end loop over inner i's
1910 
1911  }// end over loop
1912  }// end if correct boundary condition and core
1913 
1914 
1915 
1916 
1917  if(dir==X1UP && BCtype[X1UP]==RADATMBEAMINFLOW && (totalsize[1]>1) && (mycpupos[1] == ncpux1-1) ){
1918 
1919 
1920  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
1922 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
1923  OPENMPBCLOOPBLOCK{
1924  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
1925 
1926 
1927  ri=riout;
1928  rj=j;
1929  rk=k;
1930 
1931 
1932  // ptrrgeom : i.e. ref geom
1933  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
1934 
1935  FTYPE *pr;
1936  LOOPBOUND1OUT{
1937 
1938  pr = &MACP0A1(prim,i,j,k,0);
1939 
1940 
1941  //initially copying everything
1942  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
1943 
1944  if(ispstag==0){
1945  // local geom
1946  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
1947 
1948  //coordinates of the ghost cell
1949  bl_coord_ijk_2(i,j,k,CENT,X, V);
1950 
1951 
1952  FTYPE xx,yy,zz,rsq;
1953  coord(i, j, k, CENT, X);
1954  bl_coord(X, V);
1955  xx=V[1];
1956  yy=V[2];
1957  zz=V[3];
1958 
1959 
1960 
1961  FTYPE rho=pow((gamideal-1.0)/gamideal/KKK*(C3+(1.-f)*(1./xx+0.*1./xx/xx+0.*4./3./xx/xx/xx)),1./(gamideal-1.0));
1962 
1963  FTYPE pre=KKK*pow(rho,gamideal);
1964 
1965  FTYPE uint=pre/(gamideal-1.0);
1966 
1967  FTYPE Fz=0;
1968  FTYPE Fy=0.;
1969  FTYPE Fx=FLUXLEFT*(MINX/xx)*(MINX/xx);
1970 
1971  FTYPE ERAD;
1972  if(RADATM_THINRADATM){
1973  ERAD=Fx/RADATM_FERATIO;
1974  }
1975  else{
1976  ERAD=calc_LTE_EfromT(calc_PEQ_Tfromurho(uint,rho));
1977  }
1978 
1979 
1980  pr[RHO] = rho;
1981  pr[UU] = uint;
1982  pr[U1] = 0.0; // static
1983  pr[U2] = 0.0;
1984  pr[U3] = 0.0;
1985 
1986  //E, F^i in orthonormal fluid frame
1987  FTYPE pradffortho[NPR];
1988  pradffortho[PRAD0] = ERAD;
1989  pradffortho[PRAD1] = Fx;
1990  pradffortho[PRAD2] = Fy;
1991  pradffortho[PRAD3] = Fz;
1992 
1993 
1994  int whichvel=VEL4;
1995  int whichcoordfluid=MCOORD;
1996  int whichcoordrad=whichcoordfluid;
1997  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
1998 
1999  }// end if not staggered fields
2000 
2001  }// end loop over outer i's
2002 
2003  }// end over loop
2004  }// end if correct boundary condition and core
2005 
2006 
2007 
2008  }// end parallel region
2009 
2010 
2011 
2012 
2013 
2014 
2015 
2016  return(0);
2017 }
2018 
2019 
2020 
2021 
2022 
2024 int bound_radwallinflow(int dir,
2025  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
2026  int *inboundloop,
2027  int *outboundloop,
2028  int *innormalloop,
2029  int *outnormalloop,
2030  int (*inoutlohi)[NUMUPDOWN][NDIM],
2031  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
2032  int *dosetbc,
2033  int enerregion,
2034  int *localenerpos
2035  )
2036 
2037 {
2038 
2039 
2040 #pragma omp parallel // assume don't require EOS
2041  {
2042 
2043  FTYPE LTEFACTOR=1.;
2044  FTYPE URFX=100.;
2045  FTYPE URFY=-30.;
2046 
2047 
2048  int i,j,k,pl,pliter;
2049  FTYPE vcon[NDIM],X[NDIM],V[NDIM];
2050 #if(WHICHVEL==VEL3)
2051  int failreturn;
2052 #endif
2053  int ri, rj, rk; // reference i,j,k
2054  FTYPE prescale[NPR];
2055  int jj,kk;
2056  struct of_geom geomdontuse[NPR];
2057  struct of_geom *ptrgeom[NPR];
2058  struct of_geom rgeomdontuse[NPR];
2059  struct of_geom *ptrrgeom[NPR];
2060 
2061  // assign memory
2062  PALLLOOP(pl){
2063  ptrgeom[pl]=&(geomdontuse[pl]);
2064  ptrrgeom[pl]=&(rgeomdontuse[pl]);
2065  }
2066 
2067 
2068 
2069  if(dir==X1DN && BCtype[X1DN]==RADWALLINFLOW && (totalsize[1]>1) && (mycpupos[1] == 0) ){
2070 
2071 
2072 
2073  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
2075 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2076  OPENMPBCLOOPBLOCK{
2077  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
2078 
2079 
2080 
2081  ri=riin;
2082  rj=j;
2083  rk=k;
2084 
2085 
2086  // ptrrgeom : i.e. ref geom
2087  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
2088 
2089 
2090  LOOPBOUND1IN{
2091 
2092 
2093  //initially copying everything
2094  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
2095 
2096 
2097  if(ispstag==0){ // only do something special with non-field primitives
2098 
2099  // local geom
2100  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
2101 
2102  //coordinates of the ghost cell
2103  bl_coord_ijk_2(i,j,k,CENT,X, V);
2104 
2105  int whichvel;
2106  whichvel=VEL4;
2107  int whichcoord;
2108  whichcoord=CARTMINKMETRIC2;
2109 
2110  PBOUNDLOOP(pliter,pl){
2111  if(pl==RHO) MACP0A1(prim,i,j,k,RHO) = 1.0;
2112  if(pl==UU) MACP0A1(prim,i,j,k,UU) = 1.0;
2113  if(pl==U1) MACP0A1(prim,i,j,k,U1) = 0.0;
2114 
2115 
2116  if(V[2]>0.3){
2117  if(pl==URAD0) MACP0A1(prim,i,j,k,URAD0) = pow(100,4.0);
2118  if(pl==URAD1) MACP0A1(prim,i,j,k,URAD1) = URFX;
2119  if(pl==URAD2) MACP0A1(prim,i,j,k,URAD2) = URFY;
2120  if(pl==URAD3) MACP0A1(prim,i,j,k,URAD3) = 0.;
2121  }
2122  else{
2123  if(pl==URAD0) MACP0A1(prim,i,j,k,URAD0) = 1.0;
2124  if(pl==URAD1) MACP0A1(prim,i,j,k,URAD1) = 0.0;
2125  if(pl==URAD2) MACP0A1(prim,i,j,k,URAD2) = 0.;
2126  if(pl==URAD3) MACP0A1(prim,i,j,k,URAD3) = 0.;
2127  }
2128  }// end over pl's allowed for bounding
2129 
2130  // get all primitives in WHICHVEL/PRIMECOORDS value
2131  if(1) primefluid_EVrad_to_primeall(&whichvel, &whichcoord, ptrgeom[RHO],MAC(prim,i,j,k),MAC(prim,i,j,k)); // assumes ptrgeom[RHO] is same location as all other primitives (as is currently true).
2132 
2133  }// end if not staggered field
2134 
2135 
2136  }// end loop over inner i's
2137  }
2138  }
2139 
2140 
2141 
2142 
2143  if(dir==X2UP && BCtype[X2UP]==RADWALLINFLOW && (totalsize[2]>1) && (mycpupos[2] == ncpux2-1) ){
2144 
2145 
2146 
2147  OPENMPBCLOOPVARSDEFINELOOPX2DIR; OPENMPBCLOOPSETUPLOOPX2DIR;
2149 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2150  OPENMPBCLOOPBLOCK{
2151  OPENMPBCLOOPBLOCK2IJKLOOPX2DIR(i,k);
2152 
2153 
2154 
2155  ri=i;
2156  rj=rjout;
2157  rk=k;
2158 
2159 
2160  // ptrrgeom : i.e. ref geom
2161  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
2162 
2163 
2164  LOOPBOUND2OUT{
2165 
2166 
2167  //initially copying everything
2168  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
2169 
2170 
2171  if(ispstag==0){ // only do something special with non-field primitives
2172 
2173  // local geom
2174  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
2175 
2176  //coordinates of the ghost cell
2177  bl_coord_ijk_2(i,j,k,CENT,X, V);
2178 
2179  int whichvel;
2180  whichvel=VEL4;
2181  int whichcoord;
2182  whichcoord=CARTMINKMETRIC2;
2183 
2184  PBOUNDLOOP(pliter,pl){
2185  if(pl==RHO) MACP0A1(prim,i,j,k,RHO) = 1.0;
2186  if(pl==UU) MACP0A1(prim,i,j,k,UU) = 1.0;
2187  if(pl==U1) MACP0A1(prim,i,j,k,U1) = 0.0;
2188 
2189  if(1){
2190  if(pl==URAD0) MACP0A1(prim,i,j,k,URAD0) = pow(100,4.0);
2191  if(pl==URAD1) MACP0A1(prim,i,j,k,URAD1) = URFX;
2192  if(pl==URAD2) MACP0A1(prim,i,j,k,URAD2) = URFY;
2193  if(pl==URAD3) MACP0A1(prim,i,j,k,URAD3) = 0.;
2194  }
2195 
2196  }// over allowed pl's
2197 
2198 
2199  // get all primitives in WHICHVEL/PRIMECOORDS value
2200  if(1) primefluid_EVrad_to_primeall(&whichvel, &whichcoord, ptrgeom[RHO],MAC(prim,i,j,k),MAC(prim,i,j,k)); // assumes ptrgeom[RHO] is same location as all other primitives (as is currently true).
2201 
2202 
2203  }// end if not staggered field
2204 
2205 
2206  }// end loop over outer j's
2207  }
2208  }
2209 
2210 
2211 
2212 
2213 
2214 
2215 
2216 
2217  }// end parallel region
2218 
2219  return(0);
2220 }
2221 
2222 
2223 
2224 
2225 
2226 
2227 
2228 
2229 
2231 int bound_radbondiinflow(int dir,
2232  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
2233  int *inboundloop,
2234  int *outboundloop,
2235  int *innormalloop,
2236  int *outnormalloop,
2237  int (*inoutlohi)[NUMUPDOWN][NDIM],
2238  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
2239  int *dosetbc,
2240  int enerregion,
2241  int *localenerpos
2242  )
2243 
2244 {
2245 
2246 
2247 #pragma omp parallel // assume don't require EOS
2248  {
2249 
2250  int i,j,k,pl,pliter;
2251  FTYPE vcon[NDIM],X[NDIM],V[NDIM];
2252 #if(WHICHVEL==VEL3)
2253  int failreturn;
2254 #endif
2255  int ri, rj, rk; // reference i,j,k
2256  FTYPE prescale[NPR];
2257  int jj,kk;
2258  struct of_geom geomdontuse[NPR];
2259  struct of_geom *ptrgeom[NPR];
2260  struct of_geom rgeomdontuse[NPR];
2261  struct of_geom *ptrrgeom[NPR];
2262 
2263  // assign memory
2264  PALLLOOP(pl){
2265  ptrgeom[pl]=&(geomdontuse[pl]);
2266  ptrrgeom[pl]=&(rgeomdontuse[pl]);
2267  }
2268 
2269 
2270 
2271  if(dir==X1UP && BCtype[X1UP]==RADBONDIINFLOW && totalsize[1]>1 && mycpupos[1] == ncpux1-1 ){
2272 
2273 
2274  extern FTYPE RADBONDI_TESTNO;
2275  extern FTYPE RADBONDI_PRADGAS;
2276  extern FTYPE RADBONDI_TGAS0;
2277  extern FTYPE RADBONDI_MDOTPEREDD;
2278  extern FTYPE RADBONDI_MDOTEDD;
2279  extern FTYPE RADBONDI_MINX;
2280  extern FTYPE RADBONDI_MAXX;
2281 
2282 
2283  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
2285 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2286  OPENMPBCLOOPBLOCK{
2287  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
2288 
2289 
2290 
2291  ri=riout;
2292  rj=j;
2293  rk=k;
2294 
2295 
2296  // ptrrgeom : i.e. ref geom
2297  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
2298 
2299 
2300  FTYPE *pr;
2301  LOOPBOUND1OUT{
2302 
2303  pr = &MACP0A1(prim,i,j,k,0);
2304 
2305 
2306  //initially copying everything
2307  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
2308 
2309 
2310  if(ispstag==0){ // only do something special with non-field primitives
2311 
2312  // local geom
2313  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
2314 
2315  //coordinates of the ghost cell
2316  bl_coord_ijk_2(i,j,k,CENT,X, V);
2317 
2318  // Identical to init, so could use analytic, but don't for now.
2319 
2320  FTYPE xx,yy,zz,rsq;
2321  xx=V[1];
2322  yy=V[2];
2323  zz=V[3];
2324 
2325  int whichvel=VEL3; // in which vel U1-U3 set
2326  int whichcoordfluid=MCOORD; // in which coordinates U1-U3 set
2327  // get metric grid geometry for these ICs
2328  int getprim=0;
2329  struct of_geom geomrealdontuse;
2330  struct of_geom *ptrgeomreal=&geomrealdontuse;
2331  gset(getprim,whichcoordfluid,i,j,k,ptrgeomreal);
2332 
2333 
2334  FTYPE rho,ERAD,uint;
2335  FTYPE rho0,Tgas0,ur,Tgas,Trad,r,rcm,prad,pgas,vx,ut;
2336 
2337  FTYPE Fx,Fy,Fz;
2338  Fx=Fy=Fz=0;
2339 
2340  //at outern boundary
2341  r=RADBONDI_MAXX;
2342  ur=-sqrt(2./r);
2343  rho0=-RADBONDI_MDOTPEREDD*RADBONDI_MDOTEDD/(4.*Pi*r*r*ur);
2344  Tgas0=RADBONDI_TGAS0;
2345 
2346  //at given cell
2347  r=xx;
2348  ur=-sqrt(2./r);
2349  ut=sqrt((-1.-ur*ur*ptrgeomreal->gcov[GIND(1,1)])/ptrgeomreal->gcov[GIND(0,0)]);
2350  vx=ur/ut;
2351  rho=-RADBONDI_MDOTPEREDD*RADBONDI_MDOTEDD/(4.*Pi*r*r*ur);
2352  Tgas=Tgas0*pow(rho/rho0,gam-1.);
2353 
2354  uint=calc_PEQ_ufromTrho(Tgas,rho);
2355 
2356  pgas=rho*Tgas;
2357  prad=RADBONDI_PRADGAS*pgas;
2358  ERAD=prad*3.;
2359 
2360 
2361  pr[RHO] = rho ;
2362  pr[UU] = uint;
2363  pr[U1] = vx;
2364  pr[U2] = 0 ;
2365  pr[U3] = 0 ;
2366 
2367  //E, F^i in orthonormal fluid frame
2368  FTYPE pradffortho[NPR];
2369  pradffortho[PRAD0] = ERAD;
2370  pradffortho[PRAD1] = Fx;
2371  pradffortho[PRAD2] = Fy;
2372  pradffortho[PRAD3] = Fz;
2373 
2374  // dualfprintf(fail_file,"rho=%g uint=%g vx=%g ERAD=%g\n",rho,uint,vx,ERAD);
2375 
2376  int whichcoordrad=whichcoordfluid; // in which coordinates E,F are orthonormal
2377  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
2378 
2379  // PLOOP(pliter,pl) dualfprintf(fail_file,"pl=%d pr=%g\n",pl,pr[pl]);
2380 
2381  }// end if not staggered field
2382 
2383 
2384  }// end loop over inner i's
2385  }// over block
2386  }// if correct BC and core
2387  }// end parallel region
2388 
2389  return(0);
2390 }
2391 
2392 
2393 
2394 
2395 
2396 
2398 int bound_raddot(
2399  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
2400  int *inboundloop,
2401  int *outboundloop,
2402  int *innormalloop,
2403  int *outnormalloop,
2404  int (*inoutlohi)[NUMUPDOWN][NDIM],
2405  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
2406  int *dosetbc,
2407  int enerregion,
2408  int *localenerpos
2409  )
2410 
2411 {
2412 
2413 
2414 #pragma omp parallel // assume don't require EOS
2415  {
2416 
2417  extern FTYPE RADDOT_XDOT;
2418  extern FTYPE RADDOT_YDOT;
2419  extern FTYPE RADDOT_ZDOT;
2420  extern int RADDOT_IDOT;
2421  extern int RADDOT_JDOT;
2422  extern int RADDOT_KDOT;
2423  extern FTYPE RADDOT_FYDOT;
2424  extern FTYPE RADDOT_LTEFACTOR;
2425  extern FTYPE RADDOT_URFX;
2426  extern FTYPE RADDOT_F1;
2427  extern FTYPE RADDOT_F2;
2428 
2429 
2430  int i,j,k,pl,pliter;
2431  FTYPE X[NDIM],V[NDIM];
2432  int jj,kk;
2433  struct of_geom geomdontuse[NPR];
2434  struct of_geom *ptrgeom[NPR];
2435  struct of_geom rgeomdontuse[NPR];
2436  struct of_geom *ptrrgeom[NPR];
2437 
2438  // assign memory
2439  PALLLOOP(pl){
2440  ptrgeom[pl]=&(geomdontuse[pl]);
2441  ptrrgeom[pl]=&(rgeomdontuse[pl]);
2442  }
2443 
2444  OPENMP3DLOOPVARSDEFINE; OPENMP3DLOOPSETUPFULL;
2445 
2446 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize)) nowait // can nowait since each fluxvec[dir] is set separately
2447  OPENMP3DLOOPBLOCK{
2448  OPENMP3DLOOPBLOCK2IJK(i,j,k);
2450 
2451  // KORALTODO: Koral has different IC and BC for the dot
2452  if(startpos[1]+i==RADDOT_IDOT && startpos[2]+j==RADDOT_JDOT && startpos[3]+k==RADDOT_KDOT){
2453  // NOTEMARK: It's normal for inversion to fail right on the dot, but doesn't affect evolution since the dot's primitives are overwritten and never will the updated conserved quantities for the dot be used.
2454 
2455  if(ispstag==0){ // only do something special with non-field primitives
2456 
2457  FTYPE *pr=&MACP0A1(prim,i,j,k,0);
2458 
2459  // get current location
2460  bl_coord_ijk_2(i,j,k,CENT,X, V);
2461 
2462  // local geom
2463  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
2464 
2465  pr[RHO] = 1.0 ;
2466  pr[UU] = 1.0;
2467  pr[U1] = 0.0;
2468  pr[U2] = 0 ;
2469  pr[U3] = 0 ;
2470 
2471  //E, F^i in orthonormal fluid frame
2472  FTYPE pradffortho[NPR];
2473  pradffortho[PRAD0] = RADDOT_LTEFACTOR*calc_LTE_Efromurho(pr[RHO],pr[UU]);
2474  pradffortho[PRAD1] = 0;
2475  pradffortho[PRAD2] = 0;
2476  pradffortho[PRAD3] = 0;
2477 
2478  // dualfprintf(fail_file,"GOT BC DOT: nstep=%ld steppart=%d\n",nstep,steppart);
2479  if(N1==1) pradffortho[PRAD0] *= RADDOT_F1;
2480  else{
2481  pradffortho[PRAD0]*=RADDOT_F2;
2482  pradffortho[PRAD2]=RADDOT_FYDOT*pradffortho[PRAD0];
2483  }
2484 
2485  int whichvel=VEL4; // in which vel U1-U3 set
2486  int whichcoordfluid=MCOORD; // in which coordinates U1-U3 set
2487  int whichcoordrad=whichcoordfluid; // in which coordinates E,F are orthonormal
2488  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
2489 
2490  }// end if DOT
2491  }// end if not staggered field
2492  }// end loop over zones
2493 
2494 
2495 
2496  }// end parallel region
2497 
2498  return(0);
2499 }
2500 
2501 
2502 
2503 
2504 
2505 
2507 int bound_radnt(int dir,
2508  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
2509  int *inboundloop,
2510  int *outboundloop,
2511  int *innormalloop,
2512  int *outnormalloop,
2513  int (*inoutlohi)[NUMUPDOWN][NDIM],
2514  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
2515  int *dosetbc,
2516  int enerregion,
2517  int *localenerpos
2518  )
2519 
2520 {
2521 
2522 
2523 #pragma omp parallel // assume don't require EOS
2524  {
2525 
2526  extern FTYPE RADNT_MINX;
2527  extern FTYPE RADNT_MAXX;
2528  extern FTYPE RADNT_KKK;
2529  extern FTYPE RADNT_ELL;
2530  extern FTYPE RADNT_UTPOT;
2531  extern FTYPE RADNT_RHOATMMIN;
2532  extern FTYPE RADNT_RHODONUT;
2533  extern FTYPE RADNT_UINTATMMIN;
2534  extern FTYPE RADNT_ERADATMMIN;
2535  extern FTYPE RADNT_DONUTTYPE;
2536  extern FTYPE RADNT_INFLOWING;
2537  extern FTYPE RADNT_TGASATMMIN;
2538  extern FTYPE RADNT_TRADATMMIN;
2539  extern FTYPE RADNT_ROUT;
2540  extern FTYPE RADNT_OMSCALE;
2541  extern FTYPE RADNT_FULLPHI;
2542  extern FTYPE RADNT_DONUTRADPMAX;
2543  extern FTYPE RADNT_HOVERR;
2544  extern FTYPE RADNT_LPOW;
2545 
2546 
2547  int i,j,k,pl,pliter;
2548  FTYPE vcon[NDIM],X[NDIM],V[NDIM];
2549 #if(WHICHVEL==VEL3)
2550  int failreturn;
2551 #endif
2552  int ri, rj, rk; // reference i,j,k
2553  FTYPE prescale[NPR];
2554  int jj,kk;
2555  struct of_geom geomdontuse[NPR];
2556  struct of_geom *ptrgeom[NPR];
2557  struct of_geom rgeomdontuse[NPR];
2558  struct of_geom *ptrrgeom[NPR];
2559 
2560  // assign memory
2561  PALLLOOP(pl){
2562  ptrgeom[pl]=&(geomdontuse[pl]);
2563  ptrrgeom[pl]=&(rgeomdontuse[pl]);
2564  }
2565 
2566 
2567 
2568  if(dir==X1UP && BCtype[X1UP]==RADNTBC && totalsize[1]>1 && mycpupos[1] == ncpux1-1 ){
2569 
2570  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
2572 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2573  OPENMPBCLOOPBLOCK{
2574  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
2575 
2576 
2577 
2578  ri=riout;
2579  rj=j;
2580  rk=k;
2581 
2582 
2583  // ptrrgeom : i.e. ref geom
2584  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
2585 
2586 
2587  FTYPE *pr;
2588  LOOPBOUND1OUT{
2589  pr = &MACP0A1(prim,i,j,k,0);
2590 
2591  //initially copying everything
2592  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
2593 
2594 
2595  if(ispstag==0){ // only do something special with non-field primitives
2596 
2597  // local geom
2598  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
2599 
2600  FTYPE r,th,ph;
2601  coord(i, j, k, CENT, X);
2602  bl_coord(X, V);
2603  r=V[1];
2604  th=V[2];
2605  ph=V[3];
2606 
2607  // Identical to IC, except more involved check for inflow vs. outflow
2608 
2609  int whichcoord;
2610  if(WHICHPROBLEM==RADFLATDISK) whichcoord=MCOORD; // whatever else
2611  else whichcoord=BLCOORDS; // want to setup things in BLCOORDS
2612  int whichvel=VEL4;
2613 
2614  // get metric grid geometry for these ICs
2615  int getprim=0;
2616  struct of_geom geomrealdontuse;
2617  struct of_geom *ptrgeomreal=&geomrealdontuse;
2618  gset(getprim,whichcoord,i,j,k,ptrgeomreal);
2619 
2620  FTYPE uconlab[NDIM];
2621  FTYPE others[NUMOTHERSTATERESULTS];
2622 
2623  // see if fluid flow wants to go in or out
2624  // get PRIMECOORDS ucon
2625  ucon_calc(pr,ptrgeom[U1],uconlab,others);
2626  // get MCOORD
2627  metptomet(i,j,k,uconlab);
2628  // get whichcoord
2629  coordtrans(MCOORD,whichcoord,i,j,k,CENT,uconlab);
2630  if(uconlab[RR]<=0.0){ // check in whichvel whichcoord
2631  pr[RHO]=RADNT_RHOATMMIN*pow(r/RADNT_ROUT,-1.5);
2632  pr[UU]=RADNT_UINTATMMIN*pow(r/RADNT_ROUT,-2.5);
2633  set_zamo_velocity(whichvel,ptrgeomreal,pr); // only sets U1-U3 to zamo
2634  }
2635  else{
2636  uconlab[RR]=0.0;
2637  // overwrite pr[U1-U3] with this non-radially moving flow in whichvel whichcoord version
2638  ucon2pr(whichvel,uconlab,ptrgeomreal,pr);
2639  }
2640 
2641 
2642  if(pr[PRAD1]<0.0){ // if active cell (now ghost cell) inflows, set values, keep outflow value
2643  if(1){
2644  // assume radiation in fluid frame with zero flux
2645  FTYPE pradffortho[NPR];
2646 
2647  pradffortho[PRAD0] = RADNT_ERADATMMIN; // fluid frame ortho!
2648  pradffortho[PRAD1] = 0;
2649  pradffortho[PRAD2] = 0;
2650  pradffortho[PRAD3] = 0;
2651 
2652  int whichcoordfluid;
2653  whichcoordfluid=whichcoord;
2654  int whichcoordrad=whichcoordfluid; // in which coordinates E,F are orthonormal
2655  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
2656  }
2657  else{
2658  FTYPE gammamax=10.;
2659  FTYPE rout=2.; //normalize at r_BL=2
2660  // RADNT_ERADATMMIN here, as in koral, assumes in ut frame, so not really consistent with initial conditions
2661  pr[PRAD0]=RADNT_ERADATMMIN*(rout/r)*(rout/r)*(rout/r)*(rout/r);
2662 
2663  FTYPE ut[NDIM]={0.,-gammamax*pow(r/rout,1.),0.,0.}; // assume in BLCOORDS 4-vel or SPCMINKMETRIC if no gravity
2664  SLOOPA(jj) pr[URAD1+jj-1]=ut[jj];
2665 
2666  // get all primitives in WHICHVEL/PRIMECOORDS value
2667  if (bl2met2metp2v(whichvel, whichcoord,pr, i,j,k) >= 1){
2668  FAILSTATEMENT("bounds.koral.c:bound_radnt()", "bl2ks2ksp2v()", 1);
2669  }
2670  }
2671  }
2672  else{
2673  // keep outflow value already set
2674  }
2675 
2676 
2677  }// end if not staggered field
2678 
2679 
2680  }// end loop over inner i's
2681  }// over block
2682  }// if correct BC and core
2683 
2684 
2685 
2686 
2687 
2688 
2689 
2690  if(dir==X2UP && BCtype[X2UP]==RADNTBC && (totalsize[2]>1) && (mycpupos[2] == ncpux2-1) ){
2691 
2692  OPENMPBCLOOPVARSDEFINELOOPX2DIR; OPENMPBCLOOPSETUPLOOPX2DIR;
2694 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2695  OPENMPBCLOOPBLOCK{
2696  OPENMPBCLOOPBLOCK2IJKLOOPX2DIR(i,k);
2697 
2698  ri=i;
2699  rj=rjout;
2700  rk=k;
2701 
2702  // ptrrgeom : i.e. ref geom
2703  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
2704 
2705  FTYPE *pr;
2706  LOOPBOUND2OUT{
2707  pr = &MACP0A1(prim,i,j,k,0);
2708 
2709  // ASYMM already done before got here, so only change what's necessary to change
2710 
2711  if(ispstag==0){ // only do something special with non-field primitives
2712 
2713  // local geom
2714  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
2715 
2716  //coordinates of the ghost cell
2717  bl_coord_ijk_2(i,j,k,CENT,X, V);
2718  FTYPE r;
2719  r=V[1];
2720 
2721  FTYPE rin,rout;
2722  int conddisk;
2723  if(WHICHPROBLEM==RADFLATDISK){
2724  rin=15.;
2725  rout=25.;
2726  conddisk=(r<rout); // as in new koral
2727  }
2728  else{
2729  rin=6.;
2730  rout=1E10;
2731  conddisk=(r>rin);
2732  }
2733 
2734  //hot boundary
2735  if(conddisk){
2736 
2737  //E, F^i in orthonormal fluid frame
2738  FTYPE pradffortho[NPR];
2739  if(WHICHPROBLEM==RADFLATDISK) pradffortho[PRAD0] = calc_LTE_EfromT(1.e11/TEMPBAR);
2740  else pradffortho[PRAD0] = calc_LTE_EfromT(1.e11/TEMPBAR)*(1.-sqrt(rin/r))/pow(r,3.);
2741  // KORALTODO: in reality, can only constrain magnitude, not direction, of outflow away from plane.
2742  // KORALTODO: Also, in reality, can't really set vrad=0, else solution dominated by numerical diffusion.
2743  pradffortho[PRAD1] = 0;
2744  if(WHICHPROBLEM==RADFLATDISK) pradffortho[PRAD2] = 0.0; // current koral value
2745  else pradffortho[PRAD2] = -0.5*pradffortho[PRAD0];
2746  pradffortho[PRAD3] = 0;
2747 
2748 
2749  //pr[RHO] and pr[UU] remain same as from ASYMM condition as well as any field
2750  //Keplerian gas with no inflow or outflow
2751  // KORALTODO: in reality, can only constrain magnitude, not direction, of outflow away from plane. But since magnitude is chosen to be zero, then no issue here.
2752  pr[U1]=pr[U2]=0.0; // have to be careful with this for VEL3 (must have rin>>rergo).
2753  if(WHICHPROBLEM==RADFLATDISK) pr[U3]=0.0; // current koral value
2754  else pr[U3]=1./(a + pow(r,1.5));
2755 
2756  int whichvel;
2757  whichvel=VEL3; // VEL3 so can set Keplerian rotation rate
2758 
2759  int whichcoordfluid;
2760  if(WHICHPROBLEM==RADFLATDISK) whichcoordfluid=MCOORD; // whatever else
2761  else whichcoordfluid=BLCOORDS; // want to setup things in BLCOORDS
2762 
2763  // dualfprintf(fail_file,"FLAT: rho=%g uint=%g prad0=%g Eff=%g\n",pr[RHO],pr[UU]/pr[RHO],pr[PRAD0]/pr[RHO],pradffortho[PRAD0]/pr[RHO]);
2764 
2765  int whichcoordrad=whichcoordfluid; // in which coordinates E,F are orthonormal
2766  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
2767 
2768  // dualfprintf(fail_file,"FLATPOST: rho=%g uint=%g Eff=%g\n",pr[RHO],pr[UU]/pr[RHO],pr[PRAD0]/pr[RHO]);
2769  // dualfprintf(fail_file,"opacity : %g\n",pr[RHO]*KAPPA_ES_CODE(rho,T)/1E14*0.1);
2770  // dualfprintf(fail_file,"LBAR: %g\n",LBAR);
2771 
2772  } // end if actually doing something to boundary cells in "hot" boundary
2773 
2774  }// end if not staggered field
2775 
2776 
2777  }// end loop over outer j's
2778  }
2779  }
2780 
2781 
2782 
2783 
2784 
2785 
2786 
2787  }// end parallel region
2788 
2789  return(0);
2790 }
2791 
2792 
2793 
2794 
2795 
2796 
2797 
2799 int bound_x1dn_cylaxis(
2800  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
2801  int *inboundloop,
2802  int *outboundloop,
2803  int *innormalloop,
2804  int *outnormalloop,
2805  int (*inoutlohi)[NUMUPDOWN][NDIM],
2806  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
2807  int *dosetbc,
2808  int enerregion,
2809  int *localenerpos
2810  )
2811 {
2812 
2813 
2814 #pragma omp parallel // assume don't require EOS
2815  {
2816  int i,j,k,pl,pliter;
2817  FTYPE vcon[NDIM]; // coordinate basis vcon
2818 #if(WHICHVEL==VEL3)
2819  int failreturn;
2820 #endif
2821  int ri, rj, rk; // reference i,j,k
2822  FTYPE prescale[NPR];
2823  int jj,kk;
2824 
2825  extern FTYPE RADNT_OMSCALE;
2826  extern FTYPE RADNT_FULLPHI;
2827 
2828 
2829  if( (BCtype[X1DN]==CYLAXIS) ){
2830 
2831 
2832 
2833  /* inner radial BC (preserves u^t rho and u) */
2834  if ( (totalsize[1]>1) && (mycpupos[1] == 0) ) {
2836 
2837  { // start block
2838  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
2840 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2841  OPENMPBCLOOPBLOCK{
2842  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
2843 
2844  rj=j;
2845  if(RADNT_FULLPHI==2.0*Pi){
2846  rk=k+N3/2;
2847  if(rk>=N3) rk-=N3;
2848  }
2849  else rk=k;
2850  LOOPBOUND1IN{
2851  PBOUNDLOOP(pliter,pl){
2852  // SECTIONMARK: assume r=0 singularity can't move
2853  if(dirprim[pl]==FACE1 || dirprim[pl]==CORN3 || dirprim[pl]==CORN2 || dirprim[pl]==CORNT ) ri = -i; // FACE1 values
2854  else ri=-i-1; // "CENT" values for purposes of this BC
2855  MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
2856  }// over pl
2857  }// over boundary zones
2858  }
2859  }// end block
2860 
2861 
2862 
2863  if( (BCtype[X1DN]==CYLAXIS) ){
2864 
2865  /* make sure b and u are antisymmetric at the poles (preserves u^t rho and u) */
2867 
2868  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
2870 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2871  OPENMPBCLOOPBLOCK{
2872  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
2873 
2874 
2875  // SECTIONMARK: assume r=0 singularity can't move
2876  i=0;
2877  PBOUNDLOOP(pliter,pl){
2878  if(pl==U1 || pl==B1){
2879  if(dirprim[pl]==FACE1 || dirprim[pl]==CORN3 || dirprim[pl]==CORN2 || dirprim[pl]==CORNT ){
2880  MACP0A1(prim,i,j,k,pl) = 0.0;
2881  }
2882  }// else don't do this pl
2883  } // end over pl
2884 
2885  LOOPBOUND1IN {
2886  PBOUNDLOOP(pliter,pl){
2887  if(pl==U1 || pl==B1){
2888  MACP0A1(prim,i,j,k,pl) *= -1.;
2889  }// end if right pl
2890  } // end over pl
2891  } // end over boundary zones
2892  }// end loop 23
2893  }
2894  } //end if inner CPU wall
2895  }
2896  else{
2897  dualfprintf(fail_file,"Shouldn't be here in bounds\n");
2898  myexit(3946840);
2899  }
2900 
2901  } // end parallel region
2902 
2903  // if full 2pi, then for MPI would have to deal with separately.
2904 
2905  return(0);
2906 
2907 }
2908 
2909 
2910 
2911 
2912 
2913 
2914 
2915 
2917 int bound_x1up_radcylbeam(
2918  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
2919  int *inboundloop,
2920  int *outboundloop,
2921  int *innormalloop,
2922  int *outnormalloop,
2923  int (*inoutlohi)[NUMUPDOWN][NDIM],
2924  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
2925  int *dosetbc,
2926  int enerregion,
2927  int *localenerpos
2928  )
2929 
2930 {
2931 
2932 
2933 
2934 
2935 
2936 #pragma omp parallel // assume don't require EOS
2937  {
2938 
2939  int i,j,k,pl,pliter;
2940  FTYPE vcon[NDIM],X[NDIM],V[NDIM];
2941 #if(WHICHVEL==VEL3)
2942  int failreturn;
2943 #endif
2944  int ri, rj, rk; // reference i,j,k
2945  FTYPE prescale[NPR];
2946  int jj,kk;
2947  struct of_geom geomdontuse[NPR];
2948  struct of_geom *ptrgeom[NPR];
2949  struct of_geom rgeomdontuse[NPR];
2950  struct of_geom *ptrrgeom[NPR];
2951 
2952  // assign memory
2953  PALLLOOP(pl){
2954  ptrgeom[pl]=&(geomdontuse[pl]);
2955  ptrrgeom[pl]=&(rgeomdontuse[pl]);
2956  }
2957 
2958 
2959  extern FTYPE RADNT_OMSCALE;
2960  extern FTYPE RADNT_FULLPHI;
2961 
2962 
2963  if(BCtype[X1UP]==RADCYLBEAMBC && (totalsize[1]>1) && (mycpupos[1] == ncpux1-1) ){
2964 
2965 
2966  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
2968 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
2969  OPENMPBCLOOPBLOCK{
2970  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
2971 
2972 
2973  ri=riout;
2974  rj=j;
2975  rk=k;
2976 
2977 
2978  // ptrrgeom : i.e. ref geom
2979  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
2980 
2981  FTYPE *pr;
2982  LOOPBOUND1OUT{
2983 
2984  pr = &MACP0A1(prim,i,j,k,0);
2985 
2986 
2987  //initially copying everything
2988  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
2989 
2990  if(ispstag==0){
2991  // local geom
2992  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
2993 
2994  //coordinates of the ghost cell
2995  bl_coord_ijk_2(i,j,k,CENT,X, V);
2996 
2997 
2998  FTYPE xx,yy,zz,rsq;
2999  coord(i, j, k, CENT, X);
3000  bl_coord(X, V);
3001  xx=V[1];
3002  yy=V[2];
3003  zz=V[3];
3004 
3005 
3006  pr[RHO] = 1;
3007  pr[UU] = 0.1;
3008 
3009  //Keplerian gas
3010  FTYPE rCYL=V[1];
3011  FTYPE Om=RADNT_OMSCALE/(a+pow(rCYL,1.5));
3012 
3013  pr[U1] = 0.0;
3014  pr[U2] = 0.0;
3015  pr[U3] = Om;
3016 
3017  //E, F^i in orthonormal fluid frame
3018  FTYPE pradffortho[NPR];
3019  pradffortho[PRAD0] = calc_LTE_EfromT(1.e10/TEMPBAR);
3020  pradffortho[PRAD0] = 1.0;
3021  pradffortho[PRAD1] = 0;
3022  pradffortho[PRAD2] = 0;
3023  pradffortho[PRAD3] = 0;
3024 
3025  int whichvel=VEL3;
3026  int whichcoordfluid=MCOORD;
3027  int whichcoordrad=whichcoordfluid;
3028  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
3029 
3030  }// end if not staggered fields
3031 
3032  }// end loop over outer i's
3033 
3034  }// end over loop
3035  }// end if correct boundary condition and core
3036 
3037 
3038 
3039  }// end parallel region
3040 
3041 
3042 
3043 
3044 
3045 
3046 
3047  return(0);
3048 }
3049 
3051 int bound_x1up_radcyljet(
3052  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
3053  int *inboundloop,
3054  int *outboundloop,
3055  int *innormalloop,
3056  int *outnormalloop,
3057  int (*inoutlohi)[NUMUPDOWN][NDIM],
3058  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
3059  int *dosetbc,
3060  int enerregion,
3061  int *localenerpos
3062  )
3063 
3064 {
3065 
3066 
3067 
3068  extern FTYPE RADCYLJET_TYPE;
3069 
3070 
3071 #pragma omp parallel // assume don't require EOS
3072  {
3073 
3074  int i,j,k,pl,pliter;
3075  FTYPE vcon[NDIM],X[NDIM],V[NDIM];
3076 #if(WHICHVEL==VEL3)
3077  int failreturn;
3078 #endif
3079  int ri, rj, rk; // reference i,j,k
3080  FTYPE prescale[NPR];
3081  int jj,kk;
3082  struct of_geom geomdontuse[NPR];
3083  struct of_geom *ptrgeom[NPR];
3084  struct of_geom rgeomdontuse[NPR];
3085  struct of_geom *ptrrgeom[NPR];
3086 
3087  // assign memory
3088  PALLLOOP(pl){
3089  ptrgeom[pl]=&(geomdontuse[pl]);
3090  ptrrgeom[pl]=&(rgeomdontuse[pl]);
3091  }
3092 
3093 
3094 
3095  if(BCtype[X1UP]==RADCYLJETBC && (totalsize[1]>1) && (mycpupos[1] == ncpux1-1) ){
3096 
3097 
3098  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
3100 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
3101  OPENMPBCLOOPBLOCK{
3102  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
3103 
3104 
3105  ri=riout;
3106  rj=j;
3107  rk=k;
3108 
3109 
3110  // ptrrgeom : i.e. ref geom
3111  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
3112 
3113  FTYPE *pr;
3114  LOOPBOUND1OUT{
3115 
3116  pr = &MACP0A1(prim,i,j,k,0);
3117 
3118 
3119  //initially copying everything
3120  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
3121 
3122  if(ispstag==0){
3123  // local geom
3124  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
3125 
3126  //coordinates of the ghost cell
3127  bl_coord_ijk_2(i,j,k,CENT,X, V);
3128 
3129 
3130  FTYPE xx,yy,zz,rsq;
3131  coord(i, j, k, CENT, X);
3132  bl_coord(X, V);
3133  xx=V[1];
3134  yy=V[2];
3135  zz=V[3];
3136 
3137 
3138  FTYPE pradffortho[NPR];
3139  if(RADCYLJET_TYPE==2){
3140  // pr[RHO] = 1;
3141  // pr[UU] = 0.1;
3142 
3143  pr[U1] = 0.0;
3144  pr[U2] = 0.0;
3145  pr[U3] = 0.0;
3146 
3147  //E, F^i in orthonormal fluid frame
3148  pradffortho[PRAD0] = calc_LTE_EfromT(1.e10/TEMPBAR);
3149  pradffortho[PRAD0] = 1.0; // should be compared to Ehatjet in init.koral.c
3150  pradffortho[PRAD1] = 0;
3151  pradffortho[PRAD2] = 0;
3152  pradffortho[PRAD3] = 0;
3153  }
3154  if(RADCYLJET_TYPE==3){
3155  pr[RHO] = 0.1; // made same as rhojet to keep density flatish
3156  pr[UU] = 0.1; // ""
3157 
3158  pr[U1] = 0.0;
3159  pr[U2] = 0.0;
3160  pr[U3] = 0.0;
3161 
3162  //E, F^i in orthonormal fluid frame
3163  pradffortho[PRAD0] = calc_LTE_EfromT(1.e10/TEMPBAR);
3164  pradffortho[PRAD0] = 1.0; // should be compared to Ehatjet in init.koral.c
3165  pradffortho[PRAD1] = -0.1*pradffortho[PRAD0];
3166  pradffortho[PRAD2] = 0;
3167  pradffortho[PRAD3] = 0;
3168  }
3169 
3170 
3171  int whichvel=VEL3;
3172  int whichcoordfluid=MCOORD;
3173  int whichcoordrad=whichcoordfluid;
3174  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
3175 
3176  }// end if not staggered fields
3177 
3178  }// end loop over outer i's
3179 
3180  }// end over loop
3181  }// end if correct boundary condition and core
3182 
3183 
3184 
3185  }// end parallel region
3186 
3187 
3188 
3189 
3190 
3191 
3192 
3193  return(0);
3194 }
3195 
3196 
3198 int bound_x2dn_radcyljet(
3199  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
3200  int *inboundloop,
3201  int *outboundloop,
3202  int *innormalloop,
3203  int *outnormalloop,
3204  int (*inoutlohi)[NUMUPDOWN][NDIM],
3205  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
3206  int *dosetbc,
3207  int enerregion,
3208  int *localenerpos
3209  )
3210 
3211 {
3212 
3213 
3214 
3215  extern FTYPE RADCYLJET_TYPE;
3216 
3217 
3218 #pragma omp parallel // assume don't require EOS
3219  {
3220 
3221  int i,j,k,pl,pliter;
3222  FTYPE vcon[NDIM],X[NDIM],V[NDIM];
3223 #if(WHICHVEL==VEL3)
3224  int failreturn;
3225 #endif
3226  int ri, rj, rk; // reference i,j,k
3227  FTYPE prescale[NPR];
3228  int jj,kk;
3229  struct of_geom geomdontuse[NPR];
3230  struct of_geom *ptrgeom[NPR];
3231  struct of_geom rgeomdontuse[NPR];
3232  struct of_geom *ptrrgeom[NPR];
3233 
3234  // assign memory
3235  PALLLOOP(pl){
3236  ptrgeom[pl]=&(geomdontuse[pl]);
3237  ptrrgeom[pl]=&(rgeomdontuse[pl]);
3238  }
3239 
3240 
3241 
3242  if(BCtype[X2DN]==RADCYLJETBC && (totalsize[2]>1) && (mycpupos[2] == 0) ){
3243 
3244 
3245  OPENMPBCLOOPVARSDEFINELOOPX2DIR; OPENMPBCLOOPSETUPLOOPX2DIR;
3247 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
3248  OPENMPBCLOOPBLOCK{
3249  OPENMPBCLOOPBLOCK2IJKLOOPX2DIR(i,k);
3250 
3251 
3252  ri=i;
3253  rj=rjin;
3254  rk=k;
3255 
3256 
3257  // ptrrgeom : i.e. ref geom
3258  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
3259 
3260  FTYPE *pr;
3261  LOOPBOUND2IN{
3262 
3263  pr = &MACP0A1(prim,i,j,k,0);
3264 
3265 
3266  //initially copying everything
3267  PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
3268 
3269  if(ispstag==0){
3270  // local geom
3271  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
3272 
3273  extern int jetbound(int i, int j, int k, int loc, FTYPE *prin, FTYPE *prflux, FTYPE (*prim)[NSTORE2][NSTORE3][NPR]);
3274  int insidejet=jetbound(i,j,k,CENT,MAC(prim,i,j,k),MAC(prim,i,j,k),prim);
3275 
3276 
3277  }// end if not staggered fields
3278 
3279  }// end loop over outer i's
3280 
3281  }// end over loop
3282  }// end if correct boundary condition and core
3283 
3284 
3285 
3286  }// end parallel region
3287 
3288 
3289 
3290 
3291 
3292 
3293 
3294  return(0);
3295 }
3296 
3297 
3298 
3299 
3300 
3301 
3302 
3303 
3304 
3305 
3306 
3308 int bound_radbeam2dksvertbeaminflow(int dir,
3309  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
3310  int *inboundloop,
3311  int *outboundloop,
3312  int *innormalloop,
3313  int *outnormalloop,
3314  int (*inoutlohi)[NUMUPDOWN][NDIM],
3315  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
3316  int *dosetbc,
3317  int enerregion,
3318  int *localenerpos
3319  )
3320 
3321 {
3322 
3323 
3324 #pragma omp parallel // assume don't require EOS
3325  {
3326 
3327  int i,j,k,pl,pliter;
3328  FTYPE vcon[NDIM],X[NDIM],V[NDIM];
3329 #if(WHICHVEL==VEL3)
3330  int failreturn;
3331 #endif
3332  int ri, rj, rk; // reference i,j,k
3333  FTYPE prescale[NPR];
3334  int jj,kk;
3335  struct of_geom geomdontuse[NPR];
3336  struct of_geom *ptrgeom[NPR];
3337  struct of_geom rgeomdontuse[NPR];
3338  struct of_geom *ptrrgeom[NPR];
3339 
3340  // assign memory
3341  PALLLOOP(pl){
3342  ptrgeom[pl]=&(geomdontuse[pl]);
3343  ptrrgeom[pl]=&(rgeomdontuse[pl]);
3344  }
3345 
3346  if(dir==X2UP && BCtype[X2UP]==RADBEAM2DKSVERTBEAMINFLOW && totalsize[2]>1 && mycpupos[2] == ncpux2-1 ){
3347 
3348 
3349  extern int RADBEAM2DKSVERT_BEAMNO,RADBEAM2D_FLATBACKGROUND;
3350  FTYPE RHOAMB=1.e0/RHOBAR;
3351  FTYPE TAMB=1e7/TEMPBAR;
3352  FTYPE PAR_D=1./RHOBAR;
3353  FTYPE PAR_E=1e-4/RHOBAR;
3354 
3355  // BEAM PROPERTIES
3356  int IFBEAM=1; // whether to have a beam
3357  FTYPE TLEFT=1e9/TEMPBAR;
3358  // FTYPE NLEFT=0.99;
3359  FTYPE NLEFT=0.995;
3360  // FTYPE NLEFT=0.999;
3361  // FTYPE NLEFT=0.999999; // paper says this, while koral code says 0.999
3362 
3363  FTYPE BEAML,BEAMR;
3364  if (RADBEAM2DKSVERT_BEAMNO==1){
3365  BEAML=2.9;
3366  BEAMR=3.1;
3367  }
3368  else if (RADBEAM2DKSVERT_BEAMNO==2){
3369  BEAML=5.8;
3370  BEAMR=6.2;
3371  }
3372  else if (RADBEAM2DKSVERT_BEAMNO==3){
3373  BEAML=15.5;
3374  BEAMR=16.5;
3375  }
3376  else if (RADBEAM2DKSVERT_BEAMNO==4){
3377  BEAML=37;
3378  BEAMR=43;
3379  }
3380  else if (RADBEAM2DKSVERT_BEAMNO==5){
3381  BEAML=7.;
3382  BEAMR=9.;
3383  }
3384 
3385 
3386 
3387  OPENMPBCLOOPVARSDEFINELOOPX2DIR; OPENMPBCLOOPSETUPLOOPX2DIR;
3389 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
3390  OPENMPBCLOOPBLOCK{
3391  OPENMPBCLOOPBLOCK2IJKLOOPX2DIR(i,k);
3392 
3393  ri=i;
3394  rj=rjout;
3395  rk=k;
3396 
3397 
3398  // ptrrgeom : i.e. ref geom
3399  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
3400 
3401 
3402  LOOPBOUND2OUT{
3403  FTYPE *pr=&MACP0A1(prim,i,j,k,0);
3404 
3405  //initially copying everything
3406  PBOUNDLOOP(pliter,pl) pr[pl] = MACP0A1(prim,ri,rj,rk,pl);
3407 
3408  // NOTEMARK: only really makes sense near the hole if in KSCOORDS
3409  if(pr[U2]<0.0) pr[U2]=0.0; // limit so no arbitrary fluid inflow
3410  if(pr[URAD2]<0.0) pr[URAD2]=0.0; // limit so no arbitrary radiative inflow
3411 
3412 
3413  // only overwrite copy if not inside i<0 since want to keep consistent with outflow BCs used for other \theta at those i
3414  FTYPE ERADAMB;
3415  FTYPE rho,uint,Vr;
3416  if(1 && startpos[1]+i>=0 && ispstag==0){ // only do something special with non-field primitives
3417 
3418  // local geom
3419  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
3420 
3421  //coordinates of the ghost cell
3422  bl_coord_ijk_2(i,j,k,CENT,X, V);
3423 
3424  int whichcoord;
3425  whichcoord=MCOORD;
3426 
3427  // get metric grid geometry for these ICs
3428  int getprim=0;
3429  struct of_geom geomrealdontuse;
3430  struct of_geom *ptrgeomreal=&geomrealdontuse;
3431  gset(getprim,whichcoord,i,j,k,ptrgeomreal);
3432 
3433 
3434  if(RADBEAM2D_FLATBACKGROUND){
3435  Vr=0.0;
3436  rho=RHOAMB;
3437  uint=calc_PEQ_ufromTrho(TAMB,rho);
3438  ERADAMB=calc_LTE_EfromT(TAMB);
3439 
3440  // override so like outflow conditions to avoid shear at boundary
3441  // ERADAMB=pr[URAD0];
3442 
3443  }
3444  else{
3445  //zaczynam jednak od profilu analitycznego:
3446  FTYPE r=V[1];
3447  FTYPE mD=PAR_D/(r*r*sqrt(2./r*(1.-2./r)));
3448  FTYPE mE=PAR_E/(pow(r*r*sqrt(2./r),gamideal)*pow(1.-2./r,(gamideal+1.)/4.));
3449  Vr=sqrt(2./r)*(1.-2./r);
3450 
3451 
3452  FTYPE W=1./sqrt(1.-Vr*Vr*ptrgeomreal->gcov[GIND(1,1)]); // assumes RHO location is good for all these quantities
3453  rho=PAR_D/(r*r*sqrt(2./r));
3454  FTYPE T=TAMB;
3455  // FTYPE ERAD=calc_LTE_EfromT(T);
3456  uint=mE/W;
3457  ERADAMB=calc_LTE_Efromurho(uint,rho);
3458 
3459  // override so like outflow conditions to avoid shear at boundary
3460  //ERADAMB=pr[URAD0];
3461  }
3462 
3463 
3464  // beam parameters
3465  FTYPE ERADINJ=calc_LTE_EfromT(TLEFT);
3466  // original top-hat beam
3467  //FTYPE beamshape=(FTYPE)(V[1]>BEAML && V[1]<BEAMR && IFBEAM);
3468  // gaussian-like beam:
3469  FTYPE powbeam=6.0;
3470  FTYPE beamhalfwidth=0.5*(BEAMR-BEAML);
3471  FTYPE beamcenter=(BEAML+BEAMR)*0.5;
3472  FTYPE beamshape=exp(-pow(V[1]-beamcenter,powbeam)/(2.0*pow(beamhalfwidth,powbeam)))*IFBEAM;
3473 
3474 
3475  if(1){
3476  //E, F^i in orthonormal fluid frame
3477  FTYPE Fx,Fy,Fz;
3478  // default flux
3479  Fx=Fy=Fz=0;
3480  // beam flux
3481  Fy=-NLEFT*ERADINJ;
3482 
3483  // dualfprintf(fail_file,"ERADINJ=%g Fy=%g beamshape=%g\n",ERADINJ, Fy,beamshape);
3484 
3485  FTYPE pradffortho[NPR];
3486  pradffortho[PRAD0] = ERADAMB + ERADINJ*beamshape;
3487  pradffortho[PRAD1] = Fx*beamshape;
3488  pradffortho[PRAD2] = Fy*beamshape;
3489  pradffortho[PRAD3] = Fz*beamshape;
3490 
3491  int whichvel=WHICHVEL; // in which vel U1-U3 set
3492  int whichcoordfluid=PRIMECOORDS; // in which coordinates U1-U3 set
3493  int whichcoordrad=MCOORD; // in which coordinates E,F are orthonormal
3494  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
3495 
3496  // dualfprintf(fail_file,"primrad: %g %g %g %g\n",pr[PRAD0],pr[PRAD1],pr[PRAD2],pr[PRAD3]);
3497 
3498  }
3499  else{
3500  // old harm way
3501 
3502  // set radiation quantities as R^t_\nu in orthonormal fluid frame using whichvel velocity and whichcoord coordinates
3503  int whichvel;
3504  whichvel=VEL4;
3505 
3506 
3507  FTYPE uradx,urady,uradz;
3508  FTYPE uradcon[NDIM],othersrad[NUMOTHERSTATERESULTS];
3509 
3510 
3511  // get coordinate basis in VEL4 format
3512  ucon_calc(&pr[URAD1-U1],ptrgeom[URAD1],uradcon,othersrad);
3513  // get coordinate basis in MCOORD basis
3514  uradx=uradcon[1]*sqrt(fabs(ptrgeom[URAD1]->gcov[GIND(1,1)]))/sqrt(fabs(ptrgeomreal->gcov[GIND(1,1)]));
3515  urady=uradcon[2]*sqrt(fabs(ptrgeom[URAD2]->gcov[GIND(2,2)]))/sqrt(fabs(ptrgeomreal->gcov[GIND(2,2)]));
3516  uradz=uradcon[3]*sqrt(fabs(ptrgeom[URAD3]->gcov[GIND(3,3)]))/sqrt(fabs(ptrgeomreal->gcov[GIND(3,3)]));
3517  if(urady<0.0) urady=0.0; // limit so no arbitrary radiative inflow
3518 
3519 
3520  PBOUNDLOOP(pliter,pl){
3521 
3522  //primitives in whichvel,whichcoord
3523  if(V[1]>BEAML && V[1]<BEAMR && IFBEAM){//beam to be imposed
3524 
3525  // override uradz
3526  urady=-1.0/sqrt(1.0 - NLEFT*NLEFT);
3527  uradx=uradz=0.0;
3528 
3529 
3530  if(pl==URAD0) pr[URAD0] = ERADINJ;
3531  if(pl==URAD1) pr[URAD1] = uradx;
3532  if(pl==URAD2) pr[URAD2] = urady;
3533  if(pl==URAD3) pr[URAD3] = uradz;
3534 
3535  // dualfprintf(fail_file,"GOT BEAM: ijk=%d %d %d : %g %g\n",ERADINJ,urady);
3536  }
3537  else{ //no beam
3538 
3539  // dualfprintf(fail_file,"GOTNOBEAM: ijk=%d %d %d : %g %g\n",ERADAMB,urady);
3540 
3541  // pr[URAD0] = ERADAMB;
3542  if(pl==URAD0) pr[URAD0] = ERADAMB; // so matches outer radial boundary when no beam
3543  if(pl==URAD1) pr[URAD1] = uradx;
3544  if(pl==URAD2) pr[URAD2] = urady;
3545  if(pl==URAD3) pr[URAD3] = uradz;
3546  }
3547  } // over allowed pl's to bound
3548 
3549 
3550  // KORALTODO: ERADINJ is in fluid frame, need to convert, but probably ok.
3551 
3552  // get all primitives in WHICHVEL/PRIMECOORDS value
3553  primefluid_EVrad_to_primeall(&whichvel, &whichcoord, ptrgeom[RHO],pr,pr); // assumes ptrgeom[RHO] is same location as all other primitives (as is currently true).
3554  }
3555 
3556  }// end if not staggered field
3557 
3558 
3559  }// end loop over inner i's
3560  } // over block
3561  }// end if correct BC and should be doind BC for this core
3562  }// end parallel region
3563 
3564  return(0);
3565 }
3566 
3567 
3568 
3569 
3570 
3571 
3572 
3573 
3574 
3575 
3576 
3577 
3578 
3579 
3581 int bound_radcylbeamcart(int dir,
3582  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
3583  int *inboundloop,
3584  int *outboundloop,
3585  int *innormalloop,
3586  int *outnormalloop,
3587  int (*inoutlohi)[NUMUPDOWN][NDIM],
3588  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
3589  int *dosetbc,
3590  int enerregion,
3591  int *localenerpos
3592  )
3593 
3594 {
3595 
3596 
3597 
3598 
3599 
3600 #pragma omp parallel // assume don't require EOS
3601  {
3602 
3603  int i,j,k,pl,pliter;
3604  int ri, rj, rk; // reference i,j,k
3605  struct of_geom rgeomdontuse[NPR];
3606  struct of_geom *ptrrgeom[NPR];
3607  // assign memory
3608  PALLLOOP(pl) ptrrgeom[pl]=&(rgeomdontuse[pl]);
3609 
3610 
3611 
3612  int get_radcylbeamcart(int dir, int *dirprim, int ispstag, int ri, int rj, int rk, int i, int j, int k, FTYPE *prref, FTYPE *pr);
3613 
3614 
3616  // X1DN
3618  if(dir==X1DN && totalsize[1]>1 && mycpupos[1] == 0){
3619 
3620 
3621  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
3623 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
3624  OPENMPBCLOOPBLOCK{
3625  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
3626 
3627 
3628  ri=riin;
3629  rj=j;
3630  rk=k;
3631 
3632  // ptrrgeom : i.e. ref geom
3633  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
3634 
3635  LOOPBOUND1IN{
3636 
3637  FTYPE *pr = &MACP0A1(prim,i,j,k,0);
3638  FTYPE *prref = &MACP0A1(prim,ri,rj,rk,0);
3639  get_radcylbeamcart(dir, dirprim, ispstag,ri,rj,rk,i,j,k,prref,pr);
3640 
3641  }// end loop over outer i's
3642 
3643  }// end over loop
3644  }// end if correct boundary condition and core
3645 
3646 
3647 
3649  // X1UP
3651  if(dir==X1UP && totalsize[1]>1 && mycpupos[1] == ncpux1-1){
3652 
3653 
3654  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
3656 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
3657  OPENMPBCLOOPBLOCK{
3658  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
3659 
3660 
3661  ri=riout;
3662  rj=j;
3663  rk=k;
3664 
3665  // ptrrgeom : i.e. ref geom
3666  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
3667 
3668  LOOPBOUND1OUT{
3669 
3670  FTYPE *pr = &MACP0A1(prim,i,j,k,0);
3671  FTYPE *prref = &MACP0A1(prim,ri,rj,rk,0);
3672  get_radcylbeamcart(dir, dirprim, ispstag,ri,rj,rk,i,j,k,prref,pr);
3673 
3674  }// end loop over outer i's
3675 
3676  }// end over loop
3677  }// end if correct boundary condition and core
3678 
3679 
3680 
3681 
3682 
3684  // X2DN
3686  if(dir==X2DN && totalsize[2]>1 && mycpupos[2] == 0){
3687 
3688 
3689  OPENMPBCLOOPVARSDEFINELOOPX2DIR; OPENMPBCLOOPSETUPLOOPX2DIR;
3691 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
3692  OPENMPBCLOOPBLOCK{
3693  OPENMPBCLOOPBLOCK2IJKLOOPX2DIR(i,k);
3694 
3695 
3696  ri=i;
3697  rj=rjin;
3698  rk=k;
3699 
3700 
3701  // ptrrgeom : i.e. ref geom
3702  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
3703 
3704  LOOPBOUND2IN{
3705 
3706  FTYPE *pr = &MACP0A1(prim,i,j,k,0);
3707  FTYPE *prref = &MACP0A1(prim,ri,rj,rk,0);
3708  get_radcylbeamcart(dir, dirprim, ispstag,ri,rj,rk,i,j,k,prref,pr);
3709 
3710  }// end inner loop
3711 
3712  }// end over loop
3713  }// end if correct boundary condition and core
3714 
3715 
3717  // X2UP
3719  if(dir==X2UP && totalsize[2]>1 && mycpupos[2] == ncpux2-1){
3720 
3721 
3722  OPENMPBCLOOPVARSDEFINELOOPX2DIR; OPENMPBCLOOPSETUPLOOPX2DIR;
3724 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
3725  OPENMPBCLOOPBLOCK{
3726  OPENMPBCLOOPBLOCK2IJKLOOPX2DIR(i,k);
3727 
3728 
3729  ri=i;
3730  rj=rjout;
3731  rk=k;
3732 
3733 
3734  // ptrrgeom : i.e. ref geom
3735  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
3736 
3737  LOOPBOUND2OUT{
3738 
3739  FTYPE *pr = &MACP0A1(prim,i,j,k,0);
3740  FTYPE *prref = &MACP0A1(prim,ri,rj,rk,0);
3741  get_radcylbeamcart(dir, dirprim, ispstag,ri,rj,rk,i,j,k,prref,pr);
3742 
3743  }// end inner loop
3744 
3745  }// end over loop
3746  }// end if correct boundary condition and core
3747 
3749  // X3DN
3751  if(dir==X3DN && totalsize[3]>1 && mycpupos[3] == 0){
3752 
3753 
3754  OPENMPBCLOOPVARSDEFINELOOPX3DIR; OPENMPBCLOOPSETUPLOOPX3DIR;
3756 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
3757  OPENMPBCLOOPBLOCK{
3758  OPENMPBCLOOPBLOCK2IJKLOOPX3DIR(i,j);
3759 
3760 
3761  ri=i;
3762  rj=j;
3763  rk=rkin;
3764 
3765 
3766  // ptrrgeom : i.e. ref geom
3767  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
3768 
3769  LOOPBOUND3IN{
3770 
3771  FTYPE *pr = &MACP0A1(prim,i,j,k,0);
3772  FTYPE *prref = &MACP0A1(prim,ri,rj,rk,0);
3773  get_radcylbeamcart(dir, dirprim, ispstag,ri,rj,rk,i,j,k,prref,pr);
3774 
3775  }// end inner loop
3776 
3777  }// end over loop
3778  }// end if correct boundary condition and core
3779 
3780 
3782  // X3UP
3784  if(dir==X3UP && totalsize[3]>1 && mycpupos[3] == ncpux3-1){
3785 
3786 
3787  OPENMPBCLOOPVARSDEFINELOOPX3DIR; OPENMPBCLOOPSETUPLOOPX3DIR;
3789 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
3790  OPENMPBCLOOPBLOCK{
3791  OPENMPBCLOOPBLOCK2IJKLOOPX3DIR(i,j);
3792 
3793 
3794  ri=i;
3795  rj=j;
3796  rk=rkout;
3797 
3798 
3799  // ptrrgeom : i.e. ref geom
3800  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
3801 
3802  LOOPBOUND3OUT{
3803 
3804  FTYPE *pr = &MACP0A1(prim,i,j,k,0);
3805  FTYPE *prref = &MACP0A1(prim,ri,rj,rk,0);
3806  get_radcylbeamcart(dir, dirprim, ispstag,ri,rj,rk,i,j,k,prref,pr);
3807 
3808  }// end inner loop
3809 
3810  }// end over loop
3811  }// end if correct boundary condition and core
3812 
3813 
3814 
3815 
3816 
3817  }// end parallel region
3818 
3819 
3820 
3821 
3822 
3823 
3824 
3825  return(0);
3826 }
3827 
3828 
3829 
3830 int get_radcylbeamcart(int dir, int *dirprim, int ispstag, int ri, int rj, int rk, int i, int j, int k, FTYPE *prref, FTYPE *pr)
3831 {
3832  int pliter,pl;
3833 
3834 
3835  //initially copying everything
3836  PBOUNDLOOP(pliter,pl) pr[pl] = prref[pl];
3837 
3838 
3839  if(ispstag==1) return(0); // nothing else to do for now
3840 
3841 
3842  int src;
3843  if(dir==X1DN || dir==X1UP) src=(j>0.9*totalsize[2]/2. && j<1.1*totalsize[2]/2.);
3844  else if(dir==X2DN || dir==X2UP) src=(i>0.9*totalsize[1]/2. && i<1.1*totalsize[1]/2.);
3845  else if(dir==X3DN || dir==X3UP){
3846  dualfprintf(fail_file,"Shouldn't be here for X3 dir\n");
3847  myexit(87253245);
3848  }
3849 
3850  // densities
3851  pr[RHO] = 1;
3852  pr[UU] = 0.1;
3853 
3854 
3855  if(src){
3856 
3857  FTYPE X[NDIM],V[NDIM];
3858  int jj,kk;
3859  struct of_geom geomdontuse[NPR];
3860  struct of_geom *ptrgeom[NPR];
3861  // assign memory
3862  PALLLOOP(pl) ptrgeom[pl]=&(geomdontuse[pl]);
3863 
3864  // local geom
3865  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
3866 
3867  //coordinates of the ghost cell
3868  bl_coord_ijk_2(i,j,k,CENT,X, V);
3869  FTYPE xx=V[1];
3870  FTYPE yy=V[2];
3871  FTYPE zz=V[3];
3872 
3873  //Keplerian gas
3874  extern FTYPE RADNT_OMSCALE;
3875  FTYPE Omx=-RADNT_OMSCALE/(a+pow(fabs(yy),1.5))*yy;
3876  FTYPE Omy=RADNT_OMSCALE/(a+pow(fabs(xx),1.5))*xx;
3877 
3878  // dualfprintf(fail_file,"dir=%d RADNT_OMSCALE=%g xx=%g a=%g Omx=%g Omy=%g\n",dir,RADNT_OMSCALE,xx,a,Omx,Omy);
3879 
3880  if(dir==X1DN || dir==X1UP){
3881  pr[U1] = 0.0;
3882  pr[U2] = Omy;
3883  }
3884  else if(dir==X2DN || dir==X2UP){
3885  pr[U1] = Omx;
3886  pr[U2] = 0.0;
3887  }
3888 
3889  pr[U3] = 0.0;
3890 
3891  //E, F^i in orthonormal fluid frame
3892  FTYPE pradffortho[NPR];
3893  pradffortho[PRAD0] = calc_LTE_EfromT(1.e10/TEMPBAR);
3894  pradffortho[PRAD0] = 1.0;
3895  pradffortho[PRAD1] = 0;
3896  pradffortho[PRAD2] = 0;
3897  pradffortho[PRAD3] = 0;
3898 
3899  int whichvel=VEL3;
3900  int whichcoordfluid=MCOORD;
3901  int whichcoordrad=whichcoordfluid;
3902  whichfluid_ffrad_to_primeall(&whichvel, &whichcoordfluid, &whichcoordrad, ptrgeom[RHO], pradffortho, pr, pr);
3903  }
3904  else{
3905  pr[U1] = 0.0;
3906  pr[U2] = 0.0;
3907  pr[U3] = 0.0;
3908 
3909  pr[URAD0] = 0.001;
3910  pr[URAD1] = 0.0;
3911  pr[URAD2] = 0.0;
3912  pr[URAD3] = 0.0;
3913  }
3914 
3915  return(0);
3916 }
3917 
3918 
3919 
3920 
3921 
3922 
3924 int bound_staticset(int dir,
3925  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
3926  int *inboundloop,
3927  int *outboundloop,
3928  int *innormalloop,
3929  int *outnormalloop,
3930  int (*inoutlohi)[NUMUPDOWN][NDIM],
3931  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
3932  int *dosetbc,
3933  int enerregion,
3934  int *localenerpos
3935  )
3936 
3937 {
3938 
3939 
3940 
3941 
3942 
3943 #pragma omp parallel // assume don't require EOS
3944  {
3945 
3946  int i,j,k,pl,pliter;
3947  FTYPE vcon[NDIM],X[NDIM],V[NDIM];
3948 #if(WHICHVEL==VEL3)
3949  int failreturn;
3950 #endif
3951  int ri, rj, rk; // reference i,j,k
3952  FTYPE prescale[NPR];
3953  int jj,kk;
3954  struct of_geom geomdontuse[NPR];
3955  struct of_geom *ptrgeom[NPR];
3956  struct of_geom rgeomdontuse[NPR];
3957  struct of_geom *ptrrgeom[NPR];
3958 
3959  // assign memory
3960  PALLLOOP(pl){
3961  ptrgeom[pl]=&(geomdontuse[pl]);
3962  ptrrgeom[pl]=&(rgeomdontuse[pl]);
3963  }
3964 
3965 
3966  if((BCtype[X1DN]==OUTFLOWSTATIC || BCtype[X1DN]==HORIZONOUTFLOWSTATIC) && (totalsize[1]>1) && (mycpupos[1] == 0) ){
3967 
3968 
3969  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
3971 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
3972  OPENMPBCLOOPBLOCK{
3973  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
3974 
3975 
3976  ri=riin;
3977  rj=j;
3978  rk=k;
3979 
3980 
3981  // ptrrgeom : i.e. ref geom
3982  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
3983 
3984  FTYPE *pr;
3985  LOOPBOUND1IN{
3986 
3987  pr = &MACP0A1(prim,i,j,k,0);
3988 
3989 
3990  //initially copying everything
3991  // PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
3992 
3993  if(ispstag==0){
3994 
3995  //initially copying everything
3996  PBOUNDLOOP(pliter,pl) if(pl==U1) MACP0A1(prim,i,j,k,pl) = 0.0; // static
3997 
3998  }// end if not staggered fields
3999 
4000  }// end loop over outer i's
4001 
4002  }// end over loop
4003  }// end if correct boundary condition and core
4004 
4005 
4006 
4007  if((BCtype[X1UP]==OUTFLOWSTATIC || BCtype[X1UP]==HORIZONOUTFLOWSTATIC) && (totalsize[1]>1) && (mycpupos[1] == ncpux1-1) ){
4008 
4009 
4010  OPENMPBCLOOPVARSDEFINELOOPX1DIR; OPENMPBCLOOPSETUPLOOPX1DIR;
4012 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
4013  OPENMPBCLOOPBLOCK{
4014  OPENMPBCLOOPBLOCK2IJKLOOPX1DIR(j,k);
4015 
4016 
4017  ri=riout;
4018  rj=j;
4019  rk=k;
4020 
4021 
4022  // ptrrgeom : i.e. ref geom
4023  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
4024 
4025  FTYPE *pr;
4026  LOOPBOUND1OUT{
4027 
4028  pr = &MACP0A1(prim,i,j,k,0);
4029 
4030 
4031  //initially copying everything
4032  // PBOUNDLOOP(pliter,pl) MACP0A1(prim,i,j,k,pl) = MACP0A1(prim,ri,rj,rk,pl);
4033 
4034  if(ispstag==0){
4035 
4036  //initially copying everything
4037  PBOUNDLOOP(pliter,pl) if(pl==U1) MACP0A1(prim,i,j,k,pl) = 0.0; // static
4038 
4039  }// end if not staggered fields
4040 
4041  }// end loop over outer i's
4042 
4043  }// end over loop
4044  }// end if correct boundary condition and core
4045 
4046 
4047 
4048 
4049 
4050 
4051  }// end parallel region
4052 
4053 
4054 
4055 
4056 
4057 
4058 
4059  return(0);
4060 }
4061 
4062 
4063 
4064 
4065 
4067 int bound_waldmono(int dir,
4068  int boundstage, int finalstep, SFTYPE boundtime, int whichdir, int boundvartype, int *dirprim, int ispstag, FTYPE (*prim)[NSTORE2][NSTORE3][NPR],
4069  int *inboundloop,
4070  int *outboundloop,
4071  int *innormalloop,
4072  int *outnormalloop,
4073  int (*inoutlohi)[NUMUPDOWN][NDIM],
4074  int riin, int riout, int rjin, int rjout, int rkin, int rkout,
4075  int *dosetbc,
4076  int enerregion,
4077  int *localenerpos
4078  )
4079 
4080 {
4081 
4082 
4083 #pragma omp parallel // assume don't require EOS
4084  {
4085 
4086 
4087  int i,j,k,pl,pliter;
4088  FTYPE vcon[NDIM],X[NDIM],V[NDIM];
4089 #if(WHICHVEL==VEL3)
4090  int failreturn;
4091 #endif
4092  int ri, rj, rk; // reference i,j,k
4093  FTYPE prescale[NPR];
4094  int jj,kk;
4095  struct of_geom geomdontuse[NPR];
4096  struct of_geom *ptrgeom[NPR];
4097  struct of_geom rgeomdontuse[NPR];
4098  struct of_geom *ptrrgeom[NPR];
4099 
4100  // assign memory
4101  PALLLOOP(pl){
4102  ptrgeom[pl]=&(geomdontuse[pl]);
4103  ptrrgeom[pl]=&(rgeomdontuse[pl]);
4104  }
4105 
4106 
4107 
4108 
4109 
4110  if(dir==X2UP && BCtype[X2UP]==WALDMONOBC && (totalsize[2]>1) && (mycpupos[2] == ncpux2-1) ){
4111 
4112  OPENMPBCLOOPVARSDEFINELOOPX2DIR; OPENMPBCLOOPSETUPLOOPX2DIR;
4114 #pragma omp for schedule(OPENMPSCHEDULE(),OPENMPCHUNKSIZE(blocksize))
4115  OPENMPBCLOOPBLOCK{
4116  OPENMPBCLOOPBLOCK2IJKLOOPX2DIR(i,k);
4117 
4118  ri=i;
4119  rj=rjout;
4120  rk=k;
4121 
4122  // ptrrgeom : i.e. ref geom
4123  PALLLOOP(pl) get_geometry(ri, rj, rk, dirprim[pl], ptrrgeom[pl]);
4124 
4125  FTYPE *pr;
4126  LOOPBOUND2OUT{
4127  pr = &MACP0A1(prim,i,j,k,0);
4128 
4129  // ASYMM already done before got here, so only change what's necessary to change
4130 
4131  if(ispstag==0){ // only do something special with non-field primitives
4132 
4133  // local geom
4134  PALLLOOP(pl) get_geometry(i, j, k, dirprim[pl], ptrgeom[pl]);
4135 
4136  //coordinates of the ghost cell
4137  bl_coord_ijk_2(i,j,k,CENT,X, V);
4138  FTYPE r;
4139  r=V[1];
4140 
4141  FTYPE rin,rout;
4142  int conddisk;
4143  rin=15.;
4144  rout=25.;
4145  conddisk=1;//(r<rout); // as in new koral
4146 
4147  //hot boundary
4148  if(conddisk){
4149 
4150  //pr[RHO] and pr[UU] remain same as from ASYMM condition as well as any field
4151  //Keplerian gas with no inflow or outflow
4152  pr[U1]=pr[U2]=0.0; // have to be careful with this for VEL3 (must have rin>>rergo).
4153  pr[U3]=0.0; // current koral value
4154 
4155  } // end if actually doing something to boundary cells in "hot" boundary
4156 
4157  }// end if not staggered field
4158 
4159 
4160  }// end loop over outer j's
4161  }
4162  }
4163 
4164 
4165 
4166 
4167 
4168 
4169 
4170  }// end parallel region
4171 
4172  return(0);
4173 }
4174 
Generated on Fri May 20 2016 15:52:32 for HARM by   doxygen 1.8.3.1