File timestep.c

RCS Header: /cvsroot/rheoplast/timestep.c,v 1.98 2006/03/08 03:21:15 hazelsct Exp

This file includes an explicit timestepper and a semi-implicit one. The explicit timestepper can optimize for slow or fast communication; fast is really no different than traditional, but slow does asynchronous I/O to communicate the content of the interface nodes while computing the inner ones, and should work quite a bit better for bandwidth-limited Beowulf clusters. In practice, "fast" is quite a bit faster, because the short loops involved in computing the interface nodes in the asynchronous version slow things down considerably.
It also does ``constrained implicit'' timestepping, which differs from PETSc's timestepping algorithms in that where the latter requires that all functions be of the form:

du/dt = f(u), it is often necessary to have some fields described by ``constraint'' functions of the form: f(u = 0. For example, incompressible Navier-Stokes has the motion equations which can be written as the former time-derivatives, and the continuity equation which must be written in the latter constraint form. The implicit_steptime function allows for those two types of equations to be mixed.
Both of these steppers implement temporary fields, which are calculated from the solved field variables in order to simplify the programming, save time, and in the case of rheoplast, pass needed information between modules.
When I have some time, I'll consider implementing these within PETSc's TS timestepping object class, and then submitting a patch upstream.

Included Files

#include "timestep.h"
- #include </usr/lib/petsc/include/petscsnes.h>
- #include </usr/lib/petsc/include/petscda.h>
- #include </usr/lib/petsc/include/petscblaslapack.h>
- #include <stdlib.h>

Preprocessor definitions

#define __FUNCT__ "xmin_symm"

#define __FUNCT__ "xmax_symm"

#define __FUNCT__ "ymin_symm"

#define __FUNCT__ "ymax_symm"

#define __FUNCT__ "zmin_symm"

#define __FUNCT__ "zmax_symm"

#define __FUNCT__ "explicit_steptime"

#define __FUNCT__ "FuncEvaluate"

#define __FUNCT__ "ts_implicit_function"

#define __FUNCT__ "ts_implicit_jacobian"

#define __FUNCT__ "implicit_steptime"

Local Variables

implicit_da
static DA implicit_da

implicit_jack
static Mat implicit_jack

un
static Vec un

temp_local
static Vec temp_local

Fun
static Vec Fun

Fun_local
static Vec Fun_local

Funp1
static Vec Funp1

Funp1_local
static Vec Funp1_local

unm1
static Vec unm1

unm2
static Vec unm2

unm3
static Vec unm3

unm4
static Vec unm4

unm5
static Vec unm5

current_time
static PetscScalar* current_time

temp_fields
static PetscScalar* temp_fields

implicit_deltat
static PetscScalar implicit_deltat

timestyle_coefficient
static PetscScalar timestyle_coefficient

implicit_timestyle
static EqStyle* implicit_timestyle

timestep_style
static TSStyle timestep_style

mixed_constraints
static PetscTruth** mixed_constraints

dof
static int dof

gxs
static int gxs

gys
static int gys

gzs
static int gzs

gxm
static int gxm

gym
static int gym

gzm
static int gzm

implicit_tps
static int implicit_tps

implicit_bcflags
static int implicit_bcflags

implicit_symmtypes
static int* implicit_symmtypes

implicit_dadims
static int implicit_dadims

Global Function explicit_steptime()

This is the explicit function, it does the timestepping with the slow- or fast-communication optimization. Note that the slow-communication optimization assumes that calculations can be done between the start and end of global to local scattering, that is during communication; performance testing seems to be showing that this is not the case, or else the gain from doing this is minimal... Also note: neither temporary field variables nor symmetry boundary conditions are implemented in the slow communication section!

int explicit_steptime ( DA theda, Vec globalreal, Vec localreal, int* currentstep, int timesteps, PetscScalar* currenttime, PetscScalar deltat, int tps, CommStyle comm, int bcflags, int* symmtypes, void* user, FILE* logfile )

int explicit_steptime: Returns 0 or error code: -1 for malloc error, -2 for non-periodic DA (only can use fully-periodic DAs for now).
DA theda: The DA object we're working on.
Vec globalreal: Global vector of all fields.
Vec localreal: Local vector of all fields.
int* currentstep: Entering: initial timestep number, return: final timestep number.
int timesteps: Number of timesteps to run.
PetscScalar* currenttime: Entering: initial time, return: final time.
PetscScalar deltat: Timestep size.
int tps: Number of temporary parameters to store in array.
CommStyle comm: Communication optimization to use (see typedef enum above).
int bcflags: Flags controlling boundary conditions.
int* symmtypes: Array of symmetry types for each dof (if NULL then assumes all are MIRROR_PLANE).
void* user: User parameters' structure.
FILE* logfile: Log file to fprintf debugging information to (or NULL).

Global Function implicit_steptime()

This function implements constrained implicit timestepping as described above.

int implicit_steptime ( DA theda, Vec unp1, Vec unp1_local, Mat jack, int* currentstep, int timesteps, PetscScalar* currenttime, PetscScalar deltat, int tps, EqStyle* timestyle, int bcflags, int* symmtypes, void* user, FILE* logfile )

int implicit_steptime: Returns 0 or error code: -1 for malloc error.
DA theda: The DA object we're working on.
Vec unp1: Global vector of all fields.
Vec unp1_local: Local vector of all fields.
Mat jack: Matrix to use as Jacobian.
int* currentstep: Entering: initial timestep number, return: final timestep number.
int timesteps: Number of timesteps to run.
PetscScalar* currenttime: Entering: initial time, return: final time.
PetscScalar deltat: Timestep size.
int tps: Number of temporary field variables to store in array.
EqStyle* timestyle: Constraint, time derivative, or blend nature of each field equation.
int bcflags: Flags controlling boundary conditions.
int* symmtypes: Array of symmetry types for each dof (if NULL then assumes all are MIRROR_PLANE).
void* user: User parameters' structure.
FILE* logfile: Log file to fprintf debugging information to (or NULL).

First we copy arguments into static variables, get DA info, allocate temporary field variable storage, set up the vectors to store variable and function data.
Timestep style is Crank-Nicholson by default, or you can specify otherwise using: -tsstyle implicit for fully-implicit or -tsstyle gear_2 for M=2 Gear style (see Jorge's reference).
Variable timesteps are of just one variety for now: exponential growth of deltat with timestep number up to a maximum value. This is controlled by the options -ts_dt_max, which sets the maximum value, and -ts_dt_factor, which is the amount by which to multiply dt at the beginning of each timestep (default 1.1).
The symmetry flag SYMMETRY_ZERO_INSIDE requires that the residual at the symmetry plane of nodes be set to the value of the field variable, so the solver will set that field variable to zero there. This in turns require that the field variable's tymestyle be either CONSTRAINT_ONLY or TIME_CONST_BLEND so the residual can be directly set at those nodes (done in function FuncEvaluate). Equations of type TIME_CONST_BLEND are given PetscTruth arrays indicating constraint status for all mixed fields over all points in the local part of the global array (i.e. not including ghost points). Temporary fields are given one large PetscScalar array for all temporary fields over all points (including ghost points) ordered by z, y, x, then (temporary or mixed) field variable. Both of these arrays are ordered by z, y, x, then in the latter case temporary field variable, which is similar to PETSc vector array orderings.
Now we're ready to create the various function vectors, and evaluate the time derivatives in this initial state.
Next we initialize the PETSc SNES object used to solve the equations in each timestep. If -snes_mf is specified, then we run matrix-free. Otherwise, we construct a Jacobian, either analytical if jack was passed non-null, or using PETSc's finite difference coloring method as used in SNES tutorial example 14. The finite difference part was put together based on PETSc example src/snes/utils/damgsnes.c and I really don't understand it...
Then we run the timestep loop, which has three parts:

Use PETSc's SNES solver to calculate the next timestep, with the correct timestyle coefficient for multiplying timestep equation functions and Jacobian rows (deltat for fully implicit, deltat/2 for Crank-Nicholson, 2*deltat/3 for M=2 Gear, etc.).
Check convergence and, if it failed and -noconv_cuttime is set, then back up with a smaller timestep and restore previous timestep as first guess for next.
Otherwise the timestep was successful, so produce the new next timestep guess vector and copy solution and function vectors into the old, update the time, and call the monitor.

The command-line option -guess_explicit turns on "explicit initial guess" for the next timestep, which extrapolates the difference between the previous and current timesteps to estimate the next one.

Global Function ts_implicit_function()

This is the SNES callback for the constrained implicit timestepper implicit_steptime.

int ts_implicit_function ( SNES thesnes, Vec unp1, Vec snesF, void* user )

int ts_implicit_function: It returns zero or error code.
SNES thesnes: SNES object in question.
Vec unp1: Field variable values.
Vec snesF: Where to return the function value.
void* user: User parameters' structure.

Global Function ts_implicit_jacobian()

This is the SNES Jacobian callback for the constrained implicit timestepper implicit_steptime.

int ts_implicit_jacobian ( SNES thesnes, Vec unp1, Mat* jack, Mat* preck, MatStructure* flag, void* user )

int ts_implicit_jacobian: It returns zero or error code.
SNES thesnes: SNES object in question.
Vec unp1
Mat* jack: Jacobian matrix.
Mat* preck: Preconditioner matrix.
MatStructure* flag: Flag indicating preconditioner structure.
void* user: User parameters' structure.

Vex unp1 Field variable values.

Local Function FuncEvaluate()

This evaluates the set of functions which make up the time derivatives and constraints to evaluate. It does its own localization of the global vector, etc. It assumes that temp_fields is already set up, as are all of the corner parameters (xs, ys, xm, gxs, etc.).

static int FuncEvaluate ( Vec u, Vec u_local, Vec Fu, Vec Fu_local, PetscScalar time, void* user )

int FuncEvaluate: It returns zero or error code.
Vec u: Field variable values.
Vec u_local: Local vector to copy the field variables into.
Vec Fu: Global vector to put the function values into.
Vec Fu_local: Local vector for temporary function value storage.
PetscScalar time: Current simulation time.
void* user: User parameters' structure.

Local Function xmax_symm()

This takes a local array which is oversized (since we told PETSc we're using periodic BCs) and replaces the outer edges with mirror reflections across the lines just past the outer columns. It has per-dof options to put the symmetry lines at the outer columns/planes and set certain dofs to zero at or just beyond the outer columns (in rheoplast these options are set in the modules).

static void xmax_symm ( PetscScalar* localarray, int xs, int ys, int zs, int xm, int ym, int zm, int gxs, int gys, int gzs, int gxm, int gym, int gzm, int nx, int* symmtypes, int dof )

PetscScalar* localarray: ``Local'' array with shadow nodes provided by the PETSc distributed array.
int xs: Starting global x-coordinate of the interior region.
int ys: Starting global y-coordinate of the interior region.
int zs: Starting global z-coordinate of the interior region.
int xm: Width of the interior region in the x-direction.
int ym: Width of the interior region in the y-direction.
int zm: Width of the interior region in the z-direction.
int gxs: Starting global x-coordinate of the local array.
int gys: Starting global y-coordinate of the local array.
int gzs: Starting global z-coordinate of the local array.
int gxm: Full width of the local array in the x-direction.
int gym: Full width of the local array in the y-direction.
int gzm: Full width of the local array in the z-direction.
int nx: Overall width of the entire distributed array in the x-direction.
int* symmtypes: Array of symmetry types for each dof (if NULL then assumes all are MIRROR_PLANE).
int dof: Number of degrees of freedom per node.

Local Function xmin_symm()

static void xmin_symm ( PetscScalar* localarray, int xs, int ys, int zs, int xm, int ym, int zm, int gxs, int gys, int gzs, int gxm, int gym, int gzm, int* symmtypes, int dof )

PetscScalar* localarray: ``Local'' array with shadow nodes provided by the PETSc distributed array.
int xs: Starting global x-coordinate of the interior region.
int ys: Starting global y-coordinate of the interior region.
int zs: Starting global z-coordinate of the interior region.
int xm: Width of the interior region in the x-direction.
int ym: Width of the interior region in the y-direction.
int zm: Width of the interior region in the z-direction.
int gxs: Starting global x-coordinate of the local array.
int gys: Starting global y-coordinate of the local array.
int gzs: Starting global z-coordinate of the local array.
int gxm: Full width of the local array in the x-direction.
int gym: Full width of the local array in the y-direction.
int gzm: Full width of the local array in the z-direction.
int* symmtypes: Array of symmetry types for each dof (if NULL then assumes all are MIRROR_PLANE).
int dof: Number of degrees of freedom per node.

Local Function ymax_symm()

static void ymax_symm ( PetscScalar* localarray, int xs, int ys, int zs, int xm, int ym, int zm, int gxs, int gys, int gzs, int gxm, int gym, int gzm, int ny, int* symmtypes, int dof )

PetscScalar* localarray: ``Local'' array with shadow nodes provided by the PETSc distributed array.
int xs: Starting global x-coordinate of the interior region.
int ys: Starting global y-coordinate of the interior region.
int zs: Starting global z-coordinate of the interior region.
int xm: Width of the interior region in the x-direction.
int ym: Width of the interior region in the y-direction.
int zm: Width of the interior region in the z-direction.
int gxs: Starting global x-coordinate of the local array.
int gys: Starting global y-coordinate of the local array.
int gzs: Starting global z-coordinate of the local array.
int gxm: Full width of the local array in the x-direction.
int gym: Full width of the local array in the y-direction.
int gzm: Full width of the local array in the z-direction.
int ny: Overall width of the entire distributed array in the y-direction.
int* symmtypes: Array of symmetry types for each dof (if NULL then assumes all are MIRROR_PLANE).
int dof: Number of degrees of freedom per node.

Local Function ymin_symm()

static void ymin_symm ( PetscScalar* localarray, int xs, int ys, int zs, int xm, int ym, int zm, int gxs, int gys, int gzs, int gxm, int gym, int gzm, int* symmtypes, int dof )

PetscScalar* localarray: ``Local'' array with shadow nodes provided by the PETSc distributed array.
int xs: Starting global x-coordinate of the interior region.
int ys: Starting global y-coordinate of the interior region.
int zs: Starting global z-coordinate of the interior region.
int xm: Width of the interior region in the x-direction.
int ym: Width of the interior region in the y-direction.
int zm: Width of the interior region in the z-direction.
int gxs: Starting global x-coordinate of the local array.
int gys: Starting global y-coordinate of the local array.
int gzs: Starting global z-coordinate of the local array.
int gxm: Full width of the local array in the x-direction.
int gym: Full width of the local array in the y-direction.
int gzm: Full width of the local array in the z-direction.
int* symmtypes: Array of symmetry types for each dof (if NULL then assumes all are MIRROR_PLANE).
int dof: Number of degrees of freedom per node.

Local Function zmax_symm()

static void zmax_symm ( PetscScalar* localarray, int xs, int ys, int zs, int xm, int ym, int zm, int gxs, int gys, int gzs, int gxm, int gym, int gzm, int nz, int* symmtypes, int dof )

PetscScalar* localarray: ``Local'' array with shadow nodes provided by the PETSc distributed array.
int xs: Starting global x-coordinate of the interior region.
int ys: Starting global y-coordinate of the interior region.
int zs: Starting global z-coordinate of the interior region.
int xm: Width of the interior region in the x-direction.
int ym: Width of the interior region in the y-direction.
int zm: Width of the interior region in the z-direction.
int gxs: Starting global x-coordinate of the local array.
int gys: Starting global y-coordinate of the local array.
int gzs: Starting global z-coordinate of the local array.
int gxm: Full width of the local array in the x-direction.
int gym: Full width of the local array in the y-direction.
int gzm: Full width of the local array in the z-direction.
int nz: Overall width of the entire distributed array in the z-direction.
int* symmtypes: Array of symmetry types for each dof (if NULL then assumes all are MIRROR_PLANE).
int dof: Number of degrees of freedom per node.

Local Function zmin_symm()

static void zmin_symm ( PetscScalar* localarray, int xs, int ys, int zs, int xm, int ym, int zm, int gxs, int gys, int gzs, int gxm, int gym, int gzm, int* symmtypes, int dof )

PetscScalar* localarray: ``Local'' array with shadow nodes provided by the PETSc distributed array.
int xs: Starting global x-coordinate of the interior region.
int ys: Starting global y-coordinate of the interior region.
int zs: Starting global z-coordinate of the interior region.
int xm: Width of the interior region in the x-direction.
int ym: Width of the interior region in the y-direction.
int zm: Width of the interior region in the z-direction.
int gxs: Starting global x-coordinate of the local array.
int gys: Starting global y-coordinate of the local array.
int gzs: Starting global z-coordinate of the local array.
int gxm: Full width of the local array in the x-direction.
int gym: Full width of the local array in the y-direction.
int gzm: Full width of the local array in the z-direction.
int* symmtypes: Array of symmetry types for each dof (if NULL then assumes all are MIRROR_PLANE).
int dof: Number of degrees of freedom per node.