4 void Init (
double *us,
double x1,
double x2,
double x3)
11 static int first_call = 1;
12 double x, y, R, Rc, phi,
c,
s;
29 x = R*c - Rc/sqrt(2.0);
30 y = R*s - Rc/sqrt(2.0);
39 us[
VX2] = 1.0/sqrt(R);
41 arg = - (x*x + y*y)/(h*h);
46 us[
VX1] += dvx*c + dvy*
s;
47 us[
VX2] += -dvx*s + dvy*
c;
55 us[
TRC] = (exp(arg) > 1.e-2);
68 static int first_call = 1;
70 double *R, *phi, *dVR, *dVphi, *dR, *dphi, vol;
72 double E=0.0, w, pw, wmin, pwmin;
73 double Ltot, Etot, wtot, pwtot;
89 MPI_Allreduce (&voltot, &vol, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
90 MPI_Barrier (MPI_COMM_WORLD);
101 Etot = Ltot = wtot = pwtot = 0.0;
102 wmin = pwmin = 1.e20;
104 vol = dVR[
i]*dVphi[
j]/voltot;
106 for (nv =
NVAR; nv--; ) v[nv] = d->
Vc[nv][k][j][i];
108 w = ( (vphi[j][i+1] - 1.0/sqrt(R[i+1]))*R[i+1]
109 - (vphi[j][i-1] - 1.0/sqrt(R[i-1]))*R[i-1])/(2.0*R[
i]*dR[
i]);
110 w -= (vR[j+1][
i] - vR[j-1][
i])/(2.0*R[i]*dphi[j]);
121 Ltot += R[
i]*v[
RHO]*v[
VX2]*vol;
127 pwmin =
MIN(pwmin, pw);
131 MPI_Allreduce (&Ltot, &vol, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
134 MPI_Allreduce (&Etot, &vol, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
137 MPI_Allreduce (&wtot, &vol, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
140 MPI_Allreduce (&pwtot, &vol, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
143 MPI_Allreduce (&wmin, &vol, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
146 MPI_Allreduce (&pwmin, &vol, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
149 MPI_Barrier (MPI_COMM_WORLD);
161 static double tpos = -1.0;
163 fp = fopen(
"averages.dat",
"w");
164 fprintf (fp,
"#%s %12s %13s %13s %13s %14s %12s %12s\n",
165 " t",
"dt",
"<L>",
"<E>",
"<w>",
"<w/rho>",
"min(w)",
"min(pw)");
169 fp = fopen(
"averages.dat",
"r");
170 while (fgets(sline, 512, fp)) {}
171 sscanf(sline,
"%lf\n",&tpos);
174 fp = fopen(
"averages.dat",
"a");
177 fprintf (fp,
"%12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n",
178 g_time,
g_dt, Ltot, Etot, wtot, pwtot, wmin, pwmin);
223 #if ROTATING_FRAME == YES
232 }
else if (side ==
X1_END) {
243 #if ROTATING_FRAME == YES
251 #if (BODY_FORCE & VECTOR)
253 void BodyForceVector(
double *v,
double *g,
double x1,
double x2,
double x3)
260 g[
IDIR] = -1.0/(x1*x1);
266 #if (BODY_FORCE & POTENTIAL)
#define X1_BEG
Boundary region at X1 beg.
void UserDefBoundary(const Data *d, RBox *box, int side, Grid *grid)
int vpos
Location of the variable inside the cell.
#define BOX_LOOP(B, k, j, i)
double **** Vc
The main four-index data array used for cell-centered primitive variables.
double g_dt
The current integration time step.
#define X1_END
Boundary region at X1 end.
long int g_stepNumber
Gives the current integration step number.
double g_inputParam[32]
Array containing the user-defined parameters.
long int IEND
Upper grid index of the computational domain in the the X1 direction for the local processor...
double g_time
The current integration time.
double BodyForcePotential(double x, double y, double z)
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
void Analysis(const Data *d, Grid *grid)
void Init(double *v, double x1, double x2, double x3)
void BodyForceVector(double *v, double *g, double x, double y, double z)