BaseCase.cpp 18.4 KB
Newer Older
Christopher Subich's avatar
Christopher Subich committed
1
#include "BaseCase.hpp"
2
#include "Science.hpp"
Christopher Subich's avatar
Christopher Subich committed
3
4
#include "NSIntegrator.hpp"
#include "TArray.hpp"
5
#include <blitz/array.h>
6
#include <fstream>
Christopher Subich's avatar
Christopher Subich committed
7
8
9
10
11
12

//using namespace TArray;
using namespace NSIntegrator;
using blitz::Array;
using std::vector;

13
14
/* Call the source code writing function in the constructor */
extern "C" {
15
    void WriteCaseFileSource(void);
16
17
18
}
BaseCase::BaseCase(void)
{
19
    if (master()) WriteCaseFileSource();
20
}
Christopher Subich's avatar
Christopher Subich committed
21

22
/* Implementation of non-abstract methods in BaseCase */
Christopher Subich's avatar
Christopher Subich committed
23
24
25
int BaseCase::numActive() const { return 0; }
int BaseCase::numPassive() const { return 0; }
int BaseCase::numtracers() const { /* total number of tracers */
26
    return numActive() + numPassive();
Christopher Subich's avatar
Christopher Subich committed
27
28
29
}

int BaseCase::size_x() const {
30
    return size_cube();
Christopher Subich's avatar
Christopher Subich committed
31
32
}
int BaseCase::size_y() const {
33
    return size_cube();
Christopher Subich's avatar
Christopher Subich committed
34
35
}
int BaseCase::size_z() const {
36
    return size_cube();
Christopher Subich's avatar
Christopher Subich committed
37
38
}
double BaseCase::length_x() const {
39
    return length_cube();
Christopher Subich's avatar
Christopher Subich committed
40
41
}
double BaseCase::length_y() const {
42
    return length_cube();
Christopher Subich's avatar
Christopher Subich committed
43
44
}
double BaseCase::length_z() const {
45
    return length_cube();
Christopher Subich's avatar
Christopher Subich committed
46
47
48
}

DIMTYPE BaseCase::type_x() const {
49
    return type_default();
Christopher Subich's avatar
Christopher Subich committed
50
51
}
DIMTYPE BaseCase::type_y() const {
52
    return type_default();
Christopher Subich's avatar
Christopher Subich committed
53
54
}
DIMTYPE BaseCase::type_z() const {
55
    return type_default();
Christopher Subich's avatar
Christopher Subich committed
56
57
}
DIMTYPE BaseCase::type_default() const {
58
    return PERIODIC;
Christopher Subich's avatar
Christopher Subich committed
59
60
61
}

void BaseCase::tracer_bc_x(int t_num, double & dir, double & neu) const {
62
63
64
65
66
67
68
69
70
    if (!zero_tracer_boundary) {
        dir = 0; 
        neu = 1;
    }
    else {
        dir = 1;
        neu = 0;
    }
    return;
Christopher Subich's avatar
Christopher Subich committed
71
72
}
void BaseCase::tracer_bc_y(int t_num, double & dir, double & neu) const {
73
74
75
76
77
78
79
80
81
    if (!zero_tracer_boundary) {
        dir = 0; 
        neu = 1;
    }
    else {
        dir = 1;
        neu = 0;
    }
    return;
Christopher Subich's avatar
Christopher Subich committed
82
83
}
void BaseCase::tracer_bc_z(int t_num, double & dir, double & neu) const {
84
85
86
87
88
89
90
91
92
    if (!zero_tracer_boundary) {
        dir = 0; 
        neu = 1;
    }
    else {
        dir = 1;
        neu = 0;
    }
    return;
Christopher Subich's avatar
Christopher Subich committed
93
94
}
bool BaseCase::tracer_bc_forcing() const {
95
    return false;
Christopher Subich's avatar
Christopher Subich committed
96
97
}
bool BaseCase::is_mapped() const { // Whether this problem has mapped coordinates
98
    return false;
Christopher Subich's avatar
Christopher Subich committed
99
100
101
102
103
}
// Coordinate mapping proper, if is_mapped() returns true.  This features full,
// 3D arrays, but at least initially we're restricting ourselves to 2D (x,z)
// mappings
void BaseCase::do_mapping(DTArray & xgrid, DTArray & ygrid, DTArray & zgrid) {
104
    return;
Christopher Subich's avatar
Christopher Subich committed
105
106
107
}

