Commit 220ef6e5 authored by David Deepwell's avatar David Deepwell

Partial refactor of BaseCase.cpp

Update version to 2.0.2
parent 23c3647a
This diff is collapsed.
#include "../BaseCase.hpp"
void BaseCase::automatic_grid(double MinX, double MinY, double MinZ,
Array<double,1> * xx, Array<double,1> * yy, Array<double,1> * zz) {
// create arrays if needed
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());
}
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);
// Clean up
if (xxa) delete xx;
if (yya) delete yy;
if (zza) delete zz;
}
#include "../BaseCase.hpp"
/* 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, int iter,
DTArray & u, DTArray & v, DTArray & w, vector<DTArray *> & tracer) {
if (compute_time > 0) {
// default is to not dump variables
int do_dump = 0;
// check if close to end of compute time
if (master()) {
double total_run_time = clock_time - real_start_time;
double avg_clk_step = total_run_time/(iter+1); // +1 to accound for initial iter=0
if (compute_time - total_run_time < 5*(avg_write_time + avg_clk_step)) {
do_dump = 1; // true
}
}
// 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);
fprintf(stdout,"Safety dump successful\n");
}
// Die gracefully
MPI_Finalize(); exit(0);
}
}
}
#include "../BaseCase.hpp"
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);
}
}
#include "../BaseCase.hpp"
// 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);
}
}
#include "../BaseCase.hpp"
// 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);
// integral values (the first is the surface force)
double fx_tot = pssum(sum(
(*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
(*get_quad_y())(jj)*(*tx)));
double fx_abs = pssum(sum(
(*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
(*get_quad_y())(jj)*abs(*tx)));
// extrema values
double tx_max = psmax(max(*tx));
double tx_min = psmin(min(*tx));
// bottom stress ( across channel - y )
bottom_stress_y(*ty, Hprime, v, temp, gradient_op, grid_type, is_mapped(), mu);
// integral values (the first is the surface force)
double fy_tot = pssum(sum(
(*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
(*get_quad_y())(jj)*(*ty)));
double fy_abs = pssum(sum(
(*get_quad_x())(ii)*pow(1+pow(Hprime,2),0.5)*
(*get_quad_y())(jj)*abs(*ty)));
// extrema values
double ty_max = psmax(max(*ty));
double ty_min = psmin(min(*ty));
double ts_max = psmax(max(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==0 and !restarting )
fprintf(stresses_file,"Time, "
"tx_max, tx_min, ty_max, ty_min, ts_max, "
"fx_tot, fx_abs, fy_tot, fy_abs\n");
fprintf(stresses_file,"%.17f, "
"%.17g, %.17g, %.17g, %.17g, %.17g, "
"%.17g, %.17g, %.17g, %.17g\n",
time,
tx_max, tx_min, ty_max, ty_min, ts_max,
fx_tot, fx_abs, fy_tot, fy_abs);
fclose(stresses_file);
}
}
#include "../BaseCase.hpp"
// 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);
// integral values (the first is the surface force)
double fx_tot = pssum(sum(
(*get_quad_x())(ii)*
(*get_quad_y())(jj)*(*tx)));
double fx_abs = pssum(sum(
(*get_quad_x())(ii)*
(*get_quad_y())(jj)*abs(*tx)));
// extrema values
double tx_max = psmax(max(*tx));
double tx_min = psmin(min(*tx));
// top stress ( across channel - y )
top_stress_y(*ty, v, temp, gradient_op, grid_type, size_z(), mu);
// integral values (the first is the surface force)
double fy_tot = pssum(sum(
(*get_quad_x())(ii)*
(*get_quad_y())(jj)*(*ty)));
double fy_abs = pssum(sum(
(*get_quad_x())(ii)*
(*get_quad_y())(jj)*abs(*ty)));
// extrema values
double ty_max = psmax(max(*ty));
double ty_min = psmin(min(*ty));
double ts_max = psmax(max(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==0 and !restarting )
fprintf(stresses_file,"Time, "
"tx_max, tx_min, ty_max, ty_min, ts_max, "
"fx_tot, fx_abs, fy_tot, fy_abs\n");
fprintf(stresses_file,"%.17f, "
"%.17g, %.17g, %.17g, %.17g, %.17g, "
"%.17g, %.17g, %.17g, %.17g\n",
time,
tx_max, tx_min, ty_max, ty_min, ts_max,
fx_tot, fx_abs, fy_tot, fy_abs);
fclose(stresses_file);
}
}
#include "../BaseCase.hpp"
/* 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);
}
}
#include "../BaseCase.hpp"
/* 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);
}
#include "../BaseCase.hpp"
// 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 == 0 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);
}
// Explicitly instantiate common add_diagnostic templates, for integer and double
// types
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);
#include "../BaseCase.hpp"
// Write plot time information
void BaseCase::write_plot_times(double time, double clock_time, double comp_duration,
double avg_write_time, int plot_number, bool restarting) {
if (master()) {
// in log file
fprintf(stdout,"*Write time: %.6g. Average write time: %.6g.\n",
clock_time - comp_duration, avg_write_time);
// track in a file
FILE * plottimes_file = fopen("plot_times.txt","a");
assert(plottimes_file);
if ( plot_number==get_restart_sequence()+1 and !restarting )
fprintf(plottimes_file,"Output number, Simulation time (s), "
"Write time (s), Average write time (s)\n");
fprintf(plottimes_file,"%d, %.17f, %.12g, %.12g\n",
plot_number, time, clock_time - comp_duration, avg_write_time);
fclose(plottimes_file);
}
}
......@@ -60,7 +60,7 @@ all: tests/test_deriv_x tests/test_write_x tests/test_esolve_x tests/test_heat_x
.PHONY: clean
clean:
rm -f *.o tests/*.o cases/*.o cases/*.src.? tests/*.src.? Science/*.o cases/*/*.src.? cases/*/*.o
rm -f *.o tests/*.o cases/*.o cases/*.src.? tests/*.src.? Science/*.o BaseCase/*.o cases/*/*.src.? cases/*/*.o
##
## Short-names for source files
......@@ -86,16 +86,16 @@ NSIntegrator.o: NSIntegrator.cpp NSIntegrator_impl.cc
tests/test%.o: tests/test%.cpp
$(MPICXX) ${VERSION} $(CFLAGS) -o $@ -c $<
tests/test%.x: tests/test%.o tests/test%.src.o ${SCIENCE_OBJS} ${SPINS_BASE}
tests/test%.x: tests/test%.o tests/test%.src.o ${SCIENCE_OBJS} ${BASECASE_OBJS} ${SPINS_BASE}
$(LD) ${VERSION} $(LDFLAGS) -o $@ $^ $(LDLIBS)
cases/%.o: cases/%.cpp NSIntegrator_impl.cc NSIntegrator.hpp
$(MPICXX) ${VERSION} $(CFLAGS) -o $@ -c $<
cases/%.x: cases/%.o cases/%.src.o ${SCIENCE_OBJS} ${SPINS_BASE}
cases/%.x: cases/%.o cases/%.src.o ${SCIENCE_OBJS} ${BASECASE_OBJS} ${SPINS_BASE}
$(LD) ${VERSION} $(LDFLAGS) -o $@ $^ $(LDLIBS)
cases/%_x: cases/%.o cases/%.src.o ${SCIENCE_OBJS} ${SPINS_BASE}
cases/%_x: cases/%.o cases/%.src.o ${SCIENCE_OBJS} ${BASECASE_OBJS} ${SPINS_BASE}
$(LD) ${VERSION} $(LDFLAGS) -o $@ $^ $(LDLIBS)
tests/test%.src.c: tests/test%.cpp CaseFileSource.c
......
MAJOR_VERSION=2
MINOR_VERSION=0
PATCH_VERSION=1
PATCH_VERSION=2
Markdown is supported
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