#ifndef MULTIGRID_HPP
#define MULTIGRID_HPP 1
#include
#include
#include "umfpack.h"
/* Header file for our required 2D multigrid solver, operating on a (split) grid
with the numerical Jacobian entries (from the spectral model) providing varying
coefficients */
extern double coarse_solve_time, redblack_time, apply_op_time;
extern int coarse_solve_count, redblack_count, apply_op_count;
enum SYM_TYPE {
SYM_NONE=0, // No symmetry
SYM_PERIODIC, // Periodic
SYM_EVEN,
SYM_ODD
} ;
enum CYCLE_TYPE { // Multigrid cycle type
CYCLE_V, // down-up, no recursion
CYCLE_W, // down-up-down-up with full recursion
CYCLE_F, // down-up-down-up with recursion only on the first half
CYCLE_NONE // Debugging -- relax only
};
class MG_Solver {
public:
/* Constructor. Gets x and z values of the numerical (perfectly orthogonal)
grid, for construction of the 1D finite difference operators. These arrays
must also be for the full domain, since boundary points are special.
If symmetry_problem is true, then xvals is extended by two to allow for
a ghost-point type handling of the left and right boundaries */
MG_Solver(blitz::Array xvals, blitz::Array zvals,
SYM_TYPE symmetry_problem, MPI_Comm comm = MPI_COMM_WORLD);
/* Set the specific form of symmetry. This is not in the constructor because
it might change between invocations on the same grid. As an example, on
a rectanguar grid with free-slip conditions, u (orthogonal to the boundary)
is odd symmetry while w is even. Pressure likewise has even symmetry */
void set_x_symmetry(SYM_TYPE type);
/* Note that after this point all arrays specified are split by processor */
/* Problem setup, which takes the several coefficients. With the exception
of DC, these will probably not change between invocations. For
the multilevel algorithm, this can call downwards recursively with
coarsened versions of these operators */
void problem_setup(blitz::Array & Uxx, blitz::Array & Uzz,
blitz::Array & Uxz, blitz::Array & Ux,
blitz::Array & Uz);
/* Setting the DC term. This will change based on timestep, diffusivity,
or spanwise wavenumber (in the case of a reduced 3D problem). The
general case has this varying, but I much prefer not to deal with that
right now */
void helmholtz_setup(double c);
/* BC setup; so many terms are necessary because physical grid lines might not
intersect the boundary at a right angle. Indeed, they won't at the bottom
even in the case of a hill. To avoid a duplicate function spec, this
fuction gets called twice with dim=0 for x and dim=1 for z*/
void bc_setup(int dim,blitz::Array u_min,
blitz::Array uz_min,blitz::Array ux_min,
blitz::Array u_max, blitz::Array uz_max,
blitz::Array ux_max);
/* Perform a multigrid cycle, or at the coarse level do a direct solve
via UMFPACK. Takes as input the desired cycle type (F-cycle or V-cycle),
and input residual f, and returns as output the (approximate) u.
For the coarse level, this does the coarse solve regardless of cycle type.
The bonus_in and bonus_out terms are required for indefinite problems, where
u=constant is a solution to A*u=0. In that case, we regularize the problem
by requiring mean(u)==bonus_in, and in return the f applied to get that
value is really f-bonus_out (on interior terms). This seems a bit picky,
but it's required to make gmres behave */
/* This is the public facing implementation, which load-balances f and u
and calls cycle_p */
void cycle(CYCLE_TYPE cycle, blitz::Array & f, blitz::Array & u,
double bonus_in, double & bonus_out,int pre = 2, int mid = 2, int post = 2);
// protected:
/* Local array size */
int local_size_x;
int local_x_lbound, local_x_ubound;
int coarse_x_lbound, coarse_x_ubound;
int size_x, size_z;
// Parallelization parameters
int myrank, nproc;
// Parameters for the coarse-grid solving
bool coarse_symbolic_ok, // Whether the coarse grid symbolic factorizing is okay
coarse_numeric_ok, // Same with the numerical factor
// And parameters to control when a symbolic factor is not okay
any_dxz, // Whether there's any nonzero dxz terms in the domain
bc_tangent, // Whether there's any tangent derivatives in the BCs
bc_normal; // Same, normal derivatives
// numeric and symbolic factors from UMFPACK
void * numeric_factor, * symbolic_factor;
// Arrays for the sparse representation of the coarse operator, see UMFPACK
// documentation for the full information. The short version is that
// the entry A_double[A_cols[i]+j] is in column i, and row A_rows[A_cols[i]+j]
int sparse_size; // number of entries in sparse array
int * A_cols, // Array of column indices
* A_rows; // Array of row index entries
double * A_double; // Array of entries in the sparse matrix
// Builds the sparse operator, storing proper entries in the above arrays
// and allocating if necessary
void build_sparse_operator();
// Performs consistency checks on whether a symbolic coarse grid factor
// is OK based on current conditions and the coarse-grid boolean parameters
// (see above)
void check_bc_consistency();
/* One dimensional FD operators */
/* Interpretation:
Op(i,j) is the operator for point i, referring to:
j=0 -- left (minus) neighbour
j=1 -- self
j=2 -- right (plus) neighbour
*/
blitz::Array Dx, Dxx, Dz, Dzz;
/* Local coefficients for the 2D operator */
blitz::Array uxx, uzz, uxz, ux, uz;
/* Local arrays for load-balanced f and u arrays, since those
may differ from input arrays */
blitz::Array f_balance, u_balance;
/* c, for the c*U term in the equation */
double helm_parameter;
/* Local bc coefficients */
blitz::Array u_left, u_right, ux_left, ux_right, uz_left, uz_right,
u_top, u_bot, ux_top, ux_bot, uz_top, uz_bot;
/* True if the problem is indefinite, such that A*(constant) = 0.
See comments for cycle() for the implications */
bool indefinite_problem;
/* Symmetry type */
SYM_TYPE symmetry_type;
/* Perform one red-black relaxation pass */
void do_redblack(blitz::Array & f, blitz::Array & u);
/* Apply the operator "forwards" */
void apply_operator(blitz::Array & u, blitz::Array & f);
/* Coarse u, f arrays for multilevel problem */
blitz::Array coarse_u, coarse_f;
/* Local, nonbalanced coarsened arrays. Since the array gets roughly halved,
we can either rebalance and then coarsen or coarsen and then rebalance.
This takes the latter approach */
blitz::Array local_coarse_2d;
blitz::Array local_coarse_1d;
/* If we coarsen an even number of points (odd number of
intervals), then the "include every other point" strategy
doesn't work (given that we keep the boundaries). We have
to thus pick an interval to keep between coarsening levels;
keeping the largest interval on the fine grid seems to be
an excellent choice */
int kept_interval;
bool coarsest_level;
/* (semi)-coarsen a variable on a 2D grid. The extra
boolean (even) tells the coarsening operator to
treat the problem as having even symmetry (if any)
rather than whatever has been specified. This is
useful for the coefficients and boundary conditions,
since they'll never have odd symmetry */
/* This operator does not rebalance -- it will simply write
into local_coarse_2d (1d for line version) */
void coarsen_grid(blitz::Array & q, bool even=false);
/* And do the same to a line. */
void coarsen_line(blitz::Array & q, bool even=false);
/* Interpolate (from local_coarse_2d) to full-sized array */
void interpolate_grid(blitz::Array & q_fine);
// Interpolate line commented out -- I don't see a need for it
// void interpolate_line(blitz::Array & q_coarse,
// blitz::Array & q_fine, bool even=false);
/* Local cycle, which implements the public function but assumes all relevant
arrays are properly load-balanced */
void _cycle(CYCLE_TYPE cycle, blitz::Array & f, blitz::Array & u,
double bonus_in, double & bonus_out,int pre, int mid, int post);
/* Rebalances an array; an array that is orignally spit in orig is
reassigned to balance, such that the global (compressed) array
is exactly the same */
// Arrays to hold origin, destination bounds for array rebalancing (coarse)
enum BalanceGroup {
FFine = 0, UFine = 1, FCoarse = 2, UCoarse = 3
};
int *to_lbounds[4], *to_ubounds[4], *tb_lbounds[4], *tb_ubounds[4];
bool rebalance_setup[4];
void rebalance_array(blitz::Array & orig,
blitz::Array & balance, MPI_Comm c, BalanceGroup btype);
void rebalance_line(blitz::Array & o,
blitz::Array & b, MPI_Comm c, BalanceGroup btype);
MPI_Comm my_comm;
MPI_Comm coarse_comm; // Communicator for coarse grid
MG_Solver * coarse_solver; // Coarse solver
};
/* Build a 1D finite-difference operator for arbitrary x-grid, see mkfd.m
as reference in MATLAB */
void get_fd_operator(blitz::Array & x, blitz::Array & Dx,
blitz::Array & Dxx);
/* Get the local splitting of N points in the multigrid context. This might (will)
differ from the splitting in spectral context, because we want a minimum of 2 lines
per processor, and preferably each processor will have an even number of lines
to minimize communication. e.g.
4 4 4 4 6 is preferable to
4 4 4 5 5, since the latter will mix up red/black identities
(RBRB-RBRB-RBRB-RBRB-RBRBRB vs.
RBRB-RBRB-RBRB-RBRBR-BRBRB -- the latter will involve a double-stall on
processor 4 for boundary data */
void get_local_split(int extent, int rank, int nproc, int & l_lbound, int & l_ubound);
/* Perform a line solve, which is the 1d problem:
A(z)*u_zz + B(z)*u_z + C(z) = f(z), with
a0*u+b0*u_z = f0 @ z=0 (bottom)
a1*u+b1*u_z = f1 @ z=1 (top)
Unlike the 1D gmres controller, we're not going to go to
great lengths to preserve a factored operator. Firstly,
there's no chance of an indefinite 1D problem here, and
secondly most of the time we're going to be doing a -lot-
of line solves, and the only way to cache things would
be to store all of those operators. Bad for memory. */
void line_solve(blitz::Array & u, blitz::Array & f,
blitz::Array & A, blitz::Array & B,
blitz::Array & C, double a0, double a1,
double b0, double b1, blitz::Array & Dz,
blitz::Array & Dzz);
#endif