BaseCase.hpp 9.63 KB
Newer Older
Christopher Subich's avatar
Christopher Subich committed
1
2
3
4
5
6
7
8
/* BaseCase -- basic skeletal framework for a functional NSIntegrator usercode.
   Eliminates some of the boilerplate. */
#ifndef BASECASE_HPP
#define BASECASE_HPP 1

#include <blitz/array.h>
#include "TArray.hpp"
#include "NSIntegrator.hpp"
9
#include "Science.hpp"       // Science content
Christopher Subich's avatar
Christopher Subich committed
10
11
12
13
14
15
16
17
18
19
20
21
22

using namespace TArrayn;
using namespace NSIntegrator;
using blitz::Array;
using std::vector;

class BaseCase {
   /* To reduce boilerplate, wrap some of the long functions, only calling
      them if actually used by usercode.  For example, a tracer-free code
      does not need the long-form forcing function -- nor does one with
      only passive tracers */

   public:
23
24
25
      /* Constructor */
      BaseCase();

Christopher Subich's avatar
Christopher Subich committed
26
27
28
29
30
31
32
33
34
      /* Tracers */
      virtual int numActive() const; // Number of active tracers
      virtual int numPassive() const; // Number of passive tracers
      virtual int numtracers() const; // Number of tracers (total)

      /* Grid generataion */
      virtual int size_x() const; // Grid points in x
      virtual int size_y() const; // Grid points in y
      virtual int size_z() const; // Grid points in z
35
36
37
      virtual int size_cube() const { 
         assert(0 && "size_cube not implemented");
         abort(); }; // Special case -- N*N*N
Christopher Subich's avatar
Christopher Subich committed
38
39
40
41

      virtual double length_x() const; // Length in x
      virtual double length_y() const; // Length in y
      virtual double length_z() const; // Length in z
42
43
44
      virtual double length_cube() const {
         assert(0 && "length_cube not implemented");
         abort();}; // Special case -- L*L*L
Christopher Subich's avatar
Christopher Subich committed
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

      virtual DIMTYPE type_x() const; // Expansion type in x
      virtual DIMTYPE type_y() const; // Expansion type in y
      virtual DIMTYPE type_z() const; // Expansion type in z
      virtual DIMTYPE type_default() const; // Default expansion type

      /* Functions to provide for custom tracer boundaries, along any
         dimension whose expansion type supports them (generally that
         would be chebyshev-type).  The default tracer-BC behaviour
         will be to give neumann-type BCs. Dirichlet-type BCs are possible
         with a sine-based expansion rather than cosine-based expansion
         and have been hacked in with a global variable, but this mechanism
         is far more general.  One will expect an error/abort if an improper
         BC-type is requested for a boundary that doesn't actually support it
         (namely a Robin-type BC on a non-Chebyshev dimension */


      // Further note -- we want the SAME boundary condition at the minimum and
      // maximum of the domain.
      virtual void tracer_bc_x(int tracernum, double & dirichlet, double & neumann) const;
      virtual void tracer_bc_y(int tracernum, double & dirichlet, double & neumann) const;
      virtual void tracer_bc_z(int tracernum, double & dirichlet, double & neumann) const;

      // Whether ANY tracers will have nonzero boundary conditions.  If this is true
      // then the user forcing code is responsible for calculating/applying the proper
      // RHS at the boundaries.  If this is false (default), then the prior behaviour
      // of the integrator code zeroing BCs after forcing is retained.
      virtual bool tracer_bc_forcing() const;

      virtual bool is_mapped() const; // Whether this problem has mapped coordinates
      // 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
      virtual void do_mapping(DTArray & xgrid, DTArray & ygrid, DTArray & zgrid);

      /* Physical parameters */
      virtual double get_visco() const; // Physical viscosity
      virtual double get_diffusivity(int tracernum) const; // Diffusivity
83
84
      virtual double get_rot_f() const; // Physical rotation rate
      virtual int get_restart_sequence() const; // restart sequence
85
86
87
      virtual double get_dt_max() const { // maximum time step
         assert(0 && "No maximum timestep specified!");
         abort();};
88
      virtual double get_next_plot() ; // output number for the next write
Christopher Subich's avatar
Christopher Subich committed
89
90
91
92

      /* Initialization */
      virtual double init_time() const; // Initialization time
      virtual void init_tracers(vector<DTArray *> & tracers);
93
94
95
      virtual void init_vels(DTArray & u, DTArray & v, DTArray & w) { 
         assert(0 && "init_vels not implemented");
         abort();};
96
      // Initialize velocities and tracer
97
98
99
100
      virtual void init_vels_restart(DTArray & u, DTArray & v, DTArray & w); 
      virtual void init_vels_dump(DTArray & u, DTArray & v, DTArray & w); 
      virtual void init_tracer_restart(const std::string & field, DTArray & the_tracer); 
      virtual void init_tracer_dump(const std::string & field,  DTArray & the_tracer); 
101
102
      virtual void init_grid_restart(const std::string & component,
                      const std::string & filename, DTArray & grid);
103

104
105
106
107
108
109
110
      virtual void init_tracer(int t_num, DTArray & tracer) { 
         assert(0 && "init_tracer not implemented");
         abort();}; // single-tracer

      /* Write vertical chain */
      virtual void write_chain(const char *filename, DTArray & val, int Iout, int Jout, double time);
      /* dumping functions */
111
112
113
114
      virtual void 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);
      virtual void successful_dump(int plot_number, double final_time, double plot_interval);
Christopher Subich's avatar
Christopher Subich committed
115

116
117
118
119
      virtual void write_variables(DTArray & u, DTArray & v, DTArray & w, vector<DTArray *> & tracer) { 
         assert(0 && "write_variables not defined");
         abort();}; // 

Christopher Subich's avatar
Christopher Subich committed
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
      /* Numerical checks */
      virtual double check_timestep(double step, double now);

      // Get incoming gradient operator, for differentials in analysis.  This is a null-
      // op unless the user code cares, in which case it will override this to store
      // the gradient.
      virtual void set_grad(Grad * in_grad) {return;};

      /* Forcing */
      /* Big, ultra-general forcing function */
      virtual void forcing(double t, DTArray & u, DTArray & u_f,
                                 DTArray & v, DTArray & v_f,
                                 DTArray & w, DTArray & w_f,
                                 vector<DTArray *> & tracers,
                                 vector<DTArray *> & tracers_f);
      /* No-active-tracers specialization */
      virtual void passive_forcing(double t, DTArray & u, DTArray & u_f,
                                 DTArray & v, DTArray & v_f,
                                 DTArray & w, DTArray & w_f);
      /* Independent-of-current-velocity no-tracer forcing */
      virtual void stationary_forcing(double t, DTArray& u_f, DTArray& v_f, 
            DTArray& w_f);
      
      /* If there are active tracers, split V and T focing */
      virtual void vel_forcing(double t, DTArray& u_f, DTArray& v_f,
145
146
147
                           DTArray& w_f, vector<DTArray *> & tracers) {
         assert(0 && "vel_forcing not implemented");
         abort();};
Christopher Subich's avatar
Christopher Subich committed
148
149
      virtual void tracer_forcing(double t, DTArray & u,
               DTArray & v, DTArray & w,
150
151
152
               vector<DTArray *> & tracers_f) {
         assert(0 && "tracer_forcing not implemented");
         abort();};
Christopher Subich's avatar
Christopher Subich committed
153
154
155
156
157
158
159
160

      /* Analysis and writing */

      virtual void analysis(double t, DTArray & u, DTArray & v, DTArray & w,
            vector<DTArray *> tracer, DTArray & pres); // General Analysis
      virtual void analysis(double t, DTArray & u, DTArray & v, DTArray & w,
            vector<DTArray *> tracer); // Less pressure
      virtual void vel_analysis(double t, DTArray & u, DTArray & v, 
161
162
163
164
165
            DTArray & w) {
         assert(0 && "vel_analysis not implemented");
         abort();}; // Velocity analysis
      virtual void tracer_analysis(double t, int t_num, DTArray & tracer) {
         assert(0 && "tracer_analysis not implemented");
166
         abort();}; // Single-tracer analysis
167
168
169
      template <class T> void add_diagnostic(const string str, const T val,
              string & header, string & line);
      void write_diagnostics(string header, string line, int iter, bool restarting);
170
171

      // Generate an automatic grid for unmapped cases
172
      virtual void automatic_grid(double MinX, double MinY, double MinZ, 
173
              Array<double,1> *xx=0, Array<double,1> *yy=0, Array<double,1> *zz=0);
174
175
176
177
178
179
180
181

      // Surface Stresses
      void 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);
      void 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);
Christopher Subich's avatar
Christopher Subich committed
182
183
};

184
185
186
187
188
189
190
191
#include "BaseCase_impl.cc" // Include the implementation of the add_diagnostic template

// Note explicitly instantiated add_diagnostic template functions
extern template void BaseCase::add_diagnostic<int>(const string str, const int val,
        string & header, string & line);
extern template void BaseCase::add_diagnostic<double>(const string str, const double val,
        string & header, string & line);

192
193
194
195
// 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);

Christopher Subich's avatar
Christopher Subich committed
196
197
198
extern template class FluidEvolve<BaseCase>;
typedef FluidEvolve<BaseCase> EasyFlow; // Explicit template instantiation
#endif