Commit 2e379eb9 authored by Ben Storer's avatar Ben Storer
Browse files

Adding a bunch of stuff

- Added 1D and 2D outputs
- Added parallel writing and reading
parent 1ac02e05
......@@ -180,7 +180,7 @@ void BaseCase::analysis(double t, DTArray & u, DTArray & v, DTArray & w,
void BaseCase::automatic_grid(double MinX, double MinY, double MinZ,
Array<double,1> * xx, Array<double,1> * yy, Array<double,1> * zz){
//Array<double,1> xx(split_range(size_x())), yy(size_y()), zz(size_z());
bool xxa = false, yya = false, zza = false;
if (!xx) {
xxa = true; // Delete xx when returning
......@@ -216,17 +216,17 @@ void BaseCase::automatic_grid(double MinX, double MinY, double MinZ,
// Write grid/reader
grid = (*xx)(ii) + 0*jj + 0*kk;
write_array(grid,"xgrid");
write_array_old(grid,"xgrid");
write_reader(grid,"xgrid",false);
if (size_y() > 1) {
grid = 0*ii + (*yy)(jj) + 0*kk;
write_array(grid,"ygrid");
write_array_old(grid,"ygrid");
write_reader(grid,"ygrid",false);
}
grid = 0*ii + 0*jj + (*zz)(kk);
write_array(grid,"zgrid");
write_array_old(grid,"zgrid");
write_reader(grid,"zgrid",false);
// Clean up
......@@ -304,20 +304,6 @@ void BaseCase::init_tracer_dump(const std::string & field, DTArray & the_tracer)
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,
......
......@@ -98,8 +98,6 @@ class BaseCase {
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,
......
......@@ -9,6 +9,7 @@
#include "Split_reader.hpp"
#include "T_util.hpp"
#include "Parformer.hpp"
#include <stdlib.h>
// Marek's Overturning Diagnostic
......@@ -273,29 +274,285 @@ void vorticity(TArrayn::DTArray & u, TArrayn::DTArray & v,
return;
}
void dudx(TArrayn::DTArray & u, TArrayn::DTArray * & out,
double Lx, double Ly, double Lz,
int szx, int szy, int szz,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z) {
static Transformer::Trans1D trans_x(szx,szy,szz,firstDim,
(DIM_X == PERIODIC ? FOURIER : REAL)),
trans_y(szx,szy,szz,secondDim,
(DIM_Y == PERIODIC ? FOURIER : REAL)),
trans_z (szx,szy,szz,thirdDim,
(DIM_Z == PERIODIC ? FOURIER : REAL));
static blitz::TinyVector<int,3>
local_lbounds(alloc_lbound(szx,szy,szz)),
local_extent(alloc_extent(szx,szy,szz));
static blitz::GeneralArrayStorage<3>
local_storage(alloc_storage(szx,szy,szz));
static DTArray temp(local_lbounds,local_extent,local_storage);
if (DIM_X == PERIODIC) {
deriv_fft(u,trans_x,temp);
temp *= 2*M_PI/Lx;
} else if (DIM_X == FREE_SLIP) {
deriv_dst(u,trans_x,temp);
temp *= M_PI/Lx;
} else {
assert(DIM_X == NO_SLIP);
deriv_cheb(u,trans_x,temp);
temp *= -2/Lx;
}
out = &temp;
return;
}
void dudy(TArrayn::DTArray & u, TArrayn::DTArray * & out,
double Lx, double Ly, double Lz,
int szx, int szy, int szz,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z) {
static Transformer::Trans1D trans_x(szx,szy,szz,firstDim,
(DIM_X == PERIODIC ? FOURIER : REAL)),
trans_y(szx,szy,szz,secondDim,
(DIM_Y == PERIODIC ? FOURIER : REAL)),
trans_z (szx,szy,szz,thirdDim,
(DIM_Z == PERIODIC ? FOURIER : REAL));
static blitz::TinyVector<int,3>
local_lbounds(alloc_lbound(szx,szy,szz)),
local_extent(alloc_extent(szx,szy,szz));
static blitz::GeneralArrayStorage<3>
local_storage(alloc_storage(szx,szy,szz));
static DTArray temp(local_lbounds,local_extent,local_storage);
if (DIM_Y == PERIODIC) {
deriv_fft(u,trans_y,temp);
temp *= 2*M_PI/Ly;
} else if (DIM_Y == FREE_SLIP) {
deriv_dct(u,trans_y,temp);
temp *= M_PI/Ly;
} else {
assert(DIM_Y == NO_SLIP);
deriv_cheb(u,trans_y,temp);
temp *= -2/Ly;
}
out = &temp;
return;
}
void dudz(TArrayn::DTArray & u, TArrayn::DTArray * & out,
double Lx, double Ly, double Lz,
int szx, int szy, int szz,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z) {
static Transformer::Trans1D trans_x(szx,szy,szz,firstDim,
(DIM_X == PERIODIC ? FOURIER : REAL)),
trans_y(szx,szy,szz,secondDim,
(DIM_Y == PERIODIC ? FOURIER : REAL)),
trans_z (szx,szy,szz,thirdDim,
(DIM_Z == PERIODIC ? FOURIER : REAL));
static blitz::TinyVector<int,3>
local_lbounds(alloc_lbound(szx,szy,szz)),
local_extent(alloc_extent(szx,szy,szz));
static blitz::GeneralArrayStorage<3>
local_storage(alloc_storage(szx,szy,szz));
static DTArray temp(local_lbounds,local_extent,local_storage);
if (DIM_Z == PERIODIC) {
deriv_fft(u,trans_z,temp);
temp *= 2*M_PI/Lz;
} else if (DIM_Z == FREE_SLIP) {
deriv_dct(u,trans_z,temp);
temp *= M_PI/Lz;
} else {
assert(DIM_Z == NO_SLIP);
deriv_cheb(u,trans_z,temp);
temp *= -2/Lz;
}
out = &temp;
return;
}
void ertel_pv(TArrayn::DTArray & u, TArrayn::DTArray & v,
TArrayn::DTArray & w, TArrayn::DTArray & rho,
TArrayn::DTArray * & e_pv, double f0, double N0,
double Lx, double Ly, double Lz,
int szx, int szy, int szz,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z) {
static int Nx = 0, Ny = 0, Nz = 0;
static Transformer::Trans1D trans_x(szx,szy,szz,firstDim,
(DIM_X == PERIODIC ? FOURIER : REAL)),
trans_y(szx,szy,szz,secondDim,
(DIM_Y == PERIODIC ? FOURIER : REAL)),
trans_z (szx,szy,szz,thirdDim,
(DIM_Z == PERIODIC ? FOURIER : REAL));
static blitz::TinyVector<int,3>
local_lbounds(alloc_lbound(szx,szy,szz)),
local_extent(alloc_extent(szx,szy,szz));
static blitz::GeneralArrayStorage<3>
local_storage(alloc_storage(szx,szy,szz));
static DTArray temp_pv(local_lbounds,local_extent,local_storage),
temp_a(local_lbounds,local_extent,local_storage),
temp_b(local_lbounds,local_extent,local_storage);
/* Initialization */
if (Nx == 0 || Ny == 0 || Nz == 0) {
Nx = szx; Ny = szy; Nz = szz;
}
assert (Nx == szx && Ny == szy && Nz == szz);
/* x-vorticity is w_y - v_z */
/* b_x*(w_y - v_z) */
temp_pv = 0;
if (szx > 1) { // b_x
if (DIM_X == PERIODIC) {
deriv_fft(rho,trans_x,temp_b);
temp_b *= 2*M_PI/Lx;
} else if (DIM_X == FREE_SLIP) {
deriv_dct(rho,trans_x,temp_b);
temp_b *= M_PI/Lx;
} else {
assert(DIM_X == NO_SLIP);
deriv_cheb(rho,trans_x,temp_b);
temp_b *= -2/Lx;
}
}
if (szy > 1) { // w_y
if (DIM_X == PERIODIC) {
deriv_fft(w,trans_y,temp_a);
temp_pv += temp_a*temp_b*(2*M_PI/Ly);
} else if (DIM_X == FREE_SLIP) {
deriv_dct(w,trans_y,temp_a);
temp_pv += temp_a*temp_b*(M_PI/Ly);
} else {
assert(DIM_X == NO_SLIP);
deriv_cheb(w,trans_y,temp_a);
temp_pv += temp_a*temp_b*(-2/Ly);
}
}
if (szz > 1) { // v_z
if (DIM_Z == PERIODIC) {
deriv_fft(v,trans_z,temp_a);
temp_pv -= temp_a*temp_b*(2*M_PI/Lz);
} else if (DIM_Z == FREE_SLIP) {
deriv_dct(v,trans_z,temp_a);
temp_pv -= temp_a*temp_b*(M_PI/Lz);
} else {
assert(DIM_Z == NO_SLIP);
deriv_cheb(v,trans_z,temp_a);
temp_pv -= temp_a*temp_b*(-2/Lz);
}
}
// y-vorticity is u_z - w_x
/* b_y*(u_z - w_x) */
if (szy > 1) { // b_y
if (DIM_Y == PERIODIC) {
deriv_fft(rho,trans_y,temp_b);
temp_b *= 2*M_PI/Ly;
} else if (DIM_Y == FREE_SLIP) {
deriv_dct(rho,trans_y,temp_b);
temp_b *= M_PI/Ly;
} else {
assert(DIM_Y == NO_SLIP);
deriv_cheb(rho,trans_y,temp_b);
temp_b *= -2/Ly;
}
}
if (szz > 1) { // u_z
if (DIM_Z == PERIODIC) {
deriv_fft(u,trans_z,temp_a);
temp_pv += temp_a*temp_b*(2*M_PI/Lz);
} else if (DIM_Z == FREE_SLIP) {
deriv_dct(u,trans_z,temp_a);
temp_pv += temp_a*temp_b*(M_PI/Lz);
} else {
assert(DIM_Z == NO_SLIP);
deriv_cheb(u,trans_z,temp_a);
temp_pv += temp_a*temp_b*(-2/Lz);
}
}
if (szx > 1) { // w_x
if (DIM_X == PERIODIC) {
deriv_fft(w,trans_x,temp_a);
temp_pv -= temp_a*temp_b*(2*M_PI/Lx);
} else if (DIM_X == FREE_SLIP) {
deriv_dct(w,trans_x,temp_a);
temp_pv -= temp_a*temp_b*(M_PI/Lx);
} else {
assert(DIM_X == NO_SLIP);
deriv_cheb(w,trans_x,temp_a);
temp_pv -= temp_a*temp_b*(-2/Lx);
}
}
// And finally, vort_z is v_x - u_y
// b_z*(v_x - u_y)
if (szz > 1) { // b_z
if (DIM_Z == PERIODIC) {
deriv_fft(rho,trans_z,temp_b);
temp_b *= 2*M_PI/Lz;
} else if (DIM_Z == FREE_SLIP) {
deriv_dct(rho,trans_z,temp_b);
temp_b *= M_PI/Lz;
} else {
assert(DIM_Z == NO_SLIP);
deriv_cheb(rho,trans_z,temp_b);
temp_b *= -2/Lz;
}
temp_b += N0*N0;
temp_pv += f0*temp_b;
}
if (szx > 1) { // v_x
if (DIM_X == PERIODIC) {
deriv_fft(v,trans_x,temp_a);
temp_pv += temp_a*temp_b*(2*M_PI/Lx);
} else if (DIM_X == FREE_SLIP) {
deriv_dct(v,trans_x,temp_a);
temp_pv += temp_a*temp_b*(M_PI/Lx);
} else {
assert(DIM_X == NO_SLIP);
deriv_cheb(v,trans_x,temp_a);
temp_pv += temp_a*temp_b*(-2/Lx);
}
}
if (szy > 1) { // u_y
if (DIM_Y == PERIODIC) {
deriv_fft(u,trans_y,temp_a);
temp_pv -= temp_a*temp_b*(2*M_PI/Ly);
} else if (DIM_Y == FREE_SLIP) {
deriv_dct(u,trans_y,temp_a);
temp_pv -= temp_a*temp_b*(M_PI/Ly);
} else {
assert(DIM_Y == NO_SLIP);
deriv_cheb(u,trans_y,temp_a);
temp_pv -= temp_a*temp_b*(-2/Ly);
}
}
e_pv = &temp_pv;
return;
}
// Global arrays to store quadrature weights
Array<double,1> _quadw_x, _quadw_y, _quadw_z;
// Compute quadrature weights
void compute_quadweights(int szx, int szy, int szz,
double Lx, double Ly, double Lz,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z) {
_quadw_x.resize(split_range(szx));
_quadw_y.resize(szy); _quadw_z.resize(szz);
if (DIM_X == NO_SLIP) {
blitz::firstIndex ii;
_quadw_x = 1;
for (int k = 1; k <= (szx-2)/2; k++) {
// From Trefethen, Spectral Methods in MATLAB
// clenshaw-curtis quadrature weights
_quadw_x -= 2*cos(2*k*M_PI*ii/(szx-1))/(4*k*k-1);
}
if ((szx%2))
_quadw_x -= cos(M_PI*ii)/((szx-1)*(szx-1)-1);
_quadw_x = 2*_quadw_x/(szx-1);
if (_quadw_x.lbound(firstDim) == 0) {
_quadw_x(0) = 1.0/(szx-1)/(szx-1);
double Lx, double Ly, double Lz,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z) {
_quadw_x.resize(split_range(szx));
_quadw_y.resize(szy); _quadw_z.resize(szz);
if (DIM_X == NO_SLIP) {
blitz::firstIndex ii;
_quadw_x = 1;
for (int k = 1; k <= (szx-2)/2; k++) {
// From Trefethen, Spectral Methods in MATLAB
// clenshaw-curtis quadrature weights
_quadw_x -= 2*cos(2*k*M_PI*ii/(szx-1))/(4*k*k-1);
}
if ((szx%2))
_quadw_x -= cos(M_PI*ii)/((szx-1)*(szx-1)-1);
_quadw_x = 2*_quadw_x/(szx-1);
if (_quadw_x.lbound(firstDim) == 0) {
_quadw_x(0) = 1.0/(szx-1)/(szx-1);
}
if (_quadw_x.ubound(firstDim) == (szx-1)) {
_quadw_x(szx-1) = 1.0/(szx-1)/(szx-1);
......@@ -361,3 +618,4 @@ const blitz::Array<double,1> * get_quad_z() {
}
return &_quadw_z;
}
......@@ -32,6 +32,35 @@ void vorticity(TArrayn::DTArray & u, TArrayn::DTArray & v,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z);
/* Compute dudx */
void dudx(TArrayn::DTArray & u, TArrayn::DTArray * & out,
double Lx, double Ly, double Lz,
int szx, int szy, int szz,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z);
/* Compute dudy */
void dudy(TArrayn::DTArray & u, TArrayn::DTArray * & out,
double Lx, double Ly, double Lz,
int szx, int szy, int szz,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z);
/* Compute dudz */
void dudz(TArrayn::DTArray & u, TArrayn::DTArray * & out,
double Lx, double Ly, double Lz,
int szx, int szy, int szz,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z);
/* Compute ertel pv */
void ertel_pv(TArrayn::DTArray & u, TArrayn::DTArray & v,
TArrayn::DTArray & w, TArrayn::DTArray & rho,
TArrayn::DTArray * & e_pv, double f0, double N0,
double Lx, double Ly, double Lz,
int szx, int szy, int szz,
NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
NSIntegrator::DIMTYPE DIM_Z);
// Quadrature weights
void compute_quadweights(int szx, int szy, int szz,
......
......@@ -21,940 +21,1483 @@
namespace TArrayn {
using blitz::Array;
using blitz::TinyVector;
using blitz::Range;
using blitz::cast;
using std::string;
using std::ofstream;
using std::swap;
using std::cerr; using std::endl;
using namespace Transformer;
/* Compute the derivative of an array, using the Fourier
transform. This function assumes that the array is periodic over a
physical domain of length 2*pi. This is generally not true, but
correcting for this is a simple division (/L) in the calling code, and
best done there. */
void deriv_fft(DTArray & source, Trans1D & tform, Array<double, 3> & dest) {
/* Sanity checking -- the source and destination arrays should have
the same ranges, and the dimension of interest should start at 0 */
Dimension dim = tform.getdim();
assert(all(source.lbound() == dest.lbound()));
assert(all(source.ubound() == dest.ubound()));
/* Perform first FFT */
// fprintf(stderr,"&%g&\n",pvmax(source));
tform.forward_transform(&source, FOURIER);
Array<double,1> kvec = tform.wavenums();
// Get the complex temporary array from the transformer
CTArray * temp = tform.get_complex_temp();
//fprintf(stderr,"&%g&\n",psmax(max(abs(*temp))));
assert(temp->lbound(dim) == 0);
if (dest.extent(dim) % 2 == 0 && max(kvec) == tform.max_wavenum()) {
/* If we have an even number of points, the Nyquist frequency is
representable on the discrete grid. Unfortunately, the -derivative-
of the Nyquist frequency is purely imaginary, which we cannot
store in our real-valued destination. Thus, zero out this highest
frequency. */
kvec(temp->ubound(dim)) = 0;
}
/* After the Fourier transform, we can take the derivative by multiplying
each wavenumber by I*k. Since we're only working in one dimension
here, we can do this multiplication a plane at a time. */
/* The Blitz RectDomain class allows for simple specification of a slice */
//fprintf(stderr,"&%g&\n",psmax(max(abs(*temp))));
blitz::RectDomain<3> slice = temp->domain(); /* indexing */
double norm = tform.norm_factor(); // Normalization constant
for(int kk = temp->lbound(dim); kk <= temp->ubound(dim); kk++) {
slice.lbound(dim) = kk;
slice.ubound(dim) = kk;
(*temp)(slice) *= complex<double>(0,1*kvec(kk)/norm);
}
//fprintf(stderr,"&%g&\n",psmax(max(abs(*temp))));
/* And take the inverse Fourier transform to the destination array */
tform.back_transform(&dest, FOURIER);
//fprintf(stderr,"&%g&\n",pvmax(dest));
}
/* Derivatives for the cosine and sine transforms are nearly the same, so
to avoid wasted code they can both be rolled up into a single, "real-
transform" derivative function. */
namespace { /* Nested anonymous namespace for file-local linkage */
inline void deriv_real(DTArray & source, Trans1D & tform,
blitz::Array<double, 3> & dest, S_EXP t_type) {
assert(all(source.lbound() == dest.lbound()));
assert(all(source.ubound() == dest.ubound()));
/* Get transformed dimension */
Dimension dim = tform.getdim();
S_EXP b_type; // Backwards transform type
switch (t_type) {
case SINE:
b_type = COSINE; break;
case COSINE:
b_type = SINE; break;
default: // Shouldn't reach this
// assert(0 && "Invalid real derivative transform specified");
abort();
}
/* transform */
tform.forward_transform(&source, t_type);
blitz::firstIndex ii;
Array<double, 1> kvec = tform.wavenums();
DTArray & temp = *(tform.get_real_temp());
assert(temp.lbound(dim) == 0);
// std::cout << temp;
blitz::RectDomain<3> slicel=temp.domain(), slicer=temp.domain();
double norm = tform.norm_factor(); // Normalization constant
if (t_type == COSINE) {
/* The derivative of the cosine is a sine, and the constant frequency
(first cosine term) is completely absent in the sine transform.
Therefore, in the differentiation process we must shift our
wavenumbers DOWN one. */
for (int kk = 0; kk < temp.ubound(dim); kk++) {
slicel.ubound(dim) = kk; slicel.lbound(dim) = kk;
slicer.lbound(dim) = kk+1; slicer.ubound(dim) = kk+1;
temp(slicel) = -temp(slicer) * kvec(kk+1) / norm;
}
/* And zero out the highest frequency, which does not have a
corresponding cosine term */
temp(slicer) = 0;
} else {
/* Derivative of the sine is a cosine, and the problem is reversed.
The highest sine term has no equivalent cosine, so that frequency
is ignored. In the differentiation process, the frequencies are
shifted UP one (the constant term becomes zero). */
for (int kk = temp.ubound(dim); kk > 0; kk--) {
slicel.lbound(dim) = kk; slicel.ubound(dim) = kk;
slicer.lbound(dim) = kk-1; slicer.ubound(dim) = kk-1;
temp(slicel) = temp(slicer) * kvec(kk-1) / norm;
}
/* Zero out the constant frequency term */
temp(slicer) = 0;
}
/* And transform into the destination array */
// std::cout << temp;
tform.back_transform(&dest, b_type);
}
} /* End anonymous local namespace */
/* Now, the cosine and sine derivatives are just thin wrappers around the
above deriv_real. */
void deriv_dct(DTArray & source, Trans1D & tform, Array<double, 3> & dest) {
deriv_real(source,tform,dest,COSINE);
}
void deriv_dst(DTArray & source, Trans1D & tform, Array<double, 3> & dest) {
deriv_real(source,tform,dest,SINE);
}
void deriv_cheb(DTArray & source, Trans1D & tform, Array<double,3> & dest) {
/* Chebyshev differentiation */
/* Unlike other differentiations, chebyshev differentiation is defined
by a recurrance relation:
D = d/dx(F)
D_(N-1) = 0;
D_k = D_(k+2) + 2*(k+1)*F_(k+1)
modulo possible normalization at the zero frequency. This implies
that we're going to have to store a couple temporary array slices
around -- for F_(k) and F_(k+1) */
Dimension dim = tform.getdim();
assert(all(source.lbound() == dest.lbound()));
assert(all(source.ubound() == dest.ubound()));
if (source.extent(dim) == 1) {
dest = 0; return;
}
/* Transform to Chebyshev domain */
tform.forward_transform(&source, CHEBY);
/* Get wavenumbers */
// Array<double,1> kvec(tform.wavenums());
/* Note -- optimizing for 1D Chebyshevs, and getting the wavenumbers
here every derivative involves mallocing and bad things. We know
what the wavenumbers actually will be, so just substitute into
the expression */
/* And normalization factor */
double norm = tform.norm_factor();
/* Get the real temporary array */
DTArray & temparray = *(tform.get_real_temp());
// cout << temparray;
using blitz::Array;