Commit d4bb5300 authored by Christopher Subich's avatar Christopher Subich
Browse files

Merge remote-tracking branch 'david/dump'

parents 355ad2b3 c6a6c09a
......@@ -10,203 +10,341 @@ 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=0, Array<double,1> * yy=0, Array<double,1> * zz = 0){
//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;
}
/* 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);
}
/* 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);
}
}
......@@ -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,28 @@ class BaseCase {
virtual void init_vels(DTArray & u, DTArray & v, DTArray & w) {
assert(0 && "init_vels not implemented");
abort();};
// Initialize velocities and tracer
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 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 */
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 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 +159,8 @@ 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,
Array<double,1> *xx, Array<double,1> *yy, Array<double,1> *zz);
};
extern template class FluidEvolve<BaseCase>;
......
/* Generic script for two layer fluid (tanh profile)
with zero initial velocity
and no topography */
// Required headers
#include "../TArray.hpp" // Custom extensions to the library to support FFTs
#include "../NSIntegrator.hpp" // Time-integrator for the Navier-Stokes equations
#include "../BaseCase.hpp" // Support file that contains default implementations of several functions
#include "../Options.hpp" // config-file parser
#include <random/normal.h> // Blitz random number generator
#include <blitz/array.h> // Blitz++ array library
#include <mpi.h> // MPI parallel library
#include <iostream>
#include <fstream>
#include <stdlib.h>
//#include "../Science.hpp" // Additional analysis routines
using namespace std;
using namespace NSIntegrator;
using namespace ranlib;
// Tensor variables for indexing
blitz::firstIndex ii;
blitz::secondIndex jj;
blitz::thirdIndex kk;
/* ------------------ Parameters --------------------- */
// Grid scales
double Lx, Ly, Lz; // (m)
int Nx, Ny, Nz; // Points in x, y (span), z directions
double MinX, MinY, MinZ; // Minimum x/y/z points
// Grid types
DIMTYPE intype_x, intype_y, intype_z;
// Physical constants
double g, rot_f; // gravity accel (m/s^2), Coriolis frequency (s^-1)
// Stratification parameters
double rho_0; // reference density (kg/L)
double delta_rho; // density difference between top and bottom layers (kg/L)
double pyc_asym; // % of depth to shift pycnocline above the mid-depth
double h_perc; // pycnocline half-width as % of depth
double h_mix_perc; // vertical half-width transition of mixed region
// Horizontal stratification parameters
double delta_x; // horizontal transition length (m)
double Lmix; // Width of mixed region (m)
double Hmix; // Total height of mixed region (m)
// Viscosity and diffusivity of density and tracers
double VISCO;
double DIFFU_rho;
double DIFFU_dye_1;
// Temporal parameters
double plot_interval; // Time between field writes (s)
double final_time; // Final time (s)
// Vertical chain parameters
bool savechain; // (boolean) Flag to use chain or not
double chain1_start_time; // time to start saving chain (s)
double chain1_end_time; // time to stop saving chain (s)
double chain_plot_interval; // time between chain writes (s)
// Chain locations (defined in main program)
int chain1_xindex, chain1_yindex;
// Initial velocity perturbation
double perturb;
// Dump parameters
bool restart_from_dump;
double real_start_time;
double compute_time;
double total_run_time;
double avg_write_time;
// Restarting options
bool restarting = false;
double restart_time;
double initial_time = 0;
int restart_sequence;
// Iteration counter
int itercount = 0;
/* ------------------ Derived parameters --------------------- */
// Pycnocline half-width
double h_halfwidth;
double h_mix_half;
// Diffusivity information
const int N_tracers = 2; // must be self defined, and add one for density
double DIFFU[ N_tracers ];
double *DIFFU_pointer;
// Flow speed
double c0;
// Squared maximum buoyancy frequency if the initial stratification was stable
double N2_max;