Commit 428b5399 authored by David Deepwell's avatar David Deepwell
Browse files

Move time step checker to BaseCase

The empty shell of check_timestep in BaseCase has been filled
with what was essentially the same thing in each case file.
Again, this is just more case file house cleaning.

A few ancillary functions have also been written to help.
The most important might be get_dt_max() which returns dt_max.
dt_max can either be specified in spins.conf or defined in the
case file based on the buoyancy frequency. Currently, if
dt_max > 0 and in spins.conf, then that value is used. Otherwise,
0.5/sqrt(N2_max) is used.

The re-definition of check_timestep in each case file can now be
removed so long as the associated ancillary functions and
parameters are included.
parent 25b5d62f
......@@ -103,19 +103,12 @@ void BaseCase::do_mapping(DTArray & xgrid, DTArray & ygrid, DTArray & zgrid) {
}
/* Physical parameters */
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_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_dt_max() const { return 0; }
double BaseCase::get_next_plot() { return 0; }
/* Initialization */
double BaseCase::init_time() const {
......@@ -130,8 +123,29 @@ void BaseCase::init_tracers(vector<DTArray *> & tracers) {
}
}
double BaseCase::check_timestep(double step, double now) {
return step; // Default, accept the given stepsize
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);
}
}
/* Forcing */
......
......@@ -81,6 +81,8 @@ class BaseCase {
virtual double get_diffusivity(int tracernum) const; // Diffusivity
virtual double get_rot_f() const; // Physical rotation rate
virtual int get_restart_sequence() const; // restart sequence
virtual double get_dt_max() const; // maximum time step
virtual double get_next_plot() ; // output number for the next write
/* Initialization */
virtual double init_time() const; // Initialization time
......
......@@ -42,6 +42,7 @@ double Lmix; // Width of mixed region (m)
// Temporal parameters
double final_time; // Final time (s)
double plot_interval; // Time between field writes (s)
double dt_max; // maximum time step (s)
// Restarting options
bool restarting; // are you restarting?
......@@ -95,6 +96,7 @@ class userControl : public BaseCase {
DIMTYPE type_z() const { return intype_z; }
// Coriolis parameter, viscosity, and diffusivities
double get_rot_f() const { return rot_f; }
double get_visco() const { return visco; }
double get_diffusivity(int t_num) const {
return kappa_rho;
......@@ -103,33 +105,12 @@ class userControl : public BaseCase {
// Temporal parameters
double init_time() const { return initial_time; }
int get_restart_sequence() const { return restart_sequence; }
double get_dt_max() const { return dt_max; }
double get_next_plot() { return next_plot; }
// Number of tracers (the first is density)
int numtracers() const { return Num_tracers; }
/* 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 = next_plot - 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 */
void init_vels(DTArray & u, DTArray & v, DTArray & w) {
if (master()) fprintf(stdout,"Initializing velocities\n");
......@@ -384,6 +365,7 @@ int main(int argc, char ** argv) {
option_category("Temporal options");
add_option("final_time",&final_time,"Final time");
add_option("plot_interval",&plot_interval,"Time between writes");
add_option("dt_max",&dt_max,"Maximum time step");
option_category("Restart options");
add_option("restart",&restarting,false,"Restart from prior output time.");
......@@ -435,6 +417,11 @@ int main(int argc, char ** argv) {
mu = visco*rho_0;
// Maximum buoyancy frequency (squared) if the initial stratification was stable
N2_max = g*delta_rho/(2*delta_x);
// Maximum time step
if (dt_max <= 0) {
// if dt_max not given in spins.conf, use the buoyancy frequency
dt_max = 0.5/sqrt(N2_max);
}
/* ------------------ Print some parameters --------------------- */
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment