Commit 58fbdb5d authored by Christopher Subich's avatar Christopher Subich
Browse files

Initial commit

This is the initial commit of SPINS to a git repository, coming from
an essentially unmanaged environment.  While this is a *complete*
archive, this is probably not the most useful form for development
going forward.  Notably, the future management should be:

1) Pare down the case files

2) Establish branches for individual users.  Not as a long-term goal,
   but instead to keep user-specific case files from lingering in
   the main repository long after their use is obsolete.

2b) Establish semipermanent branches for typical users.  A release
   branch should include the basic documentation cases, a benchmark
   case should include benchmarking cases, and there are probably
   others that aren't coming to mind now.

3) Add helpful MATLAB processing scripts to the repository.  Too
   much of that is ad-hoc now, completely unmanaged.

4) For papers/etc, establish tagged versions.
parents
# This file contains a list of things to *not* add to the git repository
# Object files
*.o
# Executable files
*_x
*.x
# Default output files, in case a case gets run in the main directory
# Grids
xgrid*
ygrid*
zgrid*
# Matlab reader files
*reader.m
# Numbered output files
*.[0-9]*
# Any vim swap files
*.swp
#include "BaseCase.hpp"
#include "NSIntegrator.hpp"
#include "TArray.hpp"
//using namespace TArray;
using namespace NSIntegrator;
using blitz::Array;
using std::vector;
/* 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();
}
int BaseCase::size_x() const {
return size_cube();
}
int BaseCase::size_y() const {
return size_cube();
}
int BaseCase::size_z() const {
return size_cube();
}
double BaseCase::length_x() const {
return length_cube();
}
double BaseCase::length_y() const {
return length_cube();
}
double BaseCase::length_z() const {
return length_cube();
}
DIMTYPE BaseCase::type_x() const {
return type_default();
}
DIMTYPE BaseCase::type_y() const {
return type_default();
}
DIMTYPE BaseCase::type_z() const {
return type_default();
}
DIMTYPE BaseCase::type_default() const {
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;
}
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;
}
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;
}
bool BaseCase::tracer_bc_forcing() const {
return false;
}
bool BaseCase::is_mapped() const { // Whether this problem has mapped coordinates
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;
}
/* Physical parameters */
double BaseCase::get_visco() const {
return 0;
}
double BaseCase::get_diffusivity(int t) const {
return 0;
}
/* Initialization */
double BaseCase::init_time() const {
return 0;
}
void BaseCase::init_tracers(vector<DTArray *> & tracers) {
/* Initalize tracers one-by-one */
if (numtracers() == 0) return; // No tracers, do nothing
assert(numtracers() == 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
}
/* 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;
}
}
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);
}
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;
}
/* Analysis */
void BaseCase::analysis(double t, DTArray & u, DTArray & v, DTArray & w,
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]));
}
}
/* BaseCase -- basic skeletal framework for a functional NSIntegrator usercode.
Eliminates some of the boilerplate. */
#ifndef BASECASE_HPP
#define BASECASE_HPP 1
#include <blitz/array.h>
#include "TArray.hpp"
#include "NSIntegrator.hpp"
using namespace TArrayn;
using namespace NSIntegrator;
using blitz::Array;
using std::vector;
class BaseCase {
/* To reduce boilerplate, wrap some of the long functions, only calling
them if actually used by usercode. For example, a tracer-free code
does not need the long-form forcing function -- nor does one with
only passive tracers */
public:
/* Tracers */
virtual int numActive() const; // Number of active tracers
virtual int numPassive() const; // Number of passive tracers
virtual int numtracers() const; // Number of tracers (total)
/* Grid generataion */
virtual int size_x() const; // Grid points in x
virtual int size_y() const; // Grid points in y
virtual int size_z() const; // Grid points in z
virtual int size_cube() const { abort(); }; // Special case -- N*N*N
virtual double length_x() const; // Length in x
virtual double length_y() const; // Length in y
virtual double length_z() const; // Length in z
virtual double length_cube() const {abort();}; // Special case -- L*L*L
virtual DIMTYPE type_x() const; // Expansion type in x
virtual DIMTYPE type_y() const; // Expansion type in y
virtual DIMTYPE type_z() const; // Expansion type in z
virtual DIMTYPE type_default() const; // Default expansion type
/* Functions to provide for custom tracer boundaries, along any
dimension whose expansion type supports them (generally that
would be chebyshev-type). The default tracer-BC behaviour
will be to give neumann-type BCs. Dirichlet-type BCs are possible
with a sine-based expansion rather than cosine-based expansion
and have been hacked in with a global variable, but this mechanism
is far more general. One will expect an error/abort if an improper
BC-type is requested for a boundary that doesn't actually support it
(namely a Robin-type BC on a non-Chebyshev dimension */
// Further note -- we want the SAME boundary condition at the minimum and
// maximum of the domain.
virtual void tracer_bc_x(int tracernum, double & dirichlet, double & neumann) const;
virtual void tracer_bc_y(int tracernum, double & dirichlet, double & neumann) const;
virtual void tracer_bc_z(int tracernum, double & dirichlet, double & neumann) const;
// Whether ANY tracers will have nonzero boundary conditions. If this is true
// then the user forcing code is responsible for calculating/applying the proper
// RHS at the boundaries. If this is false (default), then the prior behaviour
// of the integrator code zeroing BCs after forcing is retained.
virtual bool tracer_bc_forcing() const;
virtual bool is_mapped() const; // Whether this problem has mapped coordinates
// 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
virtual void do_mapping(DTArray & xgrid, DTArray & ygrid, DTArray & zgrid);
/* Physical parameters */
virtual double get_visco() const; // Physical viscosity
virtual double get_diffusivity(int tracernum) const; // Diffusivity
/* Initialization */
virtual double init_time() const; // Initialization time
virtual void init_tracers(vector<DTArray *> & tracers);
virtual void init_vels(DTArray & u, DTArray & v, DTArray & w) { abort();};
virtual void init_tracer(int t_num, DTArray & tracer) { abort();}; // single-tracer
/* Numerical checks */
virtual double check_timestep(double step, double now);
// Get incoming gradient operator, for differentials in analysis. This is a null-
// op unless the user code cares, in which case it will override this to store
// the gradient.
virtual void set_grad(Grad * in_grad) {return;};
/* Forcing */
/* Big, ultra-general forcing function */
virtual void 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);
/* No-active-tracers specialization */
virtual void passive_forcing(double t, DTArray & u, DTArray & u_f,
DTArray & v, DTArray & v_f,
DTArray & w, DTArray & w_f);
/* Independent-of-current-velocity no-tracer forcing */
virtual void stationary_forcing(double t, DTArray& u_f, DTArray& v_f,
DTArray& w_f);
/* If there are active tracers, split V and T focing */
virtual void vel_forcing(double t, DTArray& u_f, DTArray& v_f,
DTArray& w_f, vector<DTArray *> & tracers) {abort();};
virtual void tracer_forcing(double t, DTArray & u,
DTArray & v, DTArray & w,
vector<DTArray *> & tracers_f) {abort();};
/* Analysis and writing */
virtual void analysis(double t, DTArray & u, DTArray & v, DTArray & w,
vector<DTArray *> tracer, DTArray & pres); // General Analysis
virtual void analysis(double t, DTArray & u, DTArray & v, DTArray & w,
vector<DTArray *> tracer); // Less pressure
virtual void vel_analysis(double t, DTArray & u, DTArray & v,
DTArray & w) {abort();}; // Velocity analysis
virtual void tracer_analysis(double t, int t_num, DTArray & tracer) {abort();};
// Single-tracer analysis
};
extern template class FluidEvolve<BaseCase>;
typedef FluidEvolve<BaseCase> EasyFlow; // Explicit template instantiation
#endif
This diff is collapsed.
/* ESolver.hpp -- header file for ESolver class and any associated functions.*/
#ifndef ESOLVER_HPP
#define ESOLVER_HPP 1
/* ESolver is essentially the "guts" of SPINS, designed to transparently solve
problems of the form:
f(x,y,z)*Lap(u) - G*u = h(x,y,z)
These problems arise (twice!) at each timestep in the discretization of the
Navier-Stokes equations. Firstly, under a splitting method pressure is
calculated as
Lap(pressure) = div(u,v,w),
where (u,v,w) are pressure-free projected velocities. Without interaction
with boundaries, pressure acts (through its gradient) to make the velocities
divergence-free. In this case, clearly the constant G (above) is 0, making
this a Poisson problem.
Secondly, with implicit diffusion, the momentum problem becomes:
mu(x,y,z)*Lap(u) - u/dt = f(x,y,z),
a Helmholtz problem where mu(x,y,z) is the viscosity. If there is
additionally a sponge layer, the diffusion may be anisotropic. This
is not conceptually different, but requires an equation of the form:
muf(x,y,z)*uxx + mug(x,y,z)*uyy + muh(x,y,z)*uzz - u/dt = f(x,y,z).
In the general case, ESolver needs to deal with a volume split over
multiple processors and solutions that are not separable -- it must deal
with the 3D domain as a single entity. The only way to do this is to build
a full finite-difference preconditioner (to be written, using multigrid)
for the operator and solve with an iterative scheme (GMRES).
Separable dimensions are those that reduce the differential operator to
an algebraic one under a Fourier Transform (DFT/DCT/DST). The resulting
problem is then a repeated one in 1 or 2 dimensions, which is significantly
easier to solve. Indeed, if the arallelized dimension is separable, then
the problem of parallel multigrid can be avoided entirely.
This suggests that ESolver breaks up into a couple different cases:
1) Triply-separable with constant Jacobian, so take the spectral transform
and divide by k^2+l^2+m^2
2) Doubly-separable, with (possibly mapped) Chebyshev vertical expansion.
In this case, we have only one dimension at a time to consider, so take
the spectral transform on x/y and solve the tridiagonal FD problem with
k^2+l^2 on the diagonal
3) Singly-separable (spectral along y), requires a 2D multigrid problem.
4) None-separable, requires a 3D multigrid problem. For now we won't touch
this case, since writing 3D multigrid will be a bit finicky (although not
mathematically harder)
*/
#include "TArray.hpp"
#include "Splits.hpp"
#include "Parformer.hpp"
#include <blitz/array.h>
#include "grad.hpp"
namespace ESolver {
using TArrayn::DTArray;
using TArrayn::CTArray;
using namespace Transformer;
class ElipSolver {
/* ESolver 0.1 -- triply-periodic FFT, solving:
Lap(u) - M*u = f(x,y,z) */
/* ElipSolver 0.2 -- adding sin/cos expansions */
public:
/* The mechanics of real/spectral transformation is handled by the
TransWrapper class. This means that ElipSolver doesn't have to
directly handle temporary arrays, nor does it have to calculate
wavenumbers (and normalization factors) by itself. */
/* Note that the gradient operator is now specified. We can't really
leverage the derivatives in it since we don't -take- derivatives
over the entire field, but we can use the jacobian info inside
for our calculations. (a 3D multigrid -could- use grad directly,
but we're not writing that now.) */
ElipSolver(double M,TransWrapper * spec_space, TArrayn::Grad * in_grad);
~ElipSolver(); // Destructor
/* Solve the Helmholtz/Poisson problem with the given rhs, storage
in output. */
void solve(DTArray & rhs, DTArray & output,
S_EXP type_x=FOURIER, S_EXP type_y=FOURIER,
S_EXP type_z=FOURIER, double zbc_a=0, double zbc_b = 0,
double xbc_a=0, double xbc_b=0,
double ybc_a=0, double ybc_b=0);
/* For timestepping: reset the Helmholtz parameter */
void change_m(double M);
private:
double M; // Helmholtz parameter
TArrayn::Grad * gradient;
S_EXP t_types[3]; // Transform types
/* Specialization of solve to the case when all three
dimensions are non-mapped trig expansions. The
derivatives become algebraic after transform, meaning
no GMRES required */
void threespec_solve(DTArray & rhs, DTArray & output,
S_EXP type_x, S_EXP type_y, S_EXP type_z);
/* Specialization for one Chebyshev dimension (vertical),
including (identical) BC types for the top and bottom */
void chebz_solve(DTArray & rhs, DTArray & output,
S_EXP type_x, S_EXP type_y, double zbc_a, double zbc_b);
/* Specialization for the 2D GMRES/Multigrid solve,
with identical BCs at all solid interfaces */
void twodeemg_solve(DTArray & rhs, DTArray & output,
S_EXP type_x, S_EXP type_y, double zbc_a, double zbc_b, double xbc_a, double xbc_b);
TransWrapper * spec_transform;
// Data members for the 1D (z) Chebyshev case, which requires
// a parallel transpose for data ordering
Transposer<complex<double> > * cTransposer;
Transposer<double> * rTransposer;
DTArray * rtransyz;
CTArray * ctransyz;
// Data members for the 2D multigrid solve, which operates on a
// Nx1xM 3D array
DTArray * twod_u, * twod_f;
};
} // End namespace
#endif // ESOLVER_HPP
CC = mpicc
CXX = mpicc
LD = mpicc
# -ip -ipo seem to give about 10% better speed, but they also double compile
# times. For now, disable them.
BLITZ_DEBUG = -g -DBZ_DEBUG #-ftrapuv -check-uninit
SLOW_OPTIM = -O3 -ffast-math #-ip -ipo
# Compile with all warnings, but disable some remarks. Notably:
# 1572-- "unreliable" == comparison of floating points
# 981 -- unspecified order of evaluation of operands
# 444 -- destructor is not virtual
# 383 -- value copied to temporary, reference to temporary used
# 869 -- parameter not referenced
INCLUDE_DIRS = -I../UMFPACK/Include -I../AMD/Include -I../UFconfig
CFLAGS = $(INCLUDE_DIRS) $(BLITZ_DEBUG) $(SLOW_OPTIM) #-Wall -wd981 -wd444 -wd1572 -wd383 -wd869
# disabling notifications for
# 11000 - performing multi-file optimizations
# 11005 - generating object file <temporary name>.o
LIB_DIRS = -L../UMFPACK/Lib -L../AMD/Lib #-L/opt/scsl/lib
LAPACK_LIB = -llapack
UMFPACK_LIB = -lumfpack -lamd
LDFLAGS = $(SLOW_OPTIM) #-wd11000,11005
LDLIBS = $(LIB_DIRS) -lfftw3 -lblitz -lm -lstdc++ -lmpi $(LAPACK_LIB) $(UMFPACK_LIB)
.PHONY: all
all: tests/test_deriv_x tests/test_write_x tests/test_esolve_x tests/test_heat_x tests/test_ode_x tests/test_ns_x
.PHONY: clean
clean:
rm -f *.o tests/*.o cases/*.o
objfiles: $(shell ls *.cpp | sed -e 's/cpp/o/g')
NSIntegrator.o: NSIntegrator.cpp NSIntegrator_impl.cc
tests/test%.o: tests/tests%.cpp
$(CXX) $(CFLAGS) -o $@ -c $<
tests/test%_x: tests/test%.o TArray.o Parformer.o T_util.o ESolver.o Timestep.o NSIntegrator.o BaseCase.o Science.o Splits.o Par_util.o Split_reader.o gmres.o gmres_1d_solver.o gmres_2d_solver.o grad.o multigrid.o
$(LD) $(LDFLAGS) -o $@ $^ $(LDLIBS)
cases/%.o: cases/%.cpp NSIntegrator_impl.cc NSIntegrator.hpp
$(CXX) $(CFLAGS) -o $@ -c $<
nonhydro_x: nonhydro_sw.o TArray.o T_util.o Parformer.o Splits.o Par_util.o Split_reader.o
$(LD) $(LDFLAGS) -o $@ $^ $(LDLIBS)
derek_x: derek.o TArray.o T_util.o Parformer.o ESolver.o Timestep.o Splits.o Par_util.o Split_reader.o gmres.o gmres_1d_solver.o gmres_2d_solver.o grad.o multigrid.o
$(LD) $(LDFLAGS) -o $@ $^ $(LDLIBS)
cases/%.x: cases/%.o TArray.o T_util.o Parformer.o ESolver.o Timestep.o NSIntegrator.o BaseCase.o Science.o Splits.o Par_util.o Split_reader.o gmres.o gmres_1d_solver.o gmres_2d_solver.o grad.o multigrid.o
$(LD) $(LDFLAGS) -o $@ $^ $(LDLIBS)
cases/%_x: cases/%.o TArray.o T_util.o Parformer.o ESolver.o Timestep.o NSIntegrator.o BaseCase.o Science.o Splits.o Par_util.o Split_reader.o gmres.o gmres_1d_solver.o gmres_2d_solver.o grad.o multigrid.o
$(LD) $(LDFLAGS) -o $@ $^ $(LDLIBS)
%.o : %.cpp *.hpp
$(CXX) $(CFLAGS) -o $@ -c $<
#include "NSIntegrator.hpp"
#include "BaseCase.hpp"
/* And default global values for filtering */
double f_cutoff = 0.6, f_order = 2, f_strength = 20;
// Hack -- if a density tracer is really a perturbation density, it should
// have zero boundary conditions on a cosine/free-slip grid. Otherwise,
// the zero-derivative BC is more appropriate. Short of changing the function
// interface to specify this per-tracer, use a single global boolean
bool zero_tracer_boundary = false;
namespace NSIntegrator{
S_EXP get_trans_type(DIMTYPE dtype, Sym stype) {
switch (dtype) {
case PERIODIC:
return FOURIER;
case FREE_SLIP:
if (stype == EVEN) return COSINE;
if (stype == ODD) return SINE;
abort(); break;
case NO_SLIP:
return CHEBY;
default:
abort(); return NONE;
}
}
}
template class FluidEvolve<BaseCase>; // Explicit template instantiation
/* FluidEvolve.hpp -- control class for Navier-Stokes (incompressible)
integration */
#ifndef NSINTEGRATOR_HPP
#define NSINTEGRATOR_HPP 1
/* FluidEvolve is the "physical guts" of SPINS, in the way that ESolver is
the "numerical guts". FluidEvolve controls:
1) The maximum allowable timestep, via CFL conditions.
2) Projection of velocities
3) Computation of pressure (via ESolver)
4) Constructing Helmholtz solutions for viscosity (via ESolver)
5) Stepping any tracers (to be implemented)
It's all that and a bag of chips. It's even a template class, to allow
for plug-in user control code. This code is responsible for:
1) Setting the maximum timestep, for cosmetic or plotting purposes
(That is, user code may choose a timestep lower than what FluidEvolve
gives it, if it would result in desirable timelevels being hit
exactly.)
2) Velocity forcing, /including/ gravitational and rotational forces.
(Why? Because especially gravitational forces depend on exactly how
you specify density. With a Boussinesq approximation, it's common
to have a "background" stratification, which shows up in a forcing
term in the density equation.)
3) All user-readable IO, including writing fields and any scalar
diagnostics. Indeed, what is "interesting" here varies so much
from problem to problem that this is the original motivation for
spinning control code into a user-supplied class. */
#include <vector>
#include <set>
#include "ESolver.hpp"
#include "Timestep.hpp"
#include "TArray.hpp"
#include "T_util.hpp"
#include "Parformer.hpp"
#include <stdio.h>
#include "Par_util.hpp"
#include "grad.hpp"
// Global values to control filtering
extern double f_cutoff, f_order, f_strength;
extern bool zero_tracer_boundary;
namespace NSIntegrator {
using namespace TArrayn;
using namespace Timestep;
using namespace ESolver;
using std::vector;
using std::set;
using namespace Transformer;
// List supported physical dimension types (Fourier, free slip, no slip, mapped)
enum DIMTYPE {
PERIODIC,
FREE_SLIP, // Cosine expansion
NO_SLIP // Chebyshev expansion
};
enum Sym { // Symmetry
EVEN, ODD };
template <class Control> class FluidEvolve {
public:
/* This will have to include local/global ranges for parallelization*/
int szx, szy, szz; // Sizes
DIMTYPE tx, ty, tz; // Dimension types
/* This will have to include ranges (arrays) for topography */
double Lx, Ly, Lz; // domain lengths
double visco;
Control * usercode;
/* Parameters for allocating arrays, given parallel splitting */
blitz::TinyVector<int,3> local_lbounds, local_extents;
blitz::GeneralArrayStorage<3> local_storage;
/* Arrays for velocities and their previous timelevel forcings */
Stepped<DTArray> us, us_forcing, vs, vs_forcing, ws, ws_forcing;
/* Arrays for tracers. Since we have a variable number of them,
we'll bind the Steppeds up in a vector. */