HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
boundsflux.c
Go to the documentation of this file.
1 
2 
3 #include "decs.h"
4 
26 
27 
28 
29 #define EXTRAP 0 //atch
30 
31 
32 
33 
41 
42 
43 
44 
45 int bound_flux_user(int boundstage, int finalstep, SFTYPE boundtime, int boundvartype, FTYPE (*F1)[NSTORE2][NSTORE3][NPR+NSPECIAL], FTYPE (*F2)[NSTORE2][NSTORE3][NPR+NSPECIAL], FTYPE (*F3)[NSTORE2][NSTORE3][NPR+NSPECIAL])
46 {
47  int i,j,k,pl,pliter;
48  struct of_geom geom,rgeom;
49  FTYPE vcon[NDIM]; // coordinate basis vcon
50 #if(WHICHVEL==VEL3)
51  int failreturn;
52 #endif
53  int ri, rj, rk; // reference i,j,k
54  FTYPE prescale[NPR];
55  int enerregion;
56  int *localenerpos;
57  int *doflux;
58 
59 
60 
62  //
63  // set bound loop
64  //
67  // enerregion=ACTIVEREGION; // now replaces TRUEGLOBALENERREGION
68  // localenerpos=enerposreg[enerregion];
69  // doflux=dofluxreg[enerregion];
70 
71 
72 
73  // periodic x1
74  if ( (mycpupos[1] == 0)&&(mycpupos[1] == ncpux1 - 1) ) {
75  if( (BCtype[X1DN]==PERIODIC)&&(BCtype[X1UP]==PERIODIC) ){
76  // just copy from one side to another
77 
78  LOOPX1dir{
79 
80  // copy from upper side to lower boundary zones
81  ri=riout;
82  rj=j;
83  rk=k;
84  LOOPBOUND1IN PBOUNDLOOP(pliter,pl) MACP0A1(F1,i,j,k,pl) = MACP0A1(F1,ri+1+i,rj,rk,pl); // for i=0 -> i=N1
85  if(FLUXB==FLUXCTSTAG){
86  if(N2>1) LOOPBOUND1IN PBOUNDLOOP(pliter,pl) MACP0A1(F2,i,j,k,pl) = MACP0A1(F2,ri+1+i,rj,rk,pl);
87  if(N3>1) LOOPBOUND1IN PBOUNDLOOP(pliter,pl) MACP0A1(F3,i,j,k,pl) = MACP0A1(F3,ri+1+i,rj,rk,pl);
88  }
89 
90  // copy from lower side to upper boundary zones
91  ri=riin;
92  rj=j;
93  rk=k;
94  LOOPBOUND1OUT PBOUNDLOOP(pliter,pl) MACP0A1(F1,i,j,k,pl) = MACP0A1(F1,ri+(i-N1),rj,rk,pl); // for i=N1 -> i=0
95  if(FLUXB==FLUXCTSTAG){
96  if(N2>1) LOOPBOUND1OUT PBOUNDLOOP(pliter,pl) MACP0A1(F2,i,j,k,pl) = MACP0A1(F2,ri+(i-N1),rj,rk,pl); // for i=N1 -> i=0
97  if(N3>1) LOOPBOUND1OUT PBOUNDLOOP(pliter,pl) MACP0A1(F3,i,j,k,pl) = MACP0A1(F3,ri+(i-N1),rj,rk,pl); // for i=N1 -> i=0
98  }
99  }
100  }
101  }
102 
103 
104  // periodic x2
105  if ( (mycpupos[2] == 0)&&(mycpupos[2] == ncpux2 - 1) ) {
106  if( (BCtype[X2DN]==PERIODIC)&&(BCtype[X2UP]==PERIODIC) ){
107  // just copy from one side to another
108 
109  LOOPX2dir{
110 
111  // copy from upper side to lower boundary zones
112  ri=i;
113  rj=rjout;
114  rk=k;
115  LOOPBOUND2IN PBOUNDLOOP(pliter,pl) MACP0A1(F2,i,j,k,pl) = MACP0A1(F2,ri,rj+1+j,rk,pl); // for j=0 -> j=N2
116  if(FLUXB==FLUXCTSTAG){
117  if(N1>1) LOOPBOUND2IN PBOUNDLOOP(pliter,pl) MACP0A1(F1,i,j,k,pl) = MACP0A1(F1,ri,rj+1+j,rk,pl); // for j=0 -> j=N2
118  if(N3>1) LOOPBOUND2IN PBOUNDLOOP(pliter,pl) MACP0A1(F3,i,j,k,pl) = MACP0A1(F3,ri,rj+1+j,rk,pl); // for j=0 -> j=N2
119  }
120 
121  // copy from lower side to upper boundary zones
122  ri=i;
123  rj=rjin;
124  rk=k;
125  LOOPBOUND2OUT PBOUNDLOOP(pliter,pl) MACP0A1(F2,i,j,k,pl) = MACP0A1(F2,ri,rj+(j-N2),rk,pl); // for j=N2 -> j=0
126  if(FLUXB==FLUXCTSTAG){
127  if(N1>1) LOOPBOUND2OUT PBOUNDLOOP(pliter,pl) MACP0A1(F1,i,j,k,pl) = MACP0A1(F1,ri,rj+(j-N2),rk,pl); // for j=N2 -> j=0
128  if(N3>1) LOOPBOUND2OUT PBOUNDLOOP(pliter,pl) MACP0A1(F3,i,j,k,pl) = MACP0A1(F3,ri,rj+(j-N2),rk,pl); // for j=N2 -> j=0
129  }
130  }
131  }
132  }
133 
134 
135  // periodic x3
136  if ( (mycpupos[3] == 0)&&(mycpupos[3] == ncpux3 - 1) ) {
137  if( (BCtype[X3DN]==PERIODIC)&&(BCtype[X3UP]==PERIODIC) ){
138  // just copy from one side to another
139 
140  LOOPF_12{
141 
142  // copy from upper side to lower boundary zones
143  ri=i;
144  rj=j;
145  rk=rkout;
146  LOOPBOUND3IN PBOUNDLOOP(pliter,pl) MACP0A1(F3,i,j,k,pl) = MACP0A1(F3,ri,rj,rk+1+k,pl); // for k=0 -> k=N3
147  if(FLUXB==FLUXCTSTAG){
148  if(N1>1) LOOPBOUND3IN PBOUNDLOOP(pliter,pl) MACP0A1(F1,i,j,k,pl) = MACP0A1(F1,ri,rj,rk+1+k,pl); // for k=0 -> k=N3
149  if(N2>1) LOOPBOUND3IN PBOUNDLOOP(pliter,pl) MACP0A1(F2,i,j,k,pl) = MACP0A1(F2,ri,rj,rk+1+k,pl); // for k=0 -> k=N3
150  }
151 
152  // copy from lower side to upper boundary zones
153  ri=i;
154  rj=j;
155  rk=rkin;
156  LOOPBOUND3OUT PBOUNDLOOP(pliter,pl) MACP0A1(F3,i,j,k,pl) = MACP0A1(F3,ri,rj,rk+(k-N3),pl); // for k=N3 -> k=0
157  if(FLUXB==FLUXCTSTAG){
158  if(N1>1) LOOPBOUND3OUT PBOUNDLOOP(pliter,pl) MACP0A1(F1,i,j,k,pl) = MACP0A1(F1,ri,rj,rk+(k-N3),pl); // for k=N3 -> k=0
159  if(N2>1) LOOPBOUND3OUT PBOUNDLOOP(pliter,pl) MACP0A1(F2,i,j,k,pl) = MACP0A1(F2,ri,rj,rk+(k-N3),pl); // for k=N3 -> k=0
160  }
161  }
162  }
163  }
164 
165 
166 
167 
168 
169 
171  //
172  // X1 inner OUTFLOW/FIXEDOUTFLOW
173  //
175 
176  // first allow extrapolation
177  //
178  if (mycpupos[1] == 0) {
179  if(((BCtype[X1DN]==OUTFLOW || BCtype[X1DN]==HORIZONOUTFLOW))||(BCtype[X1DN]==FIXEDOUTFLOW)||(BCtype[X1DN]==FIXED)){ //SASMARK: FIXED is not supposed to be here but need to assign sth to fluxes
180  /* inner r boundary condition: u, just copy */
181  LOOPX1dir{
182  ri=riin;
183  rj=j;
184  rk=k;
185  LOOPBOUND1IN PBOUNDLOOP(pliter,pl){
186 #if(EXTRAP==0)
187  // zeroth orde extrap
188  MACP0A1(F1,i,j,k,pl) = MACP0A1(F1,ri,rj,rk,pl);
189 #else
190  // linear extrap
191  MACP0A1(F1,i,j,k,pl) = MACP0A1(F1,ri,rj,rk,pl)+(MACP0A1(F1,ri+1,rj,rk,pl)-MACP0A1(F1,ri,rj,rk,pl))*(FTYPE)(i-ri);
192 #endif
193  // if OUTFLOW, then disallow extrapolation if flux is from ghost to active zones
194  if((pl==RHO)&&((BCtype[X1DN]==OUTFLOW || BCtype[X1DN]==HORIZONOUTFLOW))){ // only for RHO is it obvious what to do without primitives
195  if(MACP0A1(F1,i,j,k,pl)>0.0) MACP0A1(F1,i,j,k,pl)=0.0; // GODMARK: hope this is enough to shut-down inflow (what about energy flux?)
196  }
197  }// end pl and face loop
198  }// end 2 3
199  }// end if correct bound type
200  else if(BCtype[X1DN]== ASYMM) { //atch: added this; however unsure if any symmetry exists for fluxes
201  LOOPX1dir {
202  LOOPBOUND1IN PBOUNDLOOP(pliter,pl) {
203  MACP0A1(F1,i,j,k,pl) = - MACP0A1(F1,-i,j,k,pl);
204  if( U1 == pl || B1 == pl ) MACP0A1(F1,i,j,k,pl) *= -1.0; //resymm them
205  }
206  }
207  }
208  else if(BCtype[X1DN]== BCEXTRAP_VEL3 || BCtype[X1DN] == BCEXTRAP) {
209  //extrapolates fluxes with 6th order
210  LOOPX1dir{
211  LOOPBOUND1IN PBOUNDLOOP(pliter,pl) {
212  // MACP0A1(F1,i,j,k,pl) = interpn( 6, i,
213  // 0, MACP0A1(F1,0,j,k,pl),
214  // 1, MACP0A1(F1,1,j,k,pl),
215  // 2, MACP0A1(F1,2,j,k,pl),
216  // 3, MACP0A1(F1,3,j,k,pl),
217  // 4, MACP0A1(F1,4,j,k,pl),
218  // 5, MACP0A1(F1,5,j,k,pl) );
219  }
220  }
221  }
222  }// end if mycpupos[1]==0
223 
224 
226  //
227  // X1 outer OUTFLOW/FIXEDOUTFLOW
228  //
230 
231 
232  // outer r BC:
233  if (mycpupos[1] == ncpux1 - 1) {
234  if(((BCtype[X1UP]==OUTFLOW || BCtype[X1UP]==HORIZONOUTFLOW))||(BCtype[X1UP]==FIXEDOUTFLOW)||(BCtype[X1UP]==FIXED)){ //SASMARK: FIXED is not supposed to be here but need to assign sth to fluxes
235  /* outer r BC: outflow */
236 
237  LOOPX1dir{
238  ri=riout;
239  rj=j;
240  rk=k;
241  LOOPBOUND1OUT PBOUNDLOOP(pliter,pl){
242 #if(EXTRAP==0)
243  // zeroth orde extrap
244  MACP0A1(F1,i,j,k,pl) = MACP0A1(F1,ri,rj,rk,pl);
245 #else
246  // linear extrap
247  MACP0A1(F1,i,j,k,pl) = MACP0A1(F1,ri,rj,rk,pl)+(MACP0A1(F1,ri,rj,rk,pl)-MACP0A1(F1,ri-1,rj,rk,pl))*(FTYPE)(i-ri);
248 #endif
249  // if OUTFLOW, then disallow extrapolation if flux is from ghost to active zones
250  if((pl==RHO)&&(BCtype[X1UP]==FIXEDOUTFLOW)){ //SASMARK: changed OUTFLOW to FIXEDOUTFLOW here
251  if(MACP0A1(F1,i,j,k,pl)<0.0) MACP0A1(F1,i,j,k,pl)=0.0; //SASMARK: can be a problem for e.g. Noh problem where there is OUTFLOW BC and the matter actually inflows
252  }
253  }// end pl and face loop
254  }// end 2 3
255  }// end if correct bound type
256  else if(BCtype[X1UP]== ASYMM) { //atch: added this; however unsure if any symmetry exists for fluxes
257  LOOPX1dir {
258  LOOPBOUND1OUT PBOUNDLOOP(pliter,pl) {
259  MACP0A1(F1,i,j,k,pl) = - MACP0A1(F1,N1 + N1 - i,j,k,pl);
260  if( U1 == pl || B1 == pl ) MACP0A1(F1,i,j,k,pl) *= -1.0; //resymm them
261  }
262  }
263  }
264  else if(BCtype[X1UP]== BCEXTRAP_VEL3 || BCtype[X1UP] == BCEXTRAP) {
265  //extrapolates fluxes with 6th order
266  LOOPX1dir{
267  LOOPBOUND1OUT PBOUNDLOOP(pliter,pl) {
268  // MACP0A1(F1,i,j,k,pl) = interpn( 6, i,
269  // N1, MACP0A1(F1,N1,j,k,pl),
270  // N1-1, MACP0A1(F1,N1-1,j,k,pl),
271  // N1-2, MACP0A1(F1,N1-2,j,k,pl),
272  // N1-3, MACP0A1(F1,N1-3,j,k,pl),
273  // N1-4, MACP0A1(F1,N1-4,j,k,pl),
274  // N1-5, MACP0A1(F1,N1-5,j,k,pl) );
275  }
276  }
277  }
278  }// end if mycpu is correct
279 
280 
282  //
283  // X2 inner POLARAXIS
284  //
286 
287 
288  if (mycpupos[2] == 0) {
289  if(((BCtype[X2DN]==OUTFLOW || BCtype[X2DN]==HORIZONOUTFLOW))||(BCtype[X2DN]==FIXEDOUTFLOW)||(BCtype[X2DN]==FIXED)){ //SASMARK: FIXED is not supposed to be here but need to assign sth to fluxes
290  /* inner 2 BC: outflow */
291 
292  LOOPX2dir{
293  ri=i; // correct for flux
294  rj=rjin;
295  rk=k;
296  LOOPBOUND2IN PBOUNDLOOP(pliter,pl){
297 #if(EXTRAP==0)
298  // zeroth orde extrap
299  MACP0A1(F2,i,j,k,pl) = MACP0A1(F2,ri,rj,rk,pl);
300 #else
301  // linear extrap
302  MACP0A1(F2,i,j,k,pl) = MACP0A1(F2,ri,rj,rk,pl)+(MACP0A1(F2,ri,rj,rk,pl)-MACP0A1(F2,ri,rj-l,rk,pl))*(FTYPE)(j-rj);
303 #endif
304  // if OUTFLOW, then disallow extrapolation if flux is from ghost to active zones
305  if((pl==RHO)&&(BCtype[X2DN]==FIXEDOUTFLOW)){ //SASMARK: changed OUTFLOW to FIXEDOUTFLOW here
306  if(MACP0A1(F2,i,j,k,pl)<0.0) MACP0A1(F2,i,j,k,pl)=0.0; //SASMARK: can be a problem for e.g. Noh problem where there is OUTFLOW BC and the matter actually inflows
307  }
308  }// end pl and face loop
309  }// end 2 3
310  }// end if correct bound type
311  else if((BCtype[X2DN]==POLARAXIS)||(BCtype[X2DN]==SYMM)||(BCtype[X2DN]==ASYMM) ){
312  LOOPX2dir{
313  ri=i;
314  rj=rjin;
315  rk=k;
316  LOOPBOUND2IN PBOUNDLOOP(pliter,pl) MACP0A1(F2,i,j,k,pl) = MACP0A1(F2,ri,rj+(rj-j),rk,pl); // self-assigns for j=0
317  }
318  }
319  // F2 of U2/B2 is symmetric, so no need to antisymmetrize
320  // F2 antisymmetric with respect to all other quantities
321 
322  if((BCtype[X2DN]==POLARAXIS)||(BCtype[X2DN]==ASYMM) ){
323 
324  /* make sure b and u are antisymmetric at the poles (preserves u^t rho and u) */
325  LOOPX2dir{
326  LOOPBOUND2IN {
327  //if(POSDEFMETRIC==0){ //SASMARK need to asymmetrize the fluxes; don't understand why should not asymm. the fluxes if POSDEFMETRIC is 0
328  // // u^t must be symmetric across pole, which is functions of u2 and u3 as well as their squares and othe products. u2 in KS happens to be independent of sign, but in general is could be for some other metric.
329  // // for now, assume KS-like metric where u2 is antisymmetric and u^t dep only on u2^2, not u2
330  //}
331  //else{
332  PLOOP(pliter,pl) MACP0A1(F2,i,j,k,pl)*=-1; // anti-sym all
333  MACP0A1(F2,i,j,k,U2)*=-1; // re-sym U2
334  MACP0A1(F2,i,j,k,B2)*=-1; // re-sym B2
335  //}
336  }
337  }// end loop 13
338  } // end if POLARXIS or ASYMM
339  }// end if mycpupos[2]==0
340 
341 
343  //
344  // X2 outer POLARAXIS
345  //
347 
348 
349  if (mycpupos[2] == ncpux2-1) {
350  if(((BCtype[X2UP]==OUTFLOW || BCtype[X2UP]==HORIZONOUTFLOW))||(BCtype[X2UP]==FIXEDOUTFLOW)||(BCtype[X2UP]==FIXED)){ //SASMARK: FIXED is not supposed to be here but need to assign sth to fluxes
351  /* outer 2 BC: outflow */
352 
353  LOOPX2dir{
354  ri=i; // correct for flux
355  rj=rjout;
356  rk=k;
357  LOOPBOUND2OUT PBOUNDLOOP(pliter,pl){
358 #if(EXTRAP==0)
359  // zeroth orde extrap
360  MACP0A1(F2,i,j,k,pl) = MACP0A1(F2,ri,rj,rk,pl);
361 #else
362  // linear extrap
363  MACP0A1(F2,i,j,k,pl) = MACP0A1(F2,ri,rj,rk,pl)+(MACP0A1(F2,ri,rj,rk,pl)-MACP0A1(F2,ri,rj-l,rk,pl))*(FTYPE)(j-rj);
364 #endif
365  // if OUTFLOW, then disallow extrapolation if flux is from ghost to active zones
366  if((pl==RHO)&&(BCtype[X2UP]==FIXEDOUTFLOW)){ //SASMARK: changed OUTFLOW to FIXEDOUTFLOW here
367  if(MACP0A1(F2,i,j,k,pl)>0.0) MACP0A1(F2,i,j,k,pl)=0.0; //SASMARK: can be a problem for e.g. Noh problem where there is OUTFLOW BC and the matter actually inflows
368  }
369  }// end pl and face loop
370  }// end 2 3
371  }// end if correct bound type
372  else if((BCtype[X2UP]==POLARAXIS)||(BCtype[X2UP]==SYMM)||(BCtype[X2UP]==ASYMM) ){
373  LOOPX2dir{
374  ri=i;
375  rj=rjout;
376  rk=k;
377  LOOPBOUND2OUT PBOUNDLOOP(pliter,pl) MACP0A1(F2,i,j,k,pl) = MACP0A1(F2,ri,rj+(rj-j+2),rk,pl); // self-assigns for j=N
378  }
379  }
380 
381  if((BCtype[X2UP]==POLARAXIS)||(BCtype[X2UP]==ASYMM) ){
382 
383  /* make sure b and u are antisymmetric at the poles (preserves u^t rho and u) */
384  LOOPX2dir{
385  LOOPBOUND2OUT {
386  //if(POSDEFMETRIC==0){ //SASMARK need to asymmetrize the fluxes; don't understand why should not asymm. the fluxes if POSDEFMETRIC is 0
387  // // u^t must be symmetric across pole, which is functions of u2 and u3 as well as their squares and othe products. u2 in KS happens to be independent of sign, but in general is could be for some other metric.
388  // // for now, assume KS-like metric where u2 is antisymmetric and u^t dep only on u2^2, not u2
389  //}
390  //else{
391  PLOOP(pliter,pl) MACP0A1(F2,i,j,k,pl)*=-1; // anti-sym all
392  MACP0A1(F2,i,j,k,U2)*=-1; // re-sym U2
393  MACP0A1(F2,i,j,k,B2)*=-1; // re-sym B2
394  //}
395  }
396  }// end loop 13
397  } // end if POLARXIS or ASYMM
398  }// end if mycpupos[2]==ncpux2-1
399 
400 
401 
402  //x3 inner
403  if ( mycpupos[3] == 0 ) {
404  if(((BCtype[X3DN]==OUTFLOW || BCtype[X3DN]==HORIZONOUTFLOW))||(BCtype[X3DN]==FIXEDOUTFLOW)||(BCtype[X3DN]==FIXED)){ //SASMARK: FIXED is not supposed to be here but need to assign sth to fluxes
405  /* inner 3 BC: outflow */
406 
407  LOOPX3dir{
408  ri=i; // correct for flux
409  rj=j;
410  rk=rkin;
411  LOOPBOUND3IN PBOUNDLOOP(pliter,pl){
412 #if(EXTRAP==0)
413  // zeroth orde extrap
414  MACP0A1(F3,i,j,k,pl) = MACP0A1(F3,ri,rj,rk,pl);
415 #else
416  // linear extrap
417  MACP0A1(F3,i,j,k,pl) = MACP0A1(F3,ri,rj,rk,pl)+(MACP0A1(F3,ri,rj,rk,pl)-MACP0A1(F3,ri,rj,rk-l,pl))*(FTYPE)(k-rk);
418 #endif
419  // if OUTFLOW, then disallow extrapolation if flux is from ghost to active zones
420  if((pl==RHO)&&(BCtype[X3DN]==FIXEDOUTFLOW)){ //SASMARK: changed OUTFLOW to FIXEDOUTFLOW here
421  if(MACP0A1(F3,i,j,k,pl)<0.0) MACP0A1(F3,i,j,k,pl)=0.0; //SASMARK: can be a problem for e.g. Noh problem where there is OUTFLOW BC and the matter actually inflows
422  }
423  }// end pl and face loop
424  }// end 2 3
425  }// end if correct bound type
426  }
427 
428  //x3 outer
429  if ( mycpupos[3] == ncpux3 - 1 ) {
430  if(((BCtype[X3UP]==OUTFLOW || BCtype[X3UP]==HORIZONOUTFLOW))||(BCtype[X3UP]==FIXEDOUTFLOW)||(BCtype[X3UP]==FIXED)){ //SASMARK: FIXED is not supposed to be here but need to assign sth to fluxes
431  /* outer 3 BC: outflow */
432 
433  LOOPX3dir{
434  ri=i; // correct for flux
435  rj=j;
436  rk=rkout;
437  LOOPBOUND3OUT PBOUNDLOOP(pliter,pl){
438 #if(EXTRAP==0)
439  // zeroth orde extrap
440  MACP0A1(F3,i,j,k,pl) = MACP0A1(F3,ri,rj,rk,pl);
441 #else
442  // linear extrap
443  MACP0A1(F3,i,j,k,pl) = MACP0A1(F3,ri,rj,rk,pl)+(MACP0A1(F3,ri,rj,rk,pl)-MACP0A1(F3,ri,rj,rk-l,pl))*(FTYPE)(k-rk);
444 #endif
445  // if OUTFLOW, then disallow extrapolation if flux is from ghost to active zones
446  if((pl==RHO)&&(BCtype[X3UP]==FIXEDOUTFLOW)){ //SASMARK: changed OUTFLOW to FIXEDOUTFLOW here
447  if(MACP0A1(F3,i,j,k,pl)>0.0) MACP0A1(F3,i,j,k,pl)=0.0; //SASMARK: can be a problem for e.g. Noh problem where there is OUTFLOW BC and the matter actually inflows
448  }
449  }// end pl and face loop
450  }// end 2 3
451  }// end if correct bound type
452  }
453 
454  return (0);
455 }