HARM
harm and utilities
 All Data Structures Files Functions Variables Typedefs Macros Pages
jon_interp_coorduser.c
Go to the documentation of this file.
1 #define JET6LIKEUSERCOORD 0
2 #define UNIHALFUSERCOORD 1
3 
4 #define WHICHUSERCOORD JET6LIKEUSERCOORD
5 
6 
7 // for defcoord=JET6COORDS like USERCOORDS
8 static FTYPE npow,r1jet,njet1,njet,r0jet,rsjet,Qjet, ntheta,htheta,rsjet2,r0jet2,rsjet3,r0jet3, rs, r0,npow2,cpow2,rbr,x1br, h0,cpow3;
9 
10 
11 void set_coord_parms_nodeps_user(int defcoordlocal)
12 {
13  if(1){
14  // see jet3coords_checknew.nb
15  npow=1.0;
16 
18  // RADIAL GRID SETUP
20  npow=1.0; //don't change it, essentially equivalent to changing cpow2
21 
22  //radial hyperexponential grid
23 
24  //power exponent
25  npow2=6.0; // WALD: 6.0->4.0
26 
27  cpow2=1.0; //exponent prefactor (the larger it is, the more hyperexponentiation is)
28  // cpow3=0.01;
29  cpow3=1.0;
30  //radius at which hyperexponentiation kicks in
31  // rbr = 1E3;
32  rbr = 5E2; // WALD 5E2->5E7
33 
34 
35 
36  // must be same as in dxdxp()
37  // GODMARK: Note njet here is overwritten by njet later, but could have been different values if setup variable names differently.
38  if(0){ // first attempt
39  r1jet=2.8;
40  njet=0.3;
41  r0jet=7.0;
42  rsjet=21.0;
43  Qjet=1.7;
44  }
45  else if(0){ // chosen to resolve disk then resolve jet
46  r1jet=2.8;
47  njet=0.3;
48  r0jet=20.0;
49  rsjet=80.0;
50  Qjet=1.8;
51  }
52  else if(1){
53  r1jet=2.8;
54  njet=0.3;
55  r0jet=15.0;
56  rsjet=40.0;
57  Qjet=1.3; // chosen to help keep jet resolved even within disk region
58  }
59 
60  // for switches from normal theta to ramesh theta
61  rs=40.0; // shift
62  r0=20.0; // divisor
63 
64  // for theta1
65  // hslope=0.3 ; // resolve inner-radial region near equator
66  r0jet3=20.0; // divisor
67  rsjet3=0.0; // subtractor
68 
69  // for theta2
70  h0=0.3; // inner-radial "hslope" for theta2
71  //h0=0.1; // inner-radial "hslope" for theta2 // for thinner disks, change this.
72  // GODMARK: Note that this overwrites above njet!
73  // power \theta_j \propto r^{-njet}
74  njet=1.0;
75 
76 
77  // see fix_3dpoledtissue.nb
78 #if(0)
79  ntheta=21.0;
80  htheta=0.15;
81  rsjet2=5.0;
82  r0jet2=2.0;
83 #else
84  ntheta=5.0;
85  htheta=0.15;
86  rsjet2=5.0;
87  r0jet2=2.0;
88 #endif
89 
90  }
91 
92 }
93 
94 
95 void set_coord_parms_deps_user(int defcoordlocal)
96 {
97  if(1){
98  x1br = log( rbr - R0 ) / npow; //the corresponding X[1] value
99  }
100 
101 }
102 
103 void write_coord_parms_user(int defcoordlocal, FILE *out)
104 {
105  if(1){
106  fprintf(out,"%21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g %21.15g\n",npow,r1jet,njet,r0jet,rsjet,Qjet,ntheta,htheta,rsjet2,r0jet2,rsjet3,r0jet3,rs,r0,npow2,cpow2,rbr,x1br,cpow3);
107  }
108 
109 }
110 void read_coord_parms_user(int defcoordlocal, FILE *in)
111 {
112 
113  if(1){
114  fscanf(in,HEADER19IN,&npow,&r1jet,&njet,&r0jet,&rsjet,&Qjet,&ntheta,&htheta,&rsjet2,&r0jet2,&rsjet3,&r0jet3,&rs,&r0,&npow2,&cpow2,&rbr,&x1br,&cpow3);
115 
116  }
117 }
118 
119 void read_coord_parms_mpi_user(int defcoordlocal)
120 {
121  if(1){
122 #if(USEMPI)
123  MPI_Bcast(&npow, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
124  MPI_Bcast(&r1jet, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
125  MPI_Bcast(&njet, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
126  MPI_Bcast(&r0jet, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
127  MPI_Bcast(&rsjet, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
128  MPI_Bcast(&Qjet, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
129  MPI_Bcast(&ntheta, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
130  MPI_Bcast(&htheta, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
131  MPI_Bcast(&rsjet2, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
132  MPI_Bcast(&r0jet2, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
133  MPI_Bcast(&rsjet3, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
134  MPI_Bcast(&r0jet3, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
135  MPI_Bcast(&rs, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
136  MPI_Bcast(&r0, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
137  MPI_Bcast(&npow2, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
138  MPI_Bcast(&cpow2, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
139  MPI_Bcast(&rbr, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
140  MPI_Bcast(&x1br, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
141  MPI_Bcast(&cpow3, 1, MPI_FTYPE, MPIid[0], MPI_COMM_GRMHD);
142 #endif
143  }
144 
145 }
146 
147 
148 void blcoord_user(FTYPE *X, FTYPE *V)
149 {
150  extern FTYPE mysin(FTYPE th);
151 
152  if(1){
153 
154 #if(0) // no change in exponentiation
155  // JET3COORDS-like radial grid
156  V[1] = R0+exp(pow(X[1],npow)) ;
157 #elif(WHICHUSERCOORD==UNIHALFUSERCOORD)
158 
159  Rout=2000.0;
160  theexp = npow*X[1];
161  npow=1.0;
162  FTYPE gconst1=1.0;
163  FTYPE gconst2=gconst1*.000001;
164  V[1] = R0 + gconst1*X[1] + gconst2*exp(theexp);
165 
166 
167 #elif(WHICHUSERCOORD==JET6LIKEUSERCOORD)
168 
169  FTYPE theexp = npow*X[1];
170  if( X[1] > x1br ) {
171  theexp += cpow2 * pow(X[1]-x1br,npow2);
172  }
173  V[1] = R0+exp(theexp);
174 
175 
176  // FTYPE npowtrue,npowlarger=10.0;
177  // FTYPE npowrs=1E3;
178  // FTYPE npowr0=2E2;
179  // npowtrue = npow + (npowlarger-npow)*(0.5+1.0/M_PI*atan((V[1]-npowrs)/npowr0));
180  // V[1] = R0+exp(pow(X[1],npowtrue)) ;
181 #elif(0)
182  // avoid jump in grid at rbr
183  // determine switches
184  FTYPE r0rbr=rbr/2.0;
185  FTYPE switch0 = 0.5+1.0/M_PI*atan((V[1]-rbr)/r0rbr); // 1 for outer r
186 
187  FTYPE V1 = R0+exp(npow*X[1]);
188  FTYPE V2 = R0+exp(npow*X[1] + cpow2 * pow(cpow3*(X[1]-x1br*1.0),npow2));
189 
190  V[1] = V1*(1.0-switch0) + V2*switch0;
191 
192 #endif
193 
194 
195 
196  FTYPE theta1,theta2,arctan2;
197 
198 
199 #if(0)
200  // JET3COORDS-based:
201  FTYPE myhslope=2.0-Qjet*pow(V[1]/r1jet,-njet*(0.5+1.0/M_PI*atan(V[1]/r0jet-rsjet/r0jet)));
202  theta1 = M_PI * X[2] + ((1. - myhslope) * 0.5) * mysin(2. * M_PI * X[2]);
203 #else
204  // RAMESH BASED
205  // myhslope here is h2 in MCAF paper
206  // // h0 here is h3 in MCAF paper
207  //FTYPE njetvsr;
208  //if(V[1]<rbr) njetvsr=njet;
209  // else njetvsr=njet/(V[1])*rbr;
210  //else njetvsr=
211  //njetvsr=njet;
212 
213  FTYPE localrbr=Rout; //500.0; // rbr;
214  // FTYPE localrbr=rbr;
215  FTYPE localrbrr0=100.0; //MAX(100.0,localrbr/2.0);
216 
217  FTYPE switch0 = 0.5+1.0/M_PI*atan((V[1]-localrbr)/localrbrr0); // large r
218  FTYPE switch2 = 1.0-switch0; // small r
219 
220  // switch0=0.0; switch2=1.0;
221 
222  FTYPE myhslope1=h0 + pow( (V[1]-rsjet3)/r0jet3 , njet);
223  FTYPE myhslope2=h0 + pow( (localrbr-rsjet3)/r0jet3 , njet);
224  FTYPE myhslope = myhslope1*switch2 + myhslope2*switch0;
225 
226  // myhslope = pow(pow(myhslope1,-2.0) + pow(myhslope2,-2.0),-0.5);
227 
228  myhslope=myhslope1;
229 
230  // determine theta2
231  FTYPE myx2;
232  if(X[2]>1.0) myx2=2.0-X[2];
233  else if(X[2]<0.0) myx2=-X[2];
234  else myx2=X[2];
235 
236  FTYPE th2 = 0.5*M_PI*(1.0 + atan(myhslope*(myx2-0.5))/atan(myhslope*0.5));
237 
238  if(X[2]>1.0) th2=2.0*M_PI-th2;
239  else if(X[2]<0.0) th2=-th2;
240 
241  // determine theta0
242  // JET3COORDS-based:
243  myhslope1=2.0-Qjet*pow(V[1]/r1jet,-njet*(0.5+1.0/M_PI*atan(V[1]/r0jet-rsjet/r0jet)));
244  myhslope2=2.0-Qjet*pow(localrbr/r1jet,-njet*(0.5+1.0/M_PI*atan(localrbr/r0jet-rsjet/r0jet)));
245  myhslope = myhslope1*switch2 + myhslope2*switch0;
246  // myhslope here is h0 in MCAF paper
247  FTYPE th0 = M_PI * X[2] + ((1. - myhslope) * 0.5) * mysin(2. * M_PI * X[2]);
248 
249 
250  // determine switches (only function of radius and not x2 or theta)
251  switch0 = 0.5+1.0/M_PI*atan((V[1]-rs)/r0); // switch in .nb file // switch0->1 as r->infinity
252  switch2 = 1.0-switch0; // for inner radial region
253 
254  // this works because all functions are monotonic, so final result is monotonic. Also, th(x2=1)=Pi and th(x2=0)=0 as required
255  theta1 = th0*switch2 + th2*switch0; // th0 is activated for small V[1] and th2 is activated at large radii. Notice that sum of switch2+switch0=1 so normalization correct.
256  // theta1=th0;
257  theta1=th2;
258 
259 #endif
260 
261  if(0){
262  // fix_3dpoledtissue.nb based:
263  theta2 = M_PI*0.5*(htheta*(2.0*X[2]-1.0)+(1.0-htheta)*pow(2.0*X[2]-1.0,ntheta)+1.0);
264 
265  // generate interpolation factor
266  arctan2 = 0.5 + 1.0/M_PI*(atan( (V[1]-rsjet2)/r0jet2) );
267 
268  // now interpolate between them
269  V[2] = theta2 + arctan2*(theta1-theta2);
270  }
271 
272 
273  //V[2] = theta1;
274 
275  if(1){
276  FTYPE fraceq=0.3;
277  FTYPE fracpole=(1.0-fraceq)/2.0;
278  FTYPE x2p1=0.0+fracpole;
279  FTYPE x2p2=1.0-fracpole;
280  FTYPE swide=0.04; //1E-1;
281 
282  // s(x) = 0.5 + 0.5 tanh((x-a)/b)
283 
284  // FTYPE switchh0 = 0.5+1.0/M_PI*atan((X[2]-x2p1)/swide);
285  // FTYPE switchh2 = 0.5+1.0/M_PI*atan((X[2]-x2p2)/swide);
286 
287  FTYPE switchh0 = 0.5+0.5*tanh((X[2]-x2p1)/swide);
288  FTYPE switchh2 = 0.5+0.5*tanh((X[2]-x2p2)/swide);
289 
290  FTYPE eqh=0.1;
291  FTYPE theq = M_PI * X[2] + ((1. - eqh) * 0.5) * mysin(2. * M_PI * X[2]);
292 
293  FTYPE thup=switchh0*theq + (1.0-switchh0)*theta1;
294 
295  V[2]=switchh2*theta1 + (1.0-switchh2)*thup;
296 
297  }
298  else{
299  V[2]=theta1;
300  }
301 
302  if(0){ // Sam Gralla
303  FTYPE transwidth=0.06; //1E-1;
304  FTYPE xcent=0.5; // fixed
305  FTYPE transR=0.5*(1.0-tanh(+(X[2] - xcent)/transwidth));
306  FTYPE transL=0.5*(1.0-tanh(-(X[2] - xcent)/transwidth));
307 
308  h0=0.0;
309  // FTYPE wpar=h0 + pow( (V[1]-rsjet3)/r0jet3 , -njet);
310  FTYPE wpar=pow( (V[1])/1 , -njet);
311 
312  FTYPE line1 = wpar*X[2];
313  FTYPE line2 = M_PI + wpar*(X[2]-1.0);
314 
315  V[2] = line1*transR + line2*transL;
316 
317  // dualfprintf(fail_file,"X[2]=%g V[2]=%g\n",X[2],V[2]);
318 
319  }
320 
321  if(1){ // Sam Gralla 2
322 
323  h0=0.0;
324 
325 #define cr(x) (exp(-1.0/(x)))
326 #define tr(x) (cr(x)/(cr(x) + cr(1.0-(x))))
327 #define trans(x,L,R) ((x)<=(L) ? 0.0 : ((x)>=(R) ? 1.0 : tr(((x)-(L))/((R)-(L)))) )
328 #define transR(x,center,width) ( 0.5*(1.0-tanh(+((x)-(center))/(width))))
329 #define transL(x,center,width) ( 0.5*(1.0-tanh(-((x)-(center))/(width))))
330 #define transM(x,center,width) ( exp(-pow(((x)-(center))/((width)*0.5),2.0) ) )
331 #define line1(x,w) ((x)*(w))
332 #define line2(x,w) ((x)*(w)+M_PI-(w))
333 #define line3(x,w) ((x)*(w))
334  //#define wparsam(x,r) (h0 + pow( ((r)-rsjet3)/r0jet3 , -njet))
335  //#define wparsam(x,r) (h0 + pow( ((r)-0.0)/4.2 , -njet))
336 #define wparsam(x,r) (h0 + pow(0.15 + ((r)-0.0)/10.0 , -njet))
337 #define thetasam(x,r,w,xp1,xp2) (line1(x,w)*(1.0-trans(x,xp1,xp2)) + line2(x,w)*trans(x,xp1,xp2))
338 
339  V[2] = thetasam(X[2],V[1],wparsam(X[2],V[1]),0.25,0.75);
340  // V[2] = thetasam(X[2],V[1],1.0/V[1],0.2,0.8);
341 
342  // dualfprintf(fail_file,"tr=%g %g %g %g\n",tr(0.5),line1(0.5,wparsam(0.5,V[1])),trans(X[2],0.2,0.8),line2(0.5,wparsam(0.5,V[1])));
343 
344 
345  }
346 
347  if(0){ // Sam Gralla 3
348 
349  h0=0.0;
350 
351 #define plateau(x,L,R,W) (trans(x,(L)-0.5*(W),(L)+0.5*(W))*(1.0-trans(x,(R)-0.5*(W),(R)+0.5*(W))))
352 #define lineeq(x,w) ((x)*(w)+(0.5*M_PI)-(0.5*w))
353 #define linepole(x,w) (line1(x,w))
354 #define thetaL(x,wp,weq,xp1,xp2) ( linepole(x,wp)*(1.0-trans(x,xp1,xp2)) + lineeq(x,weq)*trans(x,xp1,xp2) )
355 #define thetasam2(x,wp,weq,xp1,xp2) ( x<0.5 ? thetaL(x,wp,weq,xp1,xp2) : -thetaL(1.0-x,wp,weq,xp1,xp2)+M_PI )
356 
357  FTYPE poleslope=wparsam(X[2],V[1]);
358  FTYPE xpole=0.25;
359  // FTYPE eqslope=0.1;
360  FTYPE eqslope=0.5;
361  FTYPE xeq=0.5;
362  V[2] = thetasam2(X[2], poleslope, eqslope, xpole, xeq);
363 
364 
365  }
366 
367 
368 
369 
370 
371  // default is uniform \phi grid
372  V[3]=2.0*M_PI*X[3];
373  }
374 }
375 
376 
377 void dxdxp_analytic_user(FTYPE *X, FTYPE *V, FTYPE (*dxdxp)[NDIM])
378 {
379  dualfprintf(fail_file,"Should not be computing USERCOORDS analytically\n");
380  myexit(34698346);
381  dxdxp[3][3] = 2.0*M_PI;
382 
383 }
384 void set_points_user(void)
385 {
386 
387 
388  if(WHICHUSERCOORD==UNIHALFUSERCOORD){
389  startx[1] = 0.3999985081775278946780799743777598329673;
390  startx[2] = 0.;
391  startx[3] = 0.;
392 
393  FTYPE endx1=21.40529883372801383045167738115556702610;
394  dx[1] = (endx1-startx[1]) / totalsize[1];
395  dx[2] = 1. / totalsize[2];
396  dx[3] = 1.0/totalsize[3];
397  }
398 
399  if(WHICHUSERCOORD==JET6LIKEUSERCOORD){
400  startx[1] = pow(log(Rin-R0),1.0/npow);
401  startx[2] = 0.;
402  startx[3] = 0.;
403  dx[1] = (pow(log(Rout-R0),1.0/npow)-pow(log(Rin-R0),1.0/npow)) / totalsize[1];
404  dx[2] = 1. / totalsize[2];
405  dx[3] = 1.0/totalsize[3];
406 
407 #if(1)
408  startx[1] = log(Rin-R0)/npow;
409 
410  trifprintf( "ITERATIVE dx1: Rout=%21.15g R0=%21.15g npow=%21.15g cpow2=%21.15g cpow3=%21.15g npow2=%21.15g x1br=%21.15g rbr=%21.15g\n",Rout,R0,npow,cpow2,cpow3,npow2,x1br,rbr);
411 
412  FTYPE x1max0, x1max,dxmax;
413  int iter;
414  const FTYPE RELACC = NUMEPSILON*100.0;
415  const int ITERMAX = 100;
416 
417 
418  if( Rout < rbr ) {
419  x1max = log(Rout-R0)/npow;
420  }
421  else {
422  x1max0 = 1.1*x1br;
423  x1max = 1.2*x1br;
424 
425  //find the root via iterations
426  for( iter = 0; iter < ITERMAX; iter++ ) {
427 
428  // trifprintf( "iter=%d x1max=%21.15g x2max0=%21.15g\n",iter,x1max0,x1max);
429 
430  if( fabs((x1max - x1max0)/x1max) < RELACC ) {
431  break;
432  }
433  x1max0 = x1max;
434 
435  if(1){
436  dxmax= (pow( (log(Rout-R0) - npow*x1max0)/cpow2, 1./npow2 ) + x1br*1.0) - x1max0;
437  }
438  else{
439  // f-f0 = (x-x0)*dfdx -> if f=Rout -> x = (Rout-f0)/dfdx+x0
440 
441  FTYPE dVdx1=(npow + cpow2*npow2*cpow3*pow(cpow3*(x1max0-x1br*1.0),npow2-1.0)) * exp(npow*x1max0 + cpow2 * pow(cpow3*(x1max0-x1br*1.0),npow2));
442  FTYPE V0 = R0 + exp(npow*x1max0 + cpow2 * pow(cpow3*(x1max0-x1br*1.0),npow2));
443 
444  dxmax=(Rout-V0)/dVdx1; // x-x0
445 
446  dualfprintf(fail_file,"dVdx1=%g V0=%g dxmax=%g x1max=%g x1max0=%g\n",dVdx1,V0,dxmax,x1max,x1max0);
447  }
448 
449  // need a slight damping factor
450  FTYPE dampingfactor=0.5;
451  x1max = x1max0 + dampingfactor*dxmax;
452 
453 
454  }
455 
456  if( iter == ITERMAX ) {
457  trifprintf( "Error: iteration procedure for finding x1max has not converged: x1max = %g, dx1max/x1max = %g, iter = %d\n",
458  x1max, (x1max-x1max0)/x1max, iter );
459  exit(1);
460  }
461  else {
462  trifprintf( "x1max = %g (dx1max/x1max = %g, itno = %d)\n", x1max, (x1max-x1max0)/x1max, iter );
463  }
464  }
465 
466  dx[1] = ( x1max - startx[1] ) /totalsize[1];
467 #endif
468 
469 
470  }
471 
472 }
473 FTYPE setRin_user(int ihor, FTYPE ihoradjust)
474 {
475  if(1){
476  FTYPE ftemp;
477 
478  // see jet3coords_checknew.nb (and fix_3dpolestissue.nb) to have chosen Rin and ihor and compute required R0
479  if(npow==1.0){
480  ftemp=ihoradjust/(FTYPE)totalsize[1];
481  return(R0+pow((Rhor-R0)/pow(Rout-R0,ftemp),1.0/(1.0-ftemp)));
482  }
483  else if(npow2>0){
484  return(1.2);
485  }
486  else{
487  dualfprintf(fail_file,"ihoradjust=%21.15g totalsize[1]=%d Rhor=%21.15g R0=%21.15g npow=%21.15g Rout=%21.15g\n",ihoradjust,totalsize[1],Rhor,R0,npow,Rout);
488  return(R0+exp( pow((totalsize[1]*pow(log(Rhor-R0),1.0/npow) - ihoradjust*pow(log(Rout-R0),1.0/npow))/(totalsize[1]-ihoradjust),npow)));
489  }
490  }
491 }
492 
Generated on Fri May 20 2016 15:52:34 for HARM by   doxygen 1.8.3.1