/* Physical parameters */
108
109
110
111
112
double BaseCase::get_visco()            const { return 0; }
double BaseCase::get_diffusivity(int t) const { return 0; }
double BaseCase::get_rot_f()            const { return 0; }
int BaseCase::get_restart_sequence()    const { return 0; }
double BaseCase::get_next_plot()              { return 0; }
Christopher Subich's avatar
Christopher Subich committed
113
114
115

/* Initialization */
double BaseCase::init_time() const {
116
    return 0;
Christopher Subich's avatar
Christopher Subich committed
117
118
}
void BaseCase::init_tracers(vector<DTArray *> & tracers) {
119
120
121
122
123
124
    /* Initalize tracers one-by-one */
    if (numtracers() == 0) return; // No tracers, do nothing
    assert(numtracers() == int(tracers.size())); // Sanity check
    for (int i = 0; i < numtracers(); i++) {
        init_tracer(i, *(tracers[i]));
    }
Christopher Subich's avatar
Christopher Subich committed
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
double BaseCase::check_timestep(double t_step, double now) {
    // Check time step
    if (t_step < 1e-9) {
        // Timestep's too small, somehow stuff is blowing up
        if (master()) fprintf(stderr,"Tiny timestep (%e), aborting\n",t_step);
        return -1;
    } else if (t_step > get_dt_max()) {
        // Cap the maximum timestep size
        t_step = get_dt_max();
    }

    // Now, calculate how many timesteps remain until the next writeout
    double until_plot = get_next_plot() - now;
    int steps = ceil(until_plot / t_step);
    // Where will we be after (steps) timesteps of the current size?
    double real_until_plot = steps*t_step;
    // If that's close enough to the real writeout time, that's fine.
    if (fabs(until_plot - real_until_plot) < 1e-6) {
        return t_step;
    } else {
        // Otherwise, square up the timeteps.  This will always shrink the timestep.
    return (until_plot / steps);
    }
Christopher Subich's avatar
Christopher Subich committed
150
151
152
153
}

/* Forcing */
void BaseCase::forcing(double t,  DTArray & u, DTArray & u_f,
154
155
156
157
158
159
160
161
162
163
164
165
166
167
        DTArray & v, DTArray & v_f,  DTArray & w, DTArray & w_f,
        vector<DTArray *> & tracers, vector<DTArray *> & tracers_f) {
    /* First, if no active tracers then use the simpler form */
    if (numActive() == 0) {
        passive_forcing(t, u, u_f, v, v_f, w, w_f);
    } else {
        /* Look at split velocity-tracer forcing */
        vel_forcing(t, u_f, v_f, w_f, tracers);
        tracer_forcing(t, u, v, w, tracers_f);
    }
    /* And any/all passive tracers get 0 forcing */
    for (int i = 0; i < numPassive(); i++) {
        *(tracers_f[numActive() + i]) = 0;
    }
Christopher Subich's avatar
Christopher Subich committed
168
169
}
void BaseCase::passive_forcing(double t,  DTArray & u, DTArray & u_f,
170
171
172
        DTArray & v, DTArray & v_f,  DTArray &, DTArray & w_f) {
    /* Reduce to velocity-independent case */
    stationary_forcing(t, u_f, v_f, w_f);
Christopher Subich's avatar
Christopher Subich committed
173
174
}
void BaseCase::stationary_forcing(double t, DTArray & u_f, DTArray & v_f, 
175
176
177
178
179
        DTArray & w_f) {
    /* Default case, 0 forcing */
    u_f = 0;
    v_f = 0;
    w_f = 0;
Christopher Subich's avatar
Christopher Subich committed
180
181
182
183
}

/* Analysis */
void BaseCase::analysis(double t, DTArray & u, DTArray & v, DTArray & w,
184
185
        vector<DTArray *> tracer, DTArray & pres) {
    analysis(t,u,v,w,tracer);
Christopher Subich's avatar
Christopher Subich committed
186
187
}
void BaseCase::analysis(double t, DTArray & u, DTArray & v, DTArray & w,
188
189
190
191
192
193
194
195
        vector<DTArray *> tracer) {
    /* Do velocity and tracer analysis seperately */
    vel_analysis(t, u, v, w);
    for (int i = 0; i < numtracers(); i++) {
        tracer_analysis(t, i, *(tracer[i]));
    } 
}

196
void BaseCase::automatic_grid(double MinX, double MinY, double MinZ,
197
        Array<double,1> * xx, Array<double,1> * yy, Array<double,1> * zz){
198
    //Array<double,1> xx(split_range(size_x())), yy(size_y()), zz(size_z());
199
200
201
202
203
204
205
206
207
208
209
210
211
    bool xxa = false, yya = false, zza = false;
    if (!xx) {
        xxa = true; // Delete xx when returning
        xx = new Array<double,1>(split_range(size_x()));
    }
    if (!yy) {
        yya = true;
        yy = new Array<double,1>(size_y());
    }
    if (!zz) {
        zza = true;
        zz = new Array<double,1>(size_z());
    }
212
213
214
215
216
217
218
219
220
    Array<double,3> grid(alloc_lbound(size_x(),size_y(),size_z()),
            alloc_extent(size_x(),size_y(),size_z()),
            alloc_storage(size_x(),size_y(),size_z()));
    blitz::firstIndex ii;
    blitz::secondIndex jj;
    blitz::thirdIndex kk;

    // Generate 1D arrays
    if (type_x() == NO_SLIP) {
221
        *xx = MinX+length_x()*(0.5-0.5*cos(M_PI*ii/(size_x()-1)));
222
    } else {
223
        *xx = MinX + length_x()*(ii+0.5)/size_x();
224
    }
225
    *yy = MinY + length_y()*(ii+0.5)/size_y();
226
    if (type_z() == NO_SLIP) {
227
        *zz = MinZ+length_z()*(0.5-0.5*cos(M_PI*ii/(size_z()-1)));
228
    } else {
229
        *zz = MinZ + length_z()*(0.5+ii)/size_z();
230
231
232
    }

    // Write grid/reader
233
    grid = (*xx)(ii) + 0*jj + 0*kk;
234
235
236
237
    write_array(grid,"xgrid");
    write_reader(grid,"xgrid",false);

    if (size_y() > 1) {
238
        grid = 0*ii + (*yy)(jj) + 0*kk;
239
240
241
242
        write_array(grid,"ygrid");
        write_reader(grid,"ygrid",false);
    }

243
    grid = 0*ii + 0*jj + (*zz)(kk);
244
245
    write_array(grid,"zgrid");
    write_reader(grid,"zgrid",false);
246
247
248
249
250

    // Clean up
    if (xxa) delete xx;
    if (yya) delete yy;
    if (zza) delete zz;
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
}

/* Read velocities from regular output */
void BaseCase::init_vels_restart(DTArray & u, DTArray & v, DTArray & w){
    /* Restarting, so build the proper filenames and load the data into u, v, w */
    char filename[100];

    /* u */
    snprintf(filename,100,"u.%d",get_restart_sequence());
    if (master()) fprintf(stdout,"Reading u from %s\n",filename);
    read_array(u,filename,size_x(),size_y(),size_z());

    /* v, only necessary if this is an actual 3D run or if
       rotation is noNzero */
    if (size_y() > 1 || get_rot_f() != 0) {
        snprintf(filename,100,"v.%d",get_restart_sequence());
        if (master()) fprintf(stdout,"Reading v from %s\n",filename);
        read_array(v,filename,size_x(),size_y(),size_z());
    }

    /* w */
    snprintf(filename,100,"w.%d",get_restart_sequence());
    if (master()) fprintf(stdout,"Reading w from %s\n",filename);
    read_array(w,filename,size_x(),size_y(),size_z());
    return;
}

/* Read velocities from dump output */
void BaseCase::init_vels_dump(DTArray & u, DTArray & v, DTArray & w){
    /* Restarting, so build the proper filenames and load the data into u, v, w */

    /* u */
    if (master()) fprintf(stdout,"Reading u from u.dump\n");
    read_array(u,"u.dump",size_x(),size_y(),size_z());

    /* v, only necessary if this is an actual 3D run or if
       rotation is noNzero */
    if (size_y() > 1 || get_rot_f() != 0) {
        if (master()) fprintf(stdout,"Reading v from v.dump\n");
        read_array(v,"v.dump",size_x(),size_y(),size_z());
    }

    /* w */
    if (master()) fprintf(stdout,"Reading w from w.dump\n");
    read_array(w,"w.dump",size_x(),size_y(),size_z());
    return;
}


/* Read field from regular output */
void BaseCase::init_tracer_restart(const std::string & field, DTArray & the_tracer){
    /* Restarting, so build the proper filenames and load the data */
    char filename[100];

    snprintf(filename,100,"%s.%d",field.c_str(),get_restart_sequence());
    if (master()) fprintf(stdout,"Reading %s from %s\n",field.c_str(),filename);
    read_array(the_tracer,filename,size_x(),size_y(),size_z());
    return;
}

/* Read field from dump output */
void BaseCase::init_tracer_dump(const std::string & field, DTArray & the_tracer){
    /* Restarting, so build the proper filenames and load the data */
    char filename[100];

    snprintf(filename,100,"%s.dump",field.c_str());
    if (master()) fprintf(stdout,"Reading %s from %s\n",field.c_str(),filename);
    read_array(the_tracer,filename,size_x(),size_y(),size_z());
    return;
}

322
323
324
325
326
327
328
329
330
331
332
333
334
/* write out vertical chain of data */
void BaseCase::write_chain(const char *filename, DTArray & val, int Iout, int Jout, double time) {
    FILE *fid=fopen(filename,"a");
    if (fid == NULL) {
        fprintf(stderr,"Unable to open %s for writing\n",filename);
        MPI_Finalize(); exit(1);
    }
    fprintf(fid,"%g",time);
    for (int ki=0; ki<size_z(); ki++) fprintf(fid," %g",val(Iout,Jout,ki));
    fprintf(fid,"\n");
    fclose(fid);
}

335
336
337
338
339
340
341
342
/* Read grid from regular output */
void BaseCase::init_grid_restart(const std::string & component,
        const std::string & filename, DTArray & grid) {
    if (master()) fprintf(stdout,"Reading %s from %s\n",component.c_str(),filename.c_str());
    read_array(grid,filename.c_str(),size_x(),size_y(),size_z());
    return;
}

343

344
345
346
347
/* Check and dump */
void BaseCase::check_and_dump(double clock_time, double real_start_time,
        double compute_time, double sim_time, double avg_write_time, int plot_number,
        DTArray & u, DTArray & v, DTArray & w, vector<DTArray *> & tracer){
348
349
350
    if (compute_time > 0) {
        // default is to not dump variables
        int do_dump = 0;
351
352

        // check if close to end of compute time
353
354
355
356
357
        if (master()) {
            double total_run_time = clock_time - real_start_time;
            if (compute_time - total_run_time < 3*avg_write_time) {
                do_dump = 1; // true
            }
358
359
        }

360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
        // Broadcast to all processors whether to dump or not
        MPI_Bcast(&do_dump, 1, MPI_INT, 0, MPI_COMM_WORLD);

        // Dump variables if close to end of compute time
        if (do_dump == 1) {
            if (master()){
                fprintf(stdout,"Too close to final time, dumping!\n");
            }
            write_variables(u, v, w, tracer);

            // Write the dump time to a text file for restart purposes
            if (master()){
                FILE * dump_file;
                dump_file = fopen("dump_time.txt","w");
                assert(dump_file);
                fprintf(dump_file,"The dump time was:\n%.17g\n", sim_time);
                fprintf(dump_file,"The dump index was:\n%d\n", plot_number);
                fclose(dump_file);
            }

            // Die gracefully
            MPI_Finalize(); exit(0);
382
383
        }
    }
384
385
}

386
387
388
389
390
391
392
393
394
395
396
397
/* Change dump log file for successful completion*/
void BaseCase::successful_dump(int plot_number, double final_time, double plot_interval){
    if (master() && (plot_number == final_time/plot_interval)){
        // Write the dump time to a text file for restart purposes
        FILE * dump_file; 
        dump_file = fopen("dump_time.txt","w");
        assert(dump_file);
        fprintf(dump_file,"The dump 'time' was:\n%.12g\n", 2*final_time);
        fprintf(dump_file,"The dump index was:\n%d\n", plot_number);
        fclose(dump_file);
    }
}
398

399
400
// Explicitly instantiate common add_diagnostic templates, for integer and double
// types
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
template void BaseCase::add_diagnostic<int>(const string str, const int val,
        string & header, string & line);
template void BaseCase::add_diagnostic<double>(const string str, const double val,
        string & header, string & line);

// write the diagnostic file
void BaseCase::write_diagnostics(string header, string line,
        int iter, bool restarting) {
    // remove last two elements (comma and space)
    string clean_header = header.substr(0, header.size()-2);
    string clean_line   =   line.substr(0,   line.size()-2);
    // open file
    FILE * diagnos_file = fopen("diagnostics.txt","a");
    assert(diagnos_file);
    // print header
    if (iter == 1 and !restarting) {
        fprintf(diagnos_file,"%s\n",clean_header.c_str());
    }
    // print the line of values
    fprintf(diagnos_file, "%s\n", clean_line.c_str());
    // Close file
    fclose(diagnos_file);
}
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455

// parse expansion types 
void parse_boundary_conditions(const string xgrid_type, const string ygrid_type,
        const string zgrid_type, DIMTYPE & intype_x, DIMTYPE & intype_y, DIMTYPE & intype_z) {
    // x
    if      (xgrid_type == "FOURIER")   { intype_x = PERIODIC; }
    else if (xgrid_type == "FREE_SLIP") { intype_x = FREE_SLIP; }
    else if (xgrid_type == "NO_SLIP")   { intype_x = NO_SLIP; }
    else {
        if (master())
            fprintf(stderr,"Invalid option %s received for type_x\n",xgrid_type.c_str());
        MPI_Finalize(); exit(1);
    }
    // y
    if      (ygrid_type == "FOURIER")   { intype_y = PERIODIC; }
    else if (ygrid_type == "FREE_SLIP") { intype_y = FREE_SLIP; }
    else {
        if (master())
            fprintf(stderr,"Invalid option %s received for type_y\n",ygrid_type.c_str());
        MPI_Finalize(); exit(1);
    }
    // z
    if      (zgrid_type == "FOURIER")   { intype_z = PERIODIC; }
    else if (zgrid_type == "FREE_SLIP") { intype_z = FREE_SLIP; }
    else if (zgrid_type == "NO_SLIP")   { intype_z = NO_SLIP; }
    else {
        if (master())
            fprintf(stderr,"Invalid option %s received for type_z\n",zgrid_type.c_str());
        MPI_Finalize(); exit(1);
    }
}

456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
// Top stresses
void BaseCase::stresses_top(TArrayn::DTArray & u, TArrayn::DTArray & v, TArrayn::DTArray & w,
        TArrayn::DTArray & Hprime, TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
        const string * grid_type, const double mu, double time, int itercount, bool restarting) {
    // set-up
    static DTArray *tx = alloc_array(size_x(),size_y(),1);
    static DTArray *ty = alloc_array(size_x(),size_y(),1);
    blitz::firstIndex ii;
    blitz::secondIndex jj;

    // top stress ( along channel - x )
    top_stress_x(*tx, u, temp, gradient_op, grid_type, size_z(), mu);
    double top_tx_tot = pssum(sum(
                (*get_quad_x())(ii)*
                (*get_quad_y())(jj)*(*tx)));
    double top_tx_abs = pssum(sum(
                (*get_quad_x())(ii)*
                (*get_quad_y())(jj)*abs(*tx)));
    // top stress ( across channel - y )
    top_stress_y(*ty, v, temp, gradient_op, grid_type, size_z(), mu);
    double top_ty_tot = pssum(sum(
                (*get_quad_x())(ii)*
                (*get_quad_y())(jj)*(*ty)));
    double top_ty_abs = pssum(sum(
                (*get_quad_x())(ii)*
                (*get_quad_y())(jj)*abs(*ty)));
    // total top stress
    double top_ts = pssum(sum(
                (*get_quad_x())(ii)*
                (*get_quad_y())(jj)*pow(pow(*tx,2)+pow(*ty,2),0.5)));
    // write to a stress diagnostic file
    if (master()) {
        FILE * stresses_file = fopen("stresses_top.txt","a");
        assert(stresses_file);
        if ( itercount==1 and !restarting )
            fprintf(stresses_file,"Time, "
                    "Top_tx_tot, Top_tx_abs, Top_ty_tot, Top_ty_abs, Top_ts\n");
        fprintf(stresses_file,"%.17f, "
                "%.17g, %.17g, %.17g, %.17g, %.17g\n",
                time,
                top_tx_tot, top_tx_abs, top_ty_tot, top_ty_abs, top_ts);
        fclose(stresses_file);
    }
}

// Bottom stresses
void BaseCase::stresses_bottom(TArrayn::DTArray & u, TArrayn::DTArray & v, TArrayn::DTArray & w,
        TArrayn::DTArray & Hprime, TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
        const string * grid_type, const double mu, double time, int itercount, bool restarting) {
    // set-up
    static DTArray *tx = alloc_array(size_x(),size_y(),1);
    static DTArray *ty = alloc_array(size_x(),size_y(),1);
    blitz::firstIndex ii;
    blitz::secondIndex jj;

    // bottom stress ( along channel - x )
    bottom_stress_x(*tx, Hprime, u, w, temp, gradient_op, grid_type, is_mapped(), mu);
    double bot_tx_tot = pssum(sum(
                (*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
                (*get_quad_y())(jj)*(*tx)));
    double bot_tx_abs = pssum(sum(
                (*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
                (*get_quad_y())(jj)*abs(*tx)));
    // bottom stress ( across channel - y )
    bottom_stress_y(*ty, Hprime, v, temp, gradient_op, grid_type, is_mapped(), mu);
    double bot_ty_tot = pssum(sum(
                (*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
                (*get_quad_y())(jj)*(*ty)));
    double bot_ty_abs = pssum(sum(
                (*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
                (*get_quad_y())(jj)*abs(*ty)));
    // total bottom stress
    double bot_ts = pssum(sum(
                (*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
                (*get_quad_y())(jj)*pow(pow(*tx,2)+pow(*ty,2),0.5)));

    // write to a stress diagnostic file
    if (master()) {
        FILE * stresses_file = fopen("stresses_bottom.txt","a");
        assert(stresses_file);
        if ( itercount==1 and !restarting )
            fprintf(stresses_file,"Time, "
                    "Bottom_tx_tot, Bottom_tx_abs, Bottom_ty_tot, Bottom_ty_abs, Bottom_ts\n");
        fprintf(stresses_file,"%.17f, "
                "%.17g, %.17g, %.17g, %.17g, %.17g\n",
                time,
                bot_tx_tot, bot_tx_abs, bot_ty_tot, bot_ty_abs, bot_ts);
        fclose(stresses_file);
    }
}