
/home/meimei/strain/src/pf/README


Overview
========

This directory contains the source code for the nonlinear parabolic
spectral code pfb.f . The method is described in the papers
``Spectral Methods for Nonlinear Parabolic Systems,'' Jour. Comp. Phys.

The code is written in Fortran 77 and uses single precision arithmetic.
However, the "-r8" flag of Sun f77 is used to automatically convert to double
precision. The makefile compiles the code for fast execution on the machine
where it is compiled, using the "-fast -native -O5" flags of Sun f77.
The libraries linked should probably all be compiled with the same flags set.

The main driver is pfb.f, which writes selected results to user-specified files.
Makefiles are used extensively to automate result production.

This directory contains solvers for nonlinear parabolic systems with
periodic boundary conditions in 2-D geometry, with particular emphasis 
on phase field models for solidification. 

The solvers are based on implicit ODE-type time discretizations of the 
parabolic system, yielding linear elliptic systems at each time step.

The elliptic systems are solved by analytic preconditioning: the averaged 
operator is used as a preconditioner and the linear systems are solved by 
GMRES with a spectral discretization of the resulting integral operator.

Conventions
===========

Source code files 	<filename>.f
Object code files 	<filename>.o
Executables       	<filename>.x     This is nonstandard but convenient.
Input files 		<filename>.in
Directories of results	out/<example name>.YYMMddhhmmss where YY... is the date.
Library files in directory	$HOME/lib
This directory			$HOME/src/pf

Linkages
========

The code calls the following libraries of auxiliary routines from $HOME/lib:

dps	PostScript graphics routines.
blas	BLAS.
fft	FFT by Swarztrauber.
misc	Miscellaneous subroutines.
vce	Variable-coefficient elliptic solvers with analytic preconditioning.

Usage
=====

The drivers are compiled with 
> m pfb.x
and run on a given numerical example by entering
> m wbm
where m is an alias for /usr/ccs/bin/make DATE=`/usr/bin/date +%y%m%d%H%M%S`
and the available examples currently include
	lq 	Single variable-coefficient linear equation.
	ld 	Diagonal linear system.
	ls 	Nondiagonal linear system with diagonal principal part.
	mq 	Mean curvature equation.
	md 	Mean curvature diagonal system: each component satisfies
	ms 	Mean curvature system with lower-order coupling.
	zq 	Linear equation plus zero-order cubic nonlinearity.
	zd 	Linear diagonal system plus zero-order cubic nonlinearity.
	zs 	Linear system plus zero-order cubic nonlinearity.
	cf 	Caginalp-Fife phase field model.
	wbm 	Wheeler-Boettinger-McFadden alloy phase field model.
This command runs the example in a subdirectory of ./out named wbm.YYMMddhhmmss.

Input parameters are read from a file with suffix .in, for example wbm.in:

