doc_map_const.cpp 8.04 KB
Newer Older
Christopher Subich's avatar
Christopher Subich committed
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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
/* A sample case, for illustrating Benard convection cells */

// Required headers 
#include <blitz/array.h> // Blitz++ array library
#include "../TArray.hpp" // Custom extensions to the library to support FFTs
#include "../NSIntegrator.hpp" // Time-integrator for the Navier-Stokes equations
#include <mpi.h> // MPI parallel library
#include "../BaseCase.hpp" // Support file that contains default implementations of several functions

using namespace std;
using namespace NSIntegrator;

// Tensor variables for indexing
blitz::firstIndex ii;
blitz::secondIndex jj;
blitz::thirdIndex kk;

// Physical constants
const double g = 9.81;
const double rho_0 = 1028; // Units of kg / m^3

// Pysical parameters
const double ROT_F = 0.5e-4; // 1/s, mid-latitude
const double LENGTH_X = 4e5; // 1000km
const double LENGTH_Z = 5000; // Water depth
const double N2_max = 1e-3; // Linear stratification buoyancy frequency

const double S0 = 35; // Baseline salinity
const double beta = 7.6e-4; // Density change from salt content

// Hill parameters
const double H_HEIGHT = 1500; // Hill height
const double H_LENGTH = 12e3; // Hill length

// Tide parameters
const double TIDE_PERIOD = 44712;
const double TIDE_M2 = 2*M_PI/TIDE_PERIOD; // M2 tidal frequency
const double TIDE_STRENGTH = 0.01; // Desired maximum tidal current

// Numerical parameters
int NZ = 0; // Number of vertical points.  Number of horizontal points
int NX = 0; // will be calculated based on this.
const double plot_interval = TIDE_PERIOD/32; // Time between field writes
const double final_time = 4*TIDE_PERIOD;


class mapiw : public BaseCase {
   public:
      // Variables to set the plot sequence number and time of the last writeout
      DTArray * xgrid, * zgrid;
      Array<double,1> hill;
      int plot_number; double last_plot;

      // Resolution in X, Y (1), and Z
      int size_x() const { return NX; }
      int size_y() const { return 1; }
      int size_z() const { return NZ; }

      /* Set periodic in x, free slip in z */
      DIMTYPE type_z() const {return NO_SLIP;}
      DIMTYPE type_default() const { return PERIODIC; }

      /* The grid size is governed through the #defines above */
      double length_x() const { return LENGTH_X; }
      double length_y() const { return 1; }
      double length_z() const { return LENGTH_Z; }

      bool is_mapped() const {return true;}
      void do_mapping(DTArray & xg, DTArray & yg, DTArray & zg) {
         xgrid = alloc_array(NX,1,NZ);
         zgrid = alloc_array(NX,1,NZ);
         Array<double,1> xx(split_range(NX)), zz(NZ);
         // Use periodic coordinates in horizontal
         xx = LENGTH_X*(ii+0.5)/NX; // x-coordinate
         zz = cos(ii*M_PI/(NZ-1)); // Chebyshev in vertical

         xg = xx(ii) + 0*jj + 0*kk;
         *xgrid = xg;

         hill = H_HEIGHT*exp(-pow((xx(ii)-LENGTH_X/2)/H_LENGTH,2));
         zg = -LENGTH_Z/2+LENGTH_Z/2*zz(kk) + 0.5*(1-zz(kk))*
                  hill(ii);;
         *zgrid = zg;

         yg = 0;

         write_array(xg,"xgrid");
         write_reader(xg,"xgrid",false);
         write_array(zg,"zgrid");
         write_reader(zg,"zgrid",false);
      }
      /* Use one actively-modified tracer */
      int numActive() const { return 1; }

      // Use 0 viscosity and diffusivity
      double get_visco() const { return 0; }
      double get_diffusivity(int t_num) const { return 0; }

      /* Start at t=0 */
      double init_time() const { return 0; }

      /* Modify the timestep if necessary in order to land evenly on a plot time */
      double check_timestep (double intime, double now) {
         // Firstly, the buoyancy frequency provides a timescale that is not
         // accounted for with the velocity-based CFL condition.
         if (intime > 0.5/sqrt(N2_max)) {
            intime = 0.5/sqrt(N2_max);
         }
         // Now, calculate how many timesteps remain until the next writeout
         double until_plot = last_plot + plot_interval - now;
         int steps = ceil(until_plot / intime);
         // And calculate where we will actually be after (steps) timesteps
         // of the current size
         double true_fintime = steps*intime;

         // If that's close enough to the real writeout time, that's fine.
         if (fabs(until_plot - true_fintime) < 1e-6) {
            return intime;
         } else {
            // Otherwise, square up the timeteps.  This will always shrink the timestep.
            return (until_plot / steps);
         }
      }


