Disco/analysis.c at master · chomps/Disco · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

#include "paul.h"

int num_diagnostics(void){
   return(5);
}

void planetaryForce( struct planet * , double , double , double , double * , double * , double * , int );

void get_diagnostics( double * x , double * prim , double * Q , struct domain * theDomain ){
   double r = x[0];
   double phi = x[1];

   double rho = prim[RHO];
   double vr = prim[URR];
   double omega = prim[UPP];

   Q[0] = rho;
   Q[1] = 2.*M_PI*r*rho*vr;
   double Fr,Fp,Fz;
   Fp = 0.0;
   if( theDomain->Npl > 1 ){
      struct planet * pl = theDomain->thePlanets+1;
      planetaryForce( pl , r , phi , 0.0 , &Fr , &Fp , &Fz , 0 );
   }
   Q[2] = 2.*M_PI*r*rho*(r*Fp);
   Q[3] = 2.*M_PI*r*rho*( r*r*omega*vr );
   Q[4] = omega;
}

void zero_diagnostics( struct domain * theDomain ){

   int Nr = theDomain->Nr;
   int Nq = theDomain->num_tools;
   struct diagnostic_avg * theTools = &(theDomain->theTools);

   int i,q;
   for( i=0 ; i<Nr ; ++i ){
      for( q=0 ; q<Nq ; ++q ){
         int iq = i*Nq + q;
         theTools->Qr[iq] = 0.0;
      }
   }
   theTools->t_avg = 0.0;

}

void avg_diagnostics( struct domain * theDomain ){

   int Nr = theDomain->Nr;
   int Nq = theDomain->num_tools;
   struct diagnostic_avg * theTools = &(theDomain->theTools);
   double dt = theTools->t_avg;

   int i,q;
   for( i=0 ; i<Nr ; ++i ){
      for( q=0 ; q<Nq ; ++q ){
         int iq = i*Nq + q;
         theTools->Qr[iq] /= dt;
      }
   }
   theTools->t_avg = 0.0;

}

double get_dV( double * , double * );

void add_diagnostics( struct domain * theDomain , double dt ){

   int Nr = theDomain->Nr;
   int Nz = theDomain->Nz;
   int * Np = theDomain->Np;
   int Nq = theDomain->num_tools;
   struct diagnostic_avg * theTools = &(theDomain->theTools);
   struct cell ** theCells = theDomain->theCells;
   double * r_jph = theDomain->r_jph;
   double * z_kph = theDomain->z_kph;

   double temp_sum[Nr*Nq];
   memset( temp_sum , 0 , Nr*Nq*sizeof(double) );
   int i,j,k,q;
   int kmin = 0;
   int kmax = Nz;

   for( j=0 ; j<Nr ; ++j ){
      double dV_tot = 0.0;
      for( k=kmin ; k<kmax ; ++k ){
         int jk = j+Nr*k;
         for( i=0 ; i<Np[jk] ; ++i ){
            struct cell * c = theCells[jk]+i;
            double phip = c->piph;
            double phim = phip - c->dphi;
            double xp[3] = {r_jph[j  ] , phip , z_kph[k  ]};
            double xm[3] = {r_jph[j-1] , phim , z_kph[k-1]};
            double xc[3];
            for( q=0 ; q<3 ; ++q ) xc[q] = .5*(xp[q]+xm[q]);
            double dV = get_dV(xp,xm);
            double Q[Nq];
            get_diagnostics( xc , c->prim , Q , theDomain );
            for( q=0 ; q<Nq ; ++q ) temp_sum[ j*Nq + q ] += Q[q]*dV;
            dV_tot += dV;
         }
      }
      for( q=0 ; q<Nq ; ++q ){
         theTools->Qr[ j*Nq + q ] += temp_sum[ j*Nq + q ]*dt/dV_tot;
      }
   }

   theTools->t_avg += dt;

}