***********************Input file for pfb.f************************************
  t1 	t2
  0.0	01.0
  ieqn 		nq 	neqn 	mt	mx	my		Eqn parameters
  11		2 	6	1 	1	1
  am 	ee 	wa 	wb 	aa 	ab  			peqn 
  50.0 	0.02	10.0 	10.0 	0.40 	0.60
  isol 		nsol	kt	kx	ky			Soln parameters
  1		8	1	5	5
  pmin	pmax	pshx	pshy	umn	umx	ushx	ushy	psol 
  0.49	0.51	0.10	0.20	0.49	0.51	0.30	0.40
  nt		iode		kmin		kmax	nx	Numerical stuff
  100000	4		1		6	256
  dti		dtr	dtf	dtmn	attol	rttol	faclte	factol
  1.0E-06	2.8	1.0E-02	1.0E-10	1.0E-03	1.0E-03	0.7	0.1
  actol		rctol	artol	rrtol	stall	facstl	newtmx	facnon
  1.0E-10 	1.0E-10 1.0E-10 1.0E-10	0.3	0.4	06.	0.2
  idim	iortho	irstrt	irhs	gtoli	gtolr	gtolf	facnog
  10 	10 	10	1	1.0E-10	1.0E-01	1.0E-10	0.5
  nwrt 		dtwrt 	 	 	 	 		I/O parameters
  1 		0.100
  Output files		iu iq i1 i2 n1   inc
  'standard_out'	6  1  0  0  0    0
  'wbm.tst1.out'	7  1  0  0  0    0
  Step/order plots	iu iq i1 i2 n1   inc 	ipz  bmn bmx nxt nyt icp
  'stepsizes.ps'	8  1  9  9  1    1	30   0   15  5   5   1
  'bdforders.ps'	9  1  9  9  1    1	30   0   12  5   6   1
  Grayscale plots	iu iq i1 i2 n1   inc 	ibx  icp ipz 
  'gry1.1000.ps'	10 1  5  8  1000 0 	1    1   30
  'gry2.1000.ps'	11 1  5  8  1000 0	1    1   30
  Contour plots		iu iq i1 i2 n1   inc 	ibx  icp ipz ncon
  'con1.1000.ps'	0  1  5  8  1000 0	1    1   30  10
  'con2.1000.ps'	0  1  5  8  1000 0	1    1   30  10
  Mean value plots	iu iq i1 i2 n1   inc 	ipz  bmn bmx nxt nyt icp
  'meanvalue.ps'	30 1  9  9  1    1 	30   0   1   5   4   1
  'meanvalue.ps'	31 1  9  9  2    1	30   0   1   5   4   1
  Error log plots	iu iq i1 i2 n1   inc 	ipos ipz bmn bmx nxt nyt icp
  'err1.bx.1.ps'	40 1  9  9  1    1	7    30  0   40  4   10  1
  'err2.bx.1.ps'	41 1  9  9  1    1 	7    30  0   40  4   10  1

Input parameters: Details
=========================

t1,t2:	The run starts at t=t1 and stops at t=t2.

ieqn:	The following equations are supported: The value of ieqn is given,
	then the two-letter suffix used in the program, then a description.
	All coefficients a_j, a_j, a_ijk etc. are given by random Fourier
	series evaluated by function as, with mt by mx by my terms.
	The coefficients are constant (but random) if mt=mx=my=0.
	Uniform ellipticity is however guaranteed.
	Parameter nqmx must be set to nq in both pf*.f and as.f.
	Notation: d_t and d_j are derivatives from the following table:
	j	d_j u
	1	u_xx
	2	u_xy
	3	u_yy
	4	u_x
	5	u_y
	6	u
	The inhomogeneous term f(x,y,t) is zero unless isol<0, in which
	case it is calculated to make function u=fs the exact solution:
	f=u_t-F(tx,y,u,du,ddu).

	01. (lq) Single variable-coefficient linear equation.
	d_t u =sum_j=1^6 a_j d_j u + f(x,y,t)

	02. (ld) Diagonal linear system.
	d_t u_i =sum_j=1^6 a_ij d_j u_i + f_i, 1 \le i \le nq,
	with d_j and a_ij as in 1.

	03. (ls) Nondiagonal linear system with diagonal principal part.
	d_t u_i =sum_j=1^3 a_ij d_j u_i 
		+ sum_j=4^6 sum_k=1^nq a_ijk d_j u_k + f_i

	04. (mq) Mean curvature equation.
	u_t = (1+u_y^2)u_xx + 2 u_x u_y u_xy + (1+u_x^2)u_yy 
		+ a_1 u_x + a_2 u_y + a_3 (u+u^3) + f(x,y)
	where a_i are given by random Fourier series.

	05. (md) Mean curvature diagonal system: each component satisfies
	u_t = (1+u_y^2)u_xx + 2 u_x u_y u_xy + (1+u_x^2)u_yy 
		+ a_1 u_x + a_2 u_y + a_3 (u+u^3) + f(x,y)
	where a_i are given by random Fourier series.

	06. (ms) Mean curvature system with lower-order coupling.

	07. (zq) Linear equation plus zero-order cubic nonlinearity.

	08. (zd) Linear diagonal system plus zero-order cubic nonlinearity.

	09. (zs) Linear system plus zero-order cubic nonlinearity.

	10. (cf) Caginalp-Fife phase field model.
	u_t = a_1 \Delta u + a_2 (u-u^3) + a_3 v
	v_t = - b1 \Delta u + b2 \Delta v - b3 (u-u^3)  - b4 v
	where 
	a1=peqn(1)
	a2=peqn(2)
	a3=peqn(3)
	b1=peqn(4)
	b2=peqn(5)
	b3=peqn(6)
	b4=peqn(7).

	11. (pf) Wheeler-Boettinger-McFadden alloy phase field model.
	ut = am (ee^2 \Delta u - fa v - fb (1-v) )
	vt = \Delta v 
		+ v(1-v) (fb-fa) \Delta u 
		+ (ux vx + uy vy) (1-2v) (fb-fa)
		+ (ux ux + uy uy) v(1-v) (fb' - fa')
	where 
	am=peqn(1)
	ee=peqn(2)
	wa=peqn(3)
	wb=peqn(4)
	aa=peqn(5)
	ab=peqn(6)
	and 
	fa=wa*u*(u-1)*(u-aa)
	fb=wb*u*(u-1)*(u-ab)

