Friendlier aborts in BaseCase.hpp

src/BaseCase.hpp

@@ -28,12 +28,16 @@ class BaseCase {
       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
-      virtual int size_cube() const { abort(); }; // Special case -- N*N*N
+      virtual int size_cube() const { 
+         assert(0 && "size_cube not implemented");
+         abort(); }; // Special case -- N*N*N
 
       virtual double length_x() const; // Length in x
       virtual double length_y() const; // Length in y
       virtual double length_z() const; // Length in z
-      virtual double length_cube() const {abort();}; // Special case -- L*L*L
+      virtual double length_cube() const {
+         assert(0 && "length_cube not implemented");
+         abort();}; // Special case -- L*L*L
 
       virtual DIMTYPE type_x() const; // Expansion type in x
       virtual DIMTYPE type_y() const; // Expansion type in y
@@ -76,9 +80,13 @@ class BaseCase {
       /* Initialization */
       virtual double init_time() const; // Initialization time
       virtual void init_tracers(vector<DTArray *> & tracers);
-      virtual void init_vels(DTArray & u, DTArray & v, DTArray & w) { abort();};
+      virtual void init_vels(DTArray & u, DTArray & v, DTArray & w) { 
+         assert(0 && "init_vels not implemented");
+         abort();};
 
-      virtual void init_tracer(int t_num, DTArray & tracer) { abort();}; // single-tracer
+      virtual void init_tracer(int t_num, DTArray & tracer) { 
+         assert(0 && "init_tracer not implemented");
+         abort();}; // single-tracer
 
       /* Numerical checks */
       virtual double check_timestep(double step, double now);
@@ -105,10 +113,14 @@ class BaseCase {
       
       /* If there are active tracers, split V and T focing */
       virtual void vel_forcing(double t, DTArray& u_f, DTArray& v_f,
-                           DTArray& w_f, vector<DTArray *> & tracers) {abort();};
+                           DTArray& w_f, vector<DTArray *> & tracers) {
+         assert(0 && "vel_forcing not implemented");
+         abort();};
       virtual void tracer_forcing(double t, DTArray & u,
                DTArray & v, DTArray & w,
-               vector<DTArray *> & tracers_f) {abort();};
+               vector<DTArray *> & tracers_f) {
+         assert(0 && "tracer_forcing not implemented");
+         abort();};
 
       /* Analysis and writing */
 
@@ -117,8 +129,12 @@ class BaseCase {
       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, 
-            DTArray & w) {abort();}; // Velocity analysis
-      virtual void tracer_analysis(double t, int t_num, DTArray & tracer) {abort();};
+            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");
+         abort();};
          // Single-tracer analysis
 };
 

src/cases/bench_dipole.cpp

@@ -260,9 +260,9 @@ class userControl : public BaseCase {
 //         write_array(grid,"zgrid"); write_reader(grid,"zgrid",false);
       }
       // Once initialized, this is freely-evolving flow.  No forcing is necessary