      /* Initialize velocities at the start of the run.  For this simple
         case, initialize all velocities to 0 */
      void init_vels(DTArray & u, DTArray & v, DTArray & w) {
         u = 0;
         v = -TIDE_STRENGTH*ROT_F/TIDE_M2*LENGTH_Z/(LENGTH_Z-hill(ii));
         w = 0;
         // Also, write out the (zero) initial velocities and proper M-file readers
         write_reader(u,"u",true);
         write_reader(v,"v",true);
         write_reader(w,"w",true);
         write_array(u,"u",0);
         write_array(v,"v",0);
         write_array(w,"w",0);
         return;
      }

      void init_tracer(int t_num, DTArray & salt) {
         // The primary constituent of density is salt, so that is 
         // initialize here

         salt = S0 + (N2_max*N2_max)/(-beta*g)*(*zgrid)(ii,jj,kk);
         write_array(salt,"s",0); write_reader(salt,"s",true);
      }

      // Forcing must be done generally, since both rotation and density are
      // involved
152
153
      void forcing(double t, DTArray & u, DTArray & u_f, 
            DTArray & v, DTArray & v_f, DTArray & w,
Christopher Subich's avatar
Christopher Subich committed
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
            DTArray & w_f, vector<DTArray *> & tracers,
            vector<DTArray *> & tracers_f) {
         // Rotation couples u and v, plus a source term for the tide
         u_f = -ROT_F*v + cos(t*TIDE_M2)*TIDE_STRENGTH*
            (TIDE_M2-ROT_F*ROT_F/TIDE_M2);
         v_f = ROT_F*u;
         w_f = -g*beta*(*tracers[0]-S0);
         // And since the salt content is expressed as total content rather
         // than perturbation, no forcing is necessary.
         *tracers_f[0] = 0;
      }

      /* Basic analysis, to write out the field periodically */
      void analysis(double time, DTArray & u, DTArray & v, DTArray & w,
            vector<DTArray *> & tracer, DTArray & pressure) {
         /* If it is very close to the plot time, write data fields to disk */
         if ((time - last_plot - plot_interval) > -1e-6) {
            plot_number++;
            if (master()) fprintf(stderr,"*");
            write_array(u,"u",plot_number);
            write_array(v,"v",plot_number);
            write_array(w,"w",plot_number);
            write_array(*tracer[0],"s",plot_number);
            last_plot = last_plot + plot_interval;
         }
         // Also, calculate and write out useful information: maximum u, w, and t'
         double max_u = psmax(max(abs(u)));
         double max_w = psmax(max(abs(w)));
         double max_t = psmax(max(abs(*tracer[0])));
         if (master()) fprintf(stderr,"%.2f: %.2g %.2g %.2g\n",time,max_u,max_w,max_t);
      }

      mapiw(): hill(split_range(NX))
      { // Initialize the local variables
         plot_number = 0;
         last_plot = 0;
         // Create one-dimensional arrays for the coordinates
      }

};

/* The ``main'' routine */
int main(int argc, char ** argv) {
   /* Initialize MPI.  This is required even for single-processor runs,
      since the inner routines assume some degree of parallelization,
      even if it is trivial. */
   MPI_Init(&argc, &argv);
   if (argc > 2) { // Check command line arguments
      NX = atoi(argv[1]); // Read in number of horizontal points, if specified
      NZ = atoi(argv[2]); // and vertical
   } 
   if (NX <= 0) {
      NX = 2048;
   }
   if (NZ <= 0) {
      NZ = 128;
   }
   if (master()) {
      fprintf(stderr,"Using a grid of %d x %d points\n",NX,NZ);
   }
   mapiw mycode; // Create an instantiated object of the above class
   /// Create a flow-evolver that takes its settings from the above class
   FluidEvolve<mapiw> do_mapiw(&mycode);
   // Run to a final time of 1.
   do_mapiw.initialize();
   do_mapiw.do_run(final_time);
   MPI_Finalize(); // Cleanly exit MPI
   return 0; // End the program
}