nq:	The number of components in the system or equation.
	Parameter nqmx must be set to nq in both pf*.f and as.f.
	Parameter nq also affects the structure of the input file as follows:
	the number of solution parameters nsol should be equal to 4*nq
	if either the exact solution or the random Fourier series is used.
	The number of graphics files should be equal to 2*nq since we
	must specify a contour plot and a grayscale plot for each component.

mt:	Number of Fourier modes used in the time variable for coefficients.
mx:	Number of Fourier modes used in the x variable for coefficients.
my:	Number of Fourier modes used in the y variable for coefficients.

neqn:	Number of parameters needed to specify the coefficients. Varies
	from equation to equation; see ieqn above for details.

peqn:	Parameters for coefficients. See ieqn above for details.

isol:	Solution selection parameter. If isol<0, the exact solution is 
	provided by function fs; otherwise, only initial values are known.
	 The following solutions and initial values are available:

	-1: Selects exact solution equal to exponential of three cosines
	with parameters kt,kx and ky:
	f(x,y,t)=pmn+(pmx-pmn)*(g(x-xsh,y-ysh,t)-gmn)/(gmx-gmn)
	where pmn,pmx,xsh and ysh are given by psol(4*(jq-1)+1,2,3,4) 
	respectively for the jq'th component of f and g is given by
	g(x,y,t)=exp((a+b*cos(2pi*kt*t))*(a+b*cos(2pi*kx*x))*(a+b*cos(2pi*ky*y))
	with a=0.75, b=0.25 set so that the vanishing of one cosine does
	not freeze the other variables. Here gmn and gmx are the min and max
	of g, so component jq of f lies between pmn and pmx for that jq.
	In this case, the RHS of the PDE u_t=F(x,y,t,u,du,ddu)+G(x,y,t)
	is altered to make f the exact solution, by adding 
	G(x,y,t)=f_t-F(x,y,t,f,df,ddf) as an inhomogeneous term.

	0: Same as -1 but the RHS of the PDE is not altered to make fs
	the exact solution.

	1: Pseudorandom grid values, independently chosen from a uniform
	linear congruential generator from Numerical Recipes. Values of 
	component number iq are scaled to lie between pmn=psol(4*(iq-1)+1)
	and pmx=psol(4*(iq-1)+2).

	2: Fourier sine and cosine series with random coefficients:
	u_i(x,y)=pmn+(pmx-pmn)*(g_i(x,y)-gmn)/(gmx-gmn)
	where
	g_i(x,y)=sum_ix=0^kx sum_iy=0^ky 
	 (a_ixiyi cos(2pi ix x) +b_ixiyi sin(2pi ix x))
	 (c_ixiyi cos(2pi iy y) +d_ixiyi sin(2pi iy y))
	where pmn and pmx come from psol as usual.

	3: Regular array of Gaussians. Not in use.

	4: Random array of k=kx*ky Gaussians, scaled to minimize overlapping.
	u_i(x,y)=pmn+(pmx-pmn)*(v_i(x-xsh,y-ysh)-vmn)/(vmx-vmn)
	where
	v(x,y)=sum_j=1^k g(x-x_j,y-y_j)
	with x_j and y_j pseudorandom uniform on [0,1].
	Scale ss of Gaussians is set so that k Gaussians can fit in [0,1]^2
	with no more than 1% overlap.  Thus
	ss is 1/R^2 where R is the scale length appearing in the Gaussian
	g(x,y)=exp(-(x^2+y^2)/R^2).
	The radius where the Gaussian with scale length R is less than eps
	is rho=sqrt(-log(eps))*R. To fit k of these circles in the box [0,1]^2,
	we need k*pi*rho^2 = fraction (to leave 1-fraction of the box empty).
	Thus R^2=rho^2/(-log(eps)) = 1/2*k*pi*(-log(eps)).
	Currently these parameters are hardwired as eps=0.01 and fraction=0.8.
	As usual pmn,pmx,xsh,ysh are fouind in psol(4*(iq-1)+1,2,3,4).

kt:	Number of Fourier modes used in the time variable for coefficients.
kx:	Number of Fourier modes used in the x variable for coefficients.
ky:	Number of Fourier modes used in the y variable for coefficients.

nsol:	Number of parameters in psol needed to specify the solution. Varies
	from solution to solution; see isol above for details.

psol:	Parameters for solution. See isol above for details.

nt:	Maximum allowed number of time steps to go from t1 to t2.

dti:	Minimum value of the time step for BDF; should be very small.

dtr:	Time step increase ratio for BDF: should decrease as iord increases.

dtf:	Maximum value of the time step for BDF; should not be very small.

dtol:	Tolerances for the time integration; meaningless unless iode>2.

iode:	Type of ODE solver to use. Currently supported:
	iode=0	Modified linearly-implicit Euler.
	iode=1	Extrapolated modified linearly-implicit Euler.
	iode=2	BDF/Newton corrector with polynomial extrapolation predictor.
	iode=3	BDF/Newton with adaptive step control 
		and modified divided differences.

iord:	Order of ODE solver. Should not exceed parameter 'km' in pfb.

ptol(1-5):Tolerances for Newton iteration. Here
	redmin=ptol(1), actol=ptol(2), rctol=ptol(3), artol=ptol(4),
	and rrtol=ptol(5).
	The 'a' is for absolute, 'r' for relative. The 'c' is for correction, 
	'r' for residual.
	Newton iteration is terminated when the L1-norms cnrm, rnrm, unrm, fnrm
	of the correction c, the residual r, the solution w, and the initial
	residual f satisfy
		cnrm < ac+rc*wnrm and rnrm < ar+rr*fnrm.
	However, no more than newtmx=10 Newton steps are permitted.
	Newton iteration is terminated immediately, the time step dt
	is reduced by a factor of 2, and the step is restarted if 
	the residual norm does not decrease by a factor redmin at any step.

nx:	Number of grid points per side of the domain [0,1]^2.

irhs	How to start GMRES iteration: The initial density s is set equal to
	s = random if irhs <1 (usually a terrible way to start).
	s = 0 if irhs =1.
	s = rhs of eqn to be solved if irhs=2.

idim:	Max dimension of Krylov space in GMRES.

iortho:	Number of previous vectors to orthogonalize against.

irestart:	How often to restart GMRES.

gtoli,gtolr,
gtolf:	Tolerances for stopping GMRES. Since early Newton steps live
	far from the solution, it is not necessary to solve very accurately
	for the solution to the linear system for the correction.
	Thus the first Newton step uses GMRES with tolerance tol=gtoli,
	and each succeeding step multiplies the tolerance tol by gtolr
	until tol<gtolf, whereupon tol is set to gtolf for the rest of
	the Newton steps.

nwrt:	Printing step control.  Output takes place at multiples of 
	a single basic step, determined by dtwrt and nwrt as follows:
	If nwrt>0, then the a basic output step is nwrt computational
	steps, typically of variable size. This is useful because the
	time step adjusts to follow the solution, so small steps are
	used when the solution is changing rapidly and large steps are
	used when the solution is almost constant. 
	If nwrt<=0, then the following parameter dtwrt is the basic 
	time step size and output is done at steps which cross a multiple
	of dtwrt.

dtwrt:	Basic output interval. Output to each file is done only
	at steps which cross multiples of dtwrt if nwrt is nonnegative;
	if nwrt is positive, dtwrt is ignored.
	Later should really arrange output precisely at dtwrt steps.
	Interpolation is easy with BDF; just another call to extrapolation
	routine. 

Output control parameters:
	pfb.f reads control parameters for several different types of output:

	First, concise human-readable output to the terminal and/or a file.
	Two lines are read in the following format:
	  'irrelevance'	6  1  0  0  000  0	fn iu iq i1 i2 n1 inc
	  'pf.test.out'	7  1  0  0  000  0	fn iu iq i1 i2 n1 inc
 	where

fn:	Basic file name for output or graphics. 11 characters.
	Irrelevant if the unit number is 6 indicating standard output.

iu:	Unit number for output or graphics; must be unique and 
	between 6 and 99 inclusive, with 6 reserved for terminal output.

iq:	Frequency of output in multiples of dtwrt. Setting iq=0 gives no output.
	Setting iq<0 is significant only for error plots, as discussed below.

i1,i2:	Beginning and end positions of number n1 in basic file name fn.
	If inc>0, a new file is opened at each printing step,
	with n1 embedded in its name, and closed after one frame 
	is written. If inc=0, the basic file name is used and
	successive frames are written to the same file.
	If the named file already exists, it is opened at the end
	and the new output is appended. This is useful for keeping
	a log of successive runs.
	Setting inc<0 is significant only for error plots, as discussed below.

n1:	Initial number for graphics files. Used only if inc>0.
	Incremented by inc at each frame.

inc:	Increment for file number n1.

	Second, solution plots in grayscale format to postscript files.
	A total of 2*nq lines are read in a format differing from the 
	above only as follows: The first line parameters have exactly the
	same meanings as for concise output. The second line for each 
	of the nq components reads:
  		1  1  30		ibox icap ipsz
	where

ibox:	>0 to make box around the plot, <=0 to have no box.

icap:	>0 to write a brief caption containing date, computational time,
	index, max and min of the component, etc. to the plot.

ipsz:	Point size to write captions in plot file.

	Third, solution plots in contour format to postscript files.
	A total of 2*nq lines are read in a format the same as for the
	grayscale plots except that the second line reads:
  		1  1  30 12		ibox icap ipsz ncon
	where

ncon:	Number of contour lines to plot between max and min of the function.

	Fourth, error plots when the exact solution is known.
	The quantity plotted is the truncated number of correct bits
	in each component of the solution, as discussed below.
	A total of 2*nq lines are read in format: 
  	'pf2.bx.1.ps'	30 1  8  8  1 0		fn iu iq i1 i2 n1 inc
  		6 30 0 30 4 15		ipos ipsz bmin bmax nxt nyt
	where all parameters mean the same as above except:

iu:	We use three file numbers for each component, and they must be kept
	open throughout the whole calculation, so the unit number iu should
	be spaced out. We recommend spacing them 10 apart as in the above
	example. The unit numbers actually used are iu,iu+1,iu+2 (one for 
	each norm plotted; L1, L2 and Lm (infinity) norm).

iq:	Since we no longer need a frequency (the error plots are done
	after the calculation is over), iq is used as a switch to indicate
	whether the plots are to appear all in the same file or one
	to a file, as follows
	If iq>0, each plot is written to a new file and 
	inc has its usual meaning.
	If iq=0, no plots are done at all.
	If iq<0, then inc takes on special meaning:
		inc>0 means initialize and plot,
		inc=0 means plot (initialization has been done),
		inc<0 means plot and terminate the plot.

ipos:	Location in the filename to put the norm indicator 1,2,m.
	Three plots are generated for each component, one for each of
	the L1, L2 and max-norms of the error, and each occupies a 
	separate file. Thus if the second input file name is 'pfn.bx.1.ps'
	and ipos=6 then the correct bits in u2 are plotted in files
	'pf2.b1.1.ps','pf2.b2.1.ps' and 'pf2.bm.1.ps' The 2 is added 
	to indicate the component number, while the 'b' stands for bits.

bmin:	Lower edge of the correct bits plot. 

bmax:	Upper edge of the correct bits plot. 

nxt:	Number of tick marks along the x-axis of the plot.

nyt:	Number of tick marks along the y-axis of the plot.

Fortran conventions
===================

General statement numbers are spaced at least 10 apart: 100, 110, ...
Format statements have four-digit statement numbers 9xyz.

