Skip to content
Snippets Groups Projects
Commit c71a43f0 authored by David Deepwell's avatar David Deepwell
Browse files

Added safety dump to BaseCase and tankmode2.cpp 

Safety dump estimates output write time and writes fields with enough time before the computation limit is reached
parent 0ae6d2aa
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
src/BaseCase.cpp
+ 244
− 121
View file @ c71a43f0
......@@ -10,203 +10,326 @@ using std::vector;
/* Call the source code writing function in the constructor */
extern "C" {
void WriteCaseFileSource(void);
void WriteCaseFileSource(void);
}
BaseCase::BaseCase(void)
{
if (master()) WriteCaseFileSource();
if (master()) WriteCaseFileSource();
}
/* Implementation of non-abstract methods in BaseCase */
int BaseCase::numActive() const { return 0; }
int BaseCase::numPassive() const { return 0; }
int BaseCase::numtracers() const { /* total number of tracers */
return numActive() + numPassive();
return numActive() + numPassive();
}
int BaseCase::size_x() const {
return size_cube();
return size_cube();
}
int BaseCase::size_y() const {
return size_cube();
return size_cube();
}
int BaseCase::size_z() const {
return size_cube();
return size_cube();
}
double BaseCase::length_x() const {
return length_cube();
return length_cube();
}
double BaseCase::length_y() const {
return length_cube();
return length_cube();
}
double BaseCase::length_z() const {
return length_cube();
return length_cube();
}
DIMTYPE BaseCase::type_x() const {
return type_default();
return type_default();
}
DIMTYPE BaseCase::type_y() const {
return type_default();
return type_default();
}
DIMTYPE BaseCase::type_z() const {
return type_default();
return type_default();
}
DIMTYPE BaseCase::type_default() const {
return PERIODIC;
return PERIODIC;
}
void BaseCase::tracer_bc_x(int t_num, double & dir, double & neu) const {
if (!zero_tracer_boundary) {
dir = 0;
neu = 1;
}
else {
dir = 1;
neu = 0;
}
return;
if (!zero_tracer_boundary) {
dir = 0;
neu = 1;
}
else {
dir = 1;
neu = 0;
}
return;
}
void BaseCase::tracer_bc_y(int t_num, double & dir, double & neu) const {
if (!zero_tracer_boundary) {
dir = 0;
neu = 1;
}
else {
dir = 1;
neu = 0;
}
return;
if (!zero_tracer_boundary) {
dir = 0;
neu = 1;
}
else {
dir = 1;
neu = 0;
}
return;
}
void BaseCase::tracer_bc_z(int t_num, double & dir, double & neu) const {
if (!zero_tracer_boundary) {
dir = 0;
neu = 1;
}
else {
dir = 1;
neu = 0;
}
return;
if (!zero_tracer_boundary) {
dir = 0;
neu = 1;
}
else {
dir = 1;
neu = 0;
}
return;
}
bool BaseCase::tracer_bc_forcing() const {
return false;
return false;
}
bool BaseCase::is_mapped() const { // Whether this problem has mapped coordinates
return false;
return false;
}
// 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) {
return;
return;
}
/* Physical parameters */
double BaseCase::get_visco() const {
return 0;
return 0;
}
double BaseCase::get_diffusivity(int t) const {
return 0;
return 0;
}
double BaseCase::get_rot_f() const {
return 0;
}
int BaseCase::get_restart_sequence() const {
return 0;
}
/* Initialization */
double BaseCase::init_time() const {
return 0;
return 0;
}
void BaseCase::init_tracers(vector<DTArray *> & tracers) {
/* 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]));
}
/* 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]));
}
}
double BaseCase::check_timestep(double step, double now) {
return step; // Default, accept the given stepsize
return step; // Default, accept the given stepsize
}
/* Forcing */
void BaseCase::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) {
/* 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;
}
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;
}
}
void BaseCase::passive_forcing(double t, DTArray & u, DTArray & u_f,
DTArray & v, DTArray & v_f, DTArray &, DTArray & w_f) {
/* Reduce to velocity-independent case */
stationary_forcing(t, u_f, v_f, w_f);
DTArray & v, DTArray & v_f, DTArray &, DTArray & w_f) {
/* Reduce to velocity-independent case */
stationary_forcing(t, u_f, v_f, w_f);
}
void BaseCase::stationary_forcing(double t, DTArray & u_f, DTArray & v_f,
DTArray & w_f) {
/* Default case, 0 forcing */
u_f = 0;
v_f = 0;
w_f = 0;
DTArray & w_f) {
/* Default case, 0 forcing */
u_f = 0;
v_f = 0;
w_f = 0;
}
/* Analysis */
void BaseCase::analysis(double t, DTArray & u, DTArray & v, DTArray & w,
vector<DTArray *> tracer, DTArray & pres) {
analysis(t,u,v,w,tracer);
vector<DTArray *> tracer, DTArray & pres) {
analysis(t,u,v,w,tracer);
}
void BaseCase::analysis(double t, DTArray & u, DTArray & v, DTArray & w,
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]));
}
}
void BaseCase::automatic_grid(double MinX,double MinY,double MinZ) {
Array<double,1> xx(split_range(size_x())), yy(size_y()), zz(size_z());
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) {
xx = MinX+length_x()*(0.5+0.5*cos(M_PI*ii/(size_x()-1)));
} else {
xx = MinX + length_x()*(ii+0.5)/size_x();
}
yy = MinY + length_y()*(ii+0.5)/size_y();
if (type_z() == NO_SLIP) {
zz = MinZ+length_z()*(0.5+0.5*cos(M_PI*ii/(size_z()-1)));
} else {
zz = MinZ + length_z()*(0.5+ii)/size_z();
}
// Write grid/reader
grid = xx(ii) + 0*jj + 0*kk;
write_array(grid,"xgrid");
write_reader(grid,"xgrid",false);
if (size_y() > 1) {
grid = 0*ii + yy(jj) + 0*kk;
write_array(grid,"ygrid");
write_reader(grid,"ygrid",false);
}
grid = 0*ii + 0*jj + zz(kk);
write_array(grid,"zgrid");
write_reader(grid,"zgrid",false);
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]));
}
}
void BaseCase::automatic_grid(double MinX, double MinY, double MinZ){
Array<double,1> xx(split_range(size_x())), yy(size_y()), zz(size_z());
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) {
xx = MinX+length_x()*(0.5+0.5*cos(M_PI*ii/(size_x()-1)));
} else {
xx = MinX + length_x()*(ii+0.5)/size_x();
}
yy = MinY + length_y()*(ii+0.5)/size_y();
if (type_z() == NO_SLIP) {
zz = MinZ+length_z()*(0.5+0.5*cos(M_PI*ii/(size_z()-1)));
} else {
zz = MinZ + length_z()*(0.5+ii)/size_z();
}
// Write grid/reader
grid = xx(ii) + 0*jj + 0*kk;
write_array(grid,"xgrid");
write_reader(grid,"xgrid",false);
if (size_y() > 1) {
grid = 0*ii + yy(jj) + 0*kk;
write_array(grid,"ygrid");
write_reader(grid,"ygrid",false);
}
grid = 0*ii + 0*jj + zz(kk);
write_array(grid,"zgrid");
write_reader(grid,"zgrid",false);
}
/* 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;
}
/* 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){
int do_dump = 0;
if (master()) {
double total_run_time = clock_time - real_start_time;
// check if close to end of compute time
if ((compute_time > 0) && (compute_time - total_run_time < 3*avg_write_time)){
do_dump = 1; // true
}
}
MPI_Bcast(&do_dump, 1, MPI_INT, 0, MPI_COMM_WORLD);
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%.12g\n", sim_time);
fprintf(dump_file,"The dump index was:\n%d\n", plot_number);
fclose(dump_file);
}
// Die gracefully
MPI_Finalize(); exit(0);
}
}
/* 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);
}
}
src/BaseCase.hpp
+ 18
− 1
View file @ c71a43f0
......@@ -79,6 +79,8 @@ class BaseCase {
/* Physical parameters */
virtual double get_visco() const; // Physical viscosity
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
/* Initialization */
virtual double init_time() const; // Initialization time
......@@ -86,11 +88,26 @@ class BaseCase {
virtual void init_vels(DTArray & u, DTArray & v, DTArray & w) {
assert(0 && "init_vels not implemented");
abort();};
// Initialize velocities if restarting from regular output
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);
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);
virtual void init_tracer(int t_num, DTArray & tracer) {
assert(0 && "init_tracer not implemented");
abort();}; // single-tracer
virtual void write_variables(DTArray & u, DTArray & v, DTArray & w, vector<DTArray *> & tracer) {
assert(0 && "write_variables not defined");
abort();}; //
/* Numerical checks */
virtual double check_timestep(double step, double now);
......@@ -140,7 +157,7 @@ class BaseCase {
abort();}; // Single-tracer analysis
// Generate an automatic grid for unmapped cases
virtual void automatic_grid(double MinX,double MinY,double MinZ);
virtual void automatic_grid(double MinX, double MinY, double MinZ);
};
extern template class FluidEvolve<BaseCase>;
......
src/cases/tankmode2.cpp 0 → 100644
+ 618
− 0
View file @ c71a43f0
This diff is collapsed.
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