-      void passive_forcing(double t, const DTArray & u, DTArray & u_f, 
-            const DTArray & v, DTArray & v_f, 
-            const DTArray & w, DTArray & w_f) {
+      void passive_forcing(double t, DTArray & u, DTArray & u_f, 
+            DTArray & v, DTArray & v_f, 
+            DTArray & w, DTArray & w_f) {
          u_f = 0;
          v_f = 0;
          w_f = 0;

src/cases/bench_kh_3d.cpp

@@ -143,8 +143,8 @@ class helmholtz : public BaseCase {
          w_f = -g*((*tracers[0]))/rho_0;
       }
       // Forcing of the perturbation temperature
-      void tracer_forcing(double t, const DTArray & u, const DTArray & v,
-            const DTArray & w, vector<DTArray *> & tracers_f) {
+      void tracer_forcing(double t, DTArray & u, DTArray & v,
+            DTArray & w, vector<DTArray *> & tracers_f) {
          /* Since the perturbation temperature is a perturbation, its forcing is
             proportional to the background temperature gradient and the w-velocity */
          *tracers_f[0] = w(ii,jj,kk)*0.5*delta_rho/dz_rho*pow(cosh((zz(kk)-0.5)/dz_rho),-2);

src/cases/bench_kh_big.cpp

@@ -143,8 +143,8 @@ class helmholtz : public BaseCase {
          w_f = -g*((*tracers[0]))/rho_0;
       }
       // Forcing of the perturbation temperature
-      void tracer_forcing(double t, const DTArray & u, const DTArray & v,
-            const DTArray & w, vector<DTArray *> & tracers_f) {
+      void tracer_forcing(double t, DTArray & u, DTArray & v,
+            DTArray & w, vector<DTArray *> & tracers_f) {
          /* Since the perturbation temperature is a perturbation, its forcing is
             proportional to the background temperature gradient and the w-velocity */
          *tracers_f[0] = w(ii,jj,kk)*0.5*delta_rho/dz_rho*pow(cosh((zz(kk)-0.5)/dz_rho),-2);

src/cases/bench_kh_billow.cpp

@@ -142,8 +142,8 @@ class helmholtz : public BaseCase {
          w_f = -g*((*tracers[0]))/rho_0;
       }
       // Forcing of the perturbation temperature
-      void tracer_forcing(double t, const DTArray & u, const DTArray & v,
-            const DTArray & w, vector<DTArray *> & tracers_f) {
+      void tracer_forcing(double t, DTArray & u, DTArray & v,
+            DTArray & w, vector<DTArray *> & tracers_f) {
          /* Since the perturbation temperature is a perturbation, its forcing is
             proportional to the background temperature gradient and the w-velocity */
          *tracers_f[0] = w(ii,jj,kk)*0.5*delta_rho/dz_rho*pow(cosh((zz(kk)-0.5)/dz_rho),-2);

src/cases/dipole_2d.cpp

@@ -185,9 +185,9 @@ class userControl : public BaseCase {
          grid = zz(kk);
          write_array(grid,"zgrid"); write_reader(grid,"zgrid",false);
       }
-      void passive_forcing(double t, const DTArray & u, DTArray & u_f, 
-            const DTArray & v, DTArray & v_f, 
-            const DTArray & w, DTArray & w_f) {
+      void passive_forcing(double t, DTArray & u, DTArray & u_f, 
+            DTArray & v, DTArray & v_f, 
+            DTArray & w, DTArray & w_f) {
          u_f = 0;
          v_f = 0;
          w_f = 0;

src/cases/dipole_2d_25k.cpp

@@ -185,9 +185,9 @@ class userControl : public BaseCase {
          grid = zz(kk);
          write_array(grid,"zgrid"); write_reader(grid,"zgrid",false);
       }
-      void passive_forcing(double t, const DTArray & u, DTArray & u_f, 
-            const DTArray & v, DTArray & v_f, 
-            const DTArray & w, DTArray & w_f) {
+      void passive_forcing(double t, DTArray & u, DTArray & u_f, 
+            DTArray & v, DTArray & v_f, 
+            DTArray & w, DTArray & w_f) {
          u_f = 0;
          v_f = 0;
          w_f = 0;

src/cases/dipole_3d.cpp

@@ -214,9 +214,9 @@ class userControl : public BaseCase {
             write_array(grid,"zgrid"); write_reader(grid,"zgrid",false);
          }
       }
-      void passive_forcing(double t, const DTArray & u, DTArray & u_f, 
-            const DTArray & v, DTArray & v_f, 
-            const DTArray & w, DTArray & w_f) {
+      void passive_forcing(double t, DTArray & u, DTArray & u_f, 
+            DTArray & v, DTArray & v_f, 
+            DTArray & w, DTArray & w_f) {
          u_f = 0;
          v_f = 0;
          w_f = 0;

src/cases/doc_benard.cpp

@@ -156,8 +156,8 @@ class benard : public BaseCase {
          w_f = g*(alpha*(*tracers[0]))/rho_0;
       }
       // Forcing of the perturbation temperature
-      void tracer_forcing(double t, const DTArray & u, const DTArray & v,
-            const DTArray & w, vector<DTArray *> & tracers_f) {
+      void tracer_forcing(double t, DTArray & u, DTArray & v,
+            DTArray & w, vector<DTArray *> & tracers_f) {
          /* Since the perturbation temperature is a perturbation, its forcing is
             proportional to the background temperature gradient and the w-velocity */
          *tracers_f[0] = w*(DELTA_T/LENGTH_Z);

src/cases/doc_dambreak.cpp

@@ -176,8 +176,8 @@ class dambreak : public BaseCase {
          w_f = -g*((*tracers[0]))/rho_0;
       }
       // Forcing of the density (zero)
-      void tracer_forcing(double t, const DTArray & u, const DTArray & v,
-            const DTArray & w, vector<DTArray *> & tracers_f) {
+      void tracer_forcing(double t, DTArray & u, DTArray & v,
+            DTArray & w, vector<DTArray *> & tracers_f) {
          *tracers_f[0] = 0;
       }
 

src/cases/doc_dipole.cpp

@@ -239,9 +239,9 @@ class userControl : public BaseCase {
          write_array(grid,"zgrid"); write_reader(grid,"zgrid",false);
       }
       // Once initialized, this is freely-evolving flow.  No forcing is necessary
-      void passive_forcing(double t, const DTArray & u, DTArray & u_f, 
-            const DTArray & v, DTArray & v_f, 
-            const DTArray & w, DTArray & w_f) {
+      void passive_forcing(double t, DTArray & u, DTArray & u_f, 
+            DTArray & v, DTArray & v_f, 
+            DTArray & w, DTArray & w_f) {
          u_f = 0;
          v_f = 0;
          w_f = 0;

src/cases/doc_dipole_hv.cpp

@@ -239,9 +239,9 @@ class userControl : public BaseCase {
          write_array(grid,"zgrid"); write_reader(grid,"zgrid",false);
       }
       // Once initialized, this is freely-evolving flow.  No forcing is necessary
-      void passive_forcing(double t, const DTArray & u, DTArray & u_f, 
-            const DTArray & v, DTArray & v_f, 
-            const DTArray & w, DTArray & w_f) {
+      void passive_forcing(double t, DTArray & u, DTArray & u_f, 
+            DTArray & v, DTArray & v_f, 
+            DTArray & w, DTArray & w_f) {
          u_f = 0;
          v_f = 0;
          w_f = 0;

src/cases/doc_kh_billow.cpp

@@ -133,8 +133,8 @@ class helmholtz : public BaseCase {
          w_f = -g*((*tracers[0]))/rho_0;
       }
       // Forcing of the perturbation temperature
-      void tracer_forcing(double t, const DTArray & u, const DTArray & v,
-            const DTArray & w, vector<DTArray *> & tracers_f) {
+      void tracer_forcing(double t, DTArray & u, DTArray & v,
+            DTArray & w, vector<DTArray *> & tracers_f) {
          /* Since the perturbation temperature is a perturbation, its forcing is
             proportional to the background temperature gradient and the w-velocity */
          *tracers_f[0] = w(ii,jj,kk)*0.5*delta_rho/dz_rho*pow(cosh((zz(kk)-0.5)/dz_rho),-2);

src/cases/doc_kh_filter.cpp

@@ -134,8 +134,8 @@ class helmholtz : public BaseCase {
          w_f = -g*((*tracers[0]))/rho_0;
       }
       // Forcing of the perturbation temperature
-      void tracer_forcing(double t, const DTArray & u, const DTArray & v,
-            const DTArray & w, vector<DTArray *> & tracers_f) {
+      void tracer_forcing(double t, DTArray & u, DTArray & v,
+            DTArray & w, vector<DTArray *> & tracers_f) {
          /* Since the perturbation temperature is a perturbation, its forcing is
             proportional to the background temperature gradient and the w-velocity */
          *tracers_f[0] = w(ii,jj,kk)*0.5*delta_rho/dz_rho*pow(cosh((zz(kk)-0.5)/dz_rho),-2);

src/cases/doc_kh_hv.cpp

@@ -134,8 +134,8 @@ class helmholtz : public BaseCase {
          w_f = -g*((*tracers[0]))/rho_0;
       }
       // Forcing of the perturbation temperature
-      void tracer_forcing(double t, const DTArray & u, const DTArray & v,
-            const DTArray & w, vector<DTArray *> & tracers_f) {
+      void tracer_forcing(double t, DTArray & u, DTArray & v,
+            DTArray & w, vector<DTArray *> & tracers_f) {
          /* Since the perturbation temperature is a perturbation, its forcing is
             proportional to the background temperature gradient and the w-velocity */
          *tracers_f[0] = w(ii,jj,kk)*0.5*delta_rho/dz_rho*pow(cosh((zz(kk)-0.5)/dz_rho),-2);

src/cases/doc_map_const.cpp

@@ -149,8 +149,8 @@ class mapiw : public BaseCase {
 
       // Forcing must be done generally, since both rotation and density are
       // involved
-      void forcing(double t, const DTArray & u, DTArray & u_f, 
-            const DTArray & v, DTArray & v_f, const DTArray & w,
+      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) {
          // Rotation couples u and v, plus a source term for the tide

src/cases/doc_map_iwave.cpp

@@ -149,8 +149,8 @@ class mapiw : public BaseCase {
 
       // Forcing must be done generally, since both rotation and density are
       // involved
-      void forcing(double t, const DTArray & u, DTArray & u_f, 
-            const DTArray & v, DTArray & v_f, const DTArray & w,
+      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) {
          // Rotation couples u and v, plus a source term for the tide

src/cases/mapseich.cpp

@@ -176,8 +176,8 @@ class userControl : public BaseCase {
          u_f = 0; v_f = 0;
          w_f = -g*(rho_0 + *tracers[0]);
       }
-      void tracer_forcing(double t, const DTArray & u, const DTArray & v,
-            const DTArray & w, vector<DTArray *> & tracers_f) {
+      void tracer_forcing(double t, DTArray & u, DTArray & v,
+            DTArray & w, vector<DTArray *> & tracers_f) {
          /* Forcing on the tracers themselves.  Since rho is a passive density,
             there is none. */
          *tracers_f[0] = 0;

src/cases/mpemba.cpp

@@ -182,8 +182,8 @@ class userControl : public BaseCase {
       
       // Forcing must be done generally, since both rotation and density are
       // involved
-      void forcing(double t, const DTArray & u, DTArray & u_f, 
-            const DTArray & v, DTArray & v_f, const DTArray & w,
+      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) {
          u_f = g*(*(tracers[0]))/rho0;

src/cases/mpemba_spinup.cpp

@@ -180,8 +180,8 @@ class userControl : public BaseCase {
       
       // Forcing must be done generally, since both rotation and density are
       // involved
-      void forcing(double t, const DTArray & u, DTArray & u_f, 
-            const DTArray & v, DTArray & v_f, const DTArray & w,
+      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) {
          u_f = g*(*(tracers[0]))/rho0;