
NAME
vemp02  solves a nonsteady nonlinear functional equation by mixed finite elements
SYNOPSIS
 CALL VEMP02(

T, LU, U, EEST, LIVEM, IVEM, LLVEM, LVEM, LRVEM, RVEM, LNEK, NEK,
LRPARM, RPARM, LIPARM, IPARM, LDNOD, DNOD, LRDPRM, RDPARM,
LIDPRM, IDPARM, LNODN, NODNUM, LNOD, NOD, LNOPRM, NOPARM,
LBIG, RBIG, IBIG, MASKL, MASKF, USERB, USRFUT, USRFU, USERF, VEM50X, VEM63X)
 INTEGER

LU, LIVEM, LLVEM, LRVEM, LNEK, LRPARM, LIPARM, LDNOD, LRDPRM,
LIDPRM, LNODN, LNOPRM, LBIG
 INTEGER

IVEM(LIVEM), NEK(LNEK), IPARM(LIPARM), DNOD(LDNOD),
IDPARM(LIDPRM), NODNUM(LNODN), IBIG(*)
 DOUBLE PRECISION

T, U(LU), EEST(LU), RVEM(LRVEM), RPARM(LRPARM), RDPARM(LRDPRM),
NOD(LDNOD), NOPARM(LNOPRM), RBIG(LBIG)
 LOGICAL

LVEM(LLVEM), MASKL(NK,NK,NGROUP), MASKF(NK,NGROUP)
 EXTERNAL

USERB, USRFUT, USRFU, USERF, VEM50X, VEM63X
PURPOSE
vemp02 is a subroutine for the numerical solution of a nonlinear nonsteady
functional equation on a space of smooth functions H. The problem for
the solution U: [T0,TEND] > H has to be given in the formulation:
[1] functional equation
for all V in H with V(COMPV)D(COMPV)=0 for all COMPV=1,..,NK:
F{T,UT,U}(V)=0 for all T0<T<=TEND
[2] Dirichlet conditions
for all COMPU=1,...,NK : U(T)(COMPU)D(COMPU)=B(T)(COMPU)
for all T0<T<=TEND
[3] Initial Solution
for all COMPU=1,...,NK : U(T0)(COMPU)=U0(COMPU)
where F is a linear form depending on the searched
solution U, its derivative UT with respect of time T and the time T,
see userf. B must not depend on U, see userb. U0 is a given
solution in H. Using the finite element method the linear form is
discretized. This transforms the problem
to a nonlinear initial value problem. It is solved by difference formulas
with selfadapted step size and order control. In every time step
the resulting system of nonlinear equations is solved
by the NewtonRaphson method. The quality of the computed
solution is estimated by an error estimation in space and time
direction. In every NewtonRaphson
step and for the calculation of the error indicator
the program package lsolpp is used.
At the first call of vemp02, U specifies the
initial solution U0 and EEST its error (normally =0). You can use
vemu08 or vemu06 to define U0. The initial solution
is modified by vemp02, so that it fullfils the Dirichlet conditions
at T0. Starting from T0,
the integration in time direction stops if the current time
is greater than the given time mark T. If INTERP is
set, the solution is evaluated at T, else the solution is given back
at the current time step. The computed solution and error estimator is
postprocessed by the routines vemu03, vemu04 and
vemu05 and can be handed over to a postprocessor program, see
vemide, vempat or vemisv. If at the output
STEP=1,
vemp02 can be recalled with a new time mark. If at the output
STEP=0, TEND is reached and the integration in Tdirection
ends.
ARGUMENTS
 T double precision, scalar, input/output, local

Time mark. If T is reached, vemp02 returns to the
calling program. At the output, T gives the current time value
at which the solution U is returned. So
if INTERP=LVEM(12)=.false., the value of T will be changed.
 LU integer, scalar, input, local

Length of solution vector U, LU >=LM.
 U double precision, array: U(LU), input/output, local

The solution vector at the global nodes at time T. U(i) is
the value at the global node i+PTRMBK(MYPROC), see
vemdis. On the first call U defines the values
of the initial solution at the initial time T0.
 EEST double precision, array: EEST(LU), input/output, local

The vector of the error indicator at the global nodes at time
T. EEST(i) is the value at the global node
i+PTRMBK(MYPROC), see
vemdis. On the first call EEST defines the error
of the initial solution at the initial time T0, normally =0. For
a mesh adaptation
the error indicator can be interpolated to the center point
of elements (see vemu03) or to the geometrical
nodes (see vemu05). If ERREST is not set, EEST
is undefined.
 LIVEM integer, scalar, input, local

Length of the integer information vector, LIVEM>=
MESH+ NINFO+100+50*NGROUP+ NDNOD+
NK+ LM.
 IVEM integer, array: IVEM(LIVEM), input/output, local/global

Integer information vector.
 (1)=MESH, input, local

Start address of the mesh informations in IVEM,
MESH>203+ NPROC.
 (2)=ERR, output, global

Error number.
0  program terminated without error. 
1  program stops, since prescribed CPU time is over. 
2  MAXIT is reached. 
3  the NewtonRaphson iteration does not converge. 
4  the NewtonRaphson correction is too small. 
10  error in lsolpp. 
80  illegal T. 
90  LBIG is too small. 
95  L/I/RVEM arrays or solution array is too small. 
98  read/write error. 
99  fatal error.   (3)=STEP, input/output, global

Recall indicator. At the input it indicates
0  first call 
1  continuation of Tintegration 
2  restart  and at the output it indicates
0  TEND is reached. 
1  call again to continue Tintegration 
2  CPUtime reaches JOBEND. Save arrays for restart. 
vemp02 offers the possibility to split the computation of the
numerical solution into several parts. This may be useful for very
difficult problems with a high computational amount per NewtonRaphson
step. The iteration stops if SECOND is greater than
JOBEND= RVEM(1), where SECOND is a real function
in the VECFEM library which gives the running CPU time
in seconds. SECOND is checked after every call of lsolpp. So
if you want to stop the iteration after RUNTIME CPU seconds, you have to set
JOBEND= RVEM(1)= SECOND()+RUNTIME before the call of
vemp02. If the current value of SECOND is greater than
JOBEND on any processor, vemp02 stops and sets STEP=
IVEM(3)=2. Then you have to save the mesh arrays
NEK, RPARM, IPARM, DNOD, RDPARM,
IDPARM, NODNUM, NOD, NOPARM, the
solution array U and the
information arrays IVEM, RVEM, LVEM
for every individual process (e.g. see vemp02exm10(7)). Before you
restart vemp02 you have to read the mesh arrays,
solution and information arrays (keep the TIDS in IVEM!)
(e.g. see vemp02exm10) and set STEP=2 to indicate a
restart. vemdis need not be recalled.
 (5)=NIVEM, output, local

Used length of IVEM.
 (6)=NRVEM, output, local

Used length of RVEM.
 (7)=NLVEM, output, local

Used length of LVEM.
 (8)=SPACE, output, local

Pointer to available work space in RBIG.
 (9)=LSPACE, output, local

Length of available work space in RBIG.
If the Tintegration will be continued after a return of vemp02 with
STEP=1, only the elements RBIG(SPACE),
RBIG(SPACE+1), ..., RBIG(SPACE+LSPACE1)
in RBIG may be changed. Especially this storage is only
available as work space for the postprocessing of the solution and
error indicator.
 (10), input, local

Unit for paging.
 (11), input, local

Unit for paging.
 (12), input, local

Unit for paging. If it is necessary, vemp02 writes parts (pages) of
RBIG/IBIG to external data sets. For that purpose vemp02
needs three temporary files on units IVEM(10), IVEM(11) and
IVEM(12). They are only used if the storage for RBIG/IBIG is
too small. The needed lengths of these files cannot be computed
in advance because they depend on the given mesh and the functional
equation. Every process has to get its own data set, otherwise the
results will be chaotic. If the data sets are allocated on different
discs, the i/otime will be reduced, since the i/o is done in parallel.
 (40)=LOUT, input, local

Unit number of the standard output file, normally 6.
 (41)=OUTCNT, input, local

Output control flag.
0  only error messages are printed. 
>0  a protocol and every OUTCNTth iteration step in lsolpp are printed.   (45)=MSPACK, input, global

Maximal number of stripes to pack the global matrix, normally 100. The
packing of the global matrix is divided into
several steps (called stripes). Before the packing
starts, the needed number of stripes is estimated. If
this number is greater than MSPACK, the computation will
be stopped. You have to give more storage for
RBIG/IBIG or to increase MSPACK. In general a problem
needing more than 100 stripes is too large for the
given storage, or else the mesh is numbered badly.
 (46)=PCLASS, input, local

Packing limit, normally 0. The global matrix is stored in packed form into
RBIG/IBIG. The needed storage can be controlled by
PCLASS.
0  lsolpp will need minimal CPU time. The needed storage is large. 
1  compromise of needed storage and CPU time for lsolpp. 
2  The storage for the global matrix will be minimal. lsolpp will need much CPU time.   (51)=ORDER, input, global

Order of the integration formulas for the computation of the
element matrices (0<ORDER<19). ORDER gives
the maximal degree of the polynomials which will be
integrated exactly. You should select ORDER greater than
the square of the order of the used proposal functions.
 (60)=MAXIT, input, global

Maximal number of NewtonRaphson steps. If MAXIT=0,
no limit for the number of iterations is set.
 (70)=MS, input, global

Method selection in lsolpp.
 (71)=MSPREC, input, global

Normalization method in lsolpp.
 (72)=ITMAX, input, global

Permitted maximal number of matrixvector multiplications per
NewtonRaphson step in lsolpp.
 (73)=MUNIT, input, global

If MUNIT is greater than 20 and less than 100 the global
matrix of the last LINSOL call is written to unit MUNIT,
see fBlsolpp.
 (200)=NPROC, input, global

Number of processes, see combgn.
 (201)=MYPROC, input, local

Logical process id number, see combgn.
 (202)=NMSG, input/output, global

Message counter. The difference of the input and the output values
gives the number of communications during the vemp02 call.
 (204)=TIDS(1), input, global

Begin of the list TIDS which defines the mapping of the
logical process ids to the physical process ids. See combgn.
 (MESH), input, local

Start of mesh informations, see mesh.
 LRVEM integer, scalar, input, local

Length of the real information vector,
LRVEM>=40+ 14*NK+ 8*LM, and if ERREST
is set, LRVEM>=40 +14*NK+ 15*LM.
 RVEM double precision, array: RVEM(LRVEM), input/output, local/global

Real information vector.
 (1)=JOBEND, output, local

If the real function SECOND in the VECFEM library is
greater than JOBEND on any processor, vemp02 stops
the calculation. You can save the mesh arrays, the solution vector and
the information array for a restart (see STEP). If JOBEND
is equal to zero, no limit for the run time is set.
 (2)=EPS, output, local

Smallest positive number with 1.+EPS>1.
 (3)=EPSEST, input, global

Desired accuracy for the solution of the linear systems to
compute the error indicator, e.g. 1.D2.
 (10)=TOL, input, global

Desired accuracy for the relative error of the solution, e.g. 1.D7. The
step size and order will be selected to keep the estimated error
lower than TOL.
 (11)=T0, input/output, global

At the first call, T0 specifies the initial time at which the initial
solution U0 is given. At the output it is set to the current time mark. So
the arrays U and EEST always contain the
solution and its estimated error at time RVEM(11).
 (12)=TEND, input, global

Integration in Tdirection ends if
the current time step is greater than TEND.
 (13)=H, input, global

The current step size. At the first call a start step size
has to be given. It should not be selected too large to avoid
failure of the first Tstep.
 (14)=HMIN, input, global

Minimal step size in Tdirection. If the current step
size is equal to HMIN, the error is accepted even if it is too
large. Setting HMIN to a sufficiently small positive value
avoids decrease of the current step size H down to zero.
 (15)=HMAX, input, global

Maximal step size in Tdirection. In rare cases, for example when the
solution is very smooth and almost constant, but also shows oscillations in
short intervals, it may be necessary to choose HMAX according to the
problem. In general, one sets HMAX = TEND  T0.
 LLVEM integer, scalar, input, local

Length of the logical information vector,
LLVEM>=20+2*NK+
5*NGROUP*(NK*NK+ NK).
 LVEM logical, array: LVEM(LLVEM), input/output, local/global

Logical information vector.
 (1)=LSYM, input, global

LSYM=true indicates a symmetrical Frechet derivative
of F with respect of U and UT, see usrfu.
 (4)=SMLLCR, input, global

If SMLLCR=true, a NewtonRaphson correction which is too small
will be accepted, in the other case the solution processes is stopped,
normally false. If the problem is very illposed a small NewtonRaphson
correction is not an indicator for a good solution, so
SMLLCR=true should only be set, if you are sure that your
problem is wellposed or you solve a test problem.
 (6)=ERRSTP, input, global

If ERRSTP=true, the estimation of
the discretization error is considered in the stopping
criterion of the NewtonRaphson iteration and the
step size control, normally true.
 (7)=ERREST, input, global

If ERREST=true, the global error estimation is computed,
else EEST is undefined, normally true.
 (8)=USESNI, input, global

If USESNI=true, the simplified NewtonRaphson method may be
used, normally true. The use of the simplified NewtonRaphson
method will reduce the CPU time, since the mounting of a new
global matrix is saved, but it might be risky in the case
of illposed problems.
 (9)=NOSMTH, input, global

NOSMTH=true indicates that lsolpp returns the smoothed
solution, normally false. For some applications the use of the
non smoothed solution can improve the convergency of the NewtonRaphson
iteration.
 (10)=LMVM, input, global

If LMVM=true, the NewtonRaphson iteration is continued even if
lsolpp reaches the maximal number of matrixvector
multiplications, normally true.
 (11)=NORMMA, input, global

If NORMMA=true, the consistency order, the
step size in the time direction and the NewtonRaphson iteration
are controlled by the maximum error over all components, else they are
controlled by the error of each individual component, normally false.
 (15)=ALLP, input, global

If ALLP=true, the consistency order in the time direction
is controlled in the range of 1 to p+1, else in the range of
p1 to p+1, where p is the current order, normally false.
 (16)=INTERP, input, global

If INTERP=true, the solution is interpolated at the given mark
T and is returned to the calling program, else the solution is
given at the current time step if it is greater than T
(T is set to the current time step), normally true.
 LNEK integer, scalar, input, local

Length of the element array.
 NEK integer, array: NEK(LNEK), input, local

Array of the elements, see mesh.
 LRPARM integer, scalar, input, local

Length of the real parameter array.
 RPARM double precision, array: RPARM(LRPARM), input, local

Real parameter array, see mesh.
 LIPARM integer, scalar, input, local

Length of the integer parameter array.
 IPARM integer, array: IPARM(LIPARM), input, local

Integer parameter array, see mesh.
 LDNOD integer, scalar, input, local

Length of the array of the Dirichlet nodes.
 DNOD integer, array: DNOD(LDNOD), input, local

Array of the Dirichlet nodes, see mesh.
 LRDPRM integer, scalar, input, local

Length of the real Dirichlet parameter array.
 RDPARM double precision, array: RDPARM(LRDPRM), input, local

Array of the real Dirichlet parameters, see mesh.
 LIDPRM integer, scalar, input, local

Length of the integer Dirichlet parameter array.
 IDPARM integer, array: IDPARM(LIDPRM), input, local

Array of the integer Dirichlet parameters, see mesh.
 LNODN integer, scalar, input, local

Length of the array of the id numbers of the geometrical nodes.
 NODNUM integer, array: NODNUM(LNODN), input, local

Array of the id numbers of the geometrical nodes, see mesh.
 LNOD integer, scalar, input, local

Length of the array of the coordinates of the geometrical nodes.
 NOD double precision, array: NOD(LNOD), input, local

Array of the coordinates of the geometrical nodes, see mesh.
 LNOPRM integer, scalar, input, local

Length of the array of the node parameters.
 NOPARM double precision, array: NOPARM(LNOPRM), input, local

Array of the node parameters, see mesh.
 LBIG integer, scalar, input, local

Length of the real work array. It is impossible to compute the
needed length for LBIG before the first vemp02
run. It depends on the given mesh and the
functional equation. The needed length of LBIG
is controlled by the parameters MSPACK and PCLASS. It
should be as large as possible.
 RBIG double precision, array: RBIG(LBIG), work array, local

Real work array. Between two calls of
vemp02 to continue the Tintegration, only the elements
RBIG(SPACE),
RBIG(SPACE+1), ..., RBIG(SPACE+LSPACE1)
in RBIG may be changed.
 IBIG integer, array: IBIG(*), work array, local

Integer work array, RBIG and IBIG have to be defined
by the EQUIVALENCE statement. Between two calls of
vemp02 to continue the Tintegration, only the elements
RBIG((SPACE1)*2+1),
RBIG((SPACE1)*2+2), ...,
RBIG((SPACE1)*2+LSPACE*2)
in RBIG may be changed (for single precision versions (e.g. on Cray
coamputers) replace '*2' by '*1').
 MASKL logical, array: MASKL(NK,NK,NGROUP), input, global

If MASKL(COMPV,COMPU,G)=true, the coefficient of the COMPVth
component of the test function in linear form F depends on the COMPUth
component of the solution or its derivative with respect to
time at the elements in the group G, see usrfu.
 MASKF logical, array: MASKF(NK,NGROUP), input, global

If MASKF(COMPV,G)=true, the COMPVth component
of the test function gives a contribution to the linear
form F over the elements in group G, see userf.
If you select the masks in an optimal manner, the computation will use the
CPU and storage resources optimally. If you are uncertain about the correct
settings, you can set the masks equal to true or call vemfre
to check your entries.
 USERB external, local

Name of the subroutine in which the
Dirichlet conditions are described, see userb.
 USRFUT external, local

Name of the subroutine in which the Frechet derivative of the linear
form F with respect to UT is described, see usrfu.
 USRFU external, local

Name of the subroutine in which the Frechet derivative of the linear
form F with respect to U is described, see usrfu.
 USERF external, local

Name of the subroutine in which the coefficients of the
linear forms F are described, see userf.
 VEM50X external, global

Name of the subroutine in which the element matrices are computed. The
general case is VEM500. Special versions for structural analysis
are not yet available.
 VEM60X external, global

Name of the subroutine in which the storage of VEM50X is
specified. The general case is VEM630.
EXAMPLE
See vemexamples.
METHOD
The functional equation [1][3] is discretized by the finite
element method. The finite element mesh is defined by the user. The
initial value problem is integrated by the finite difference
method. The solution of the system of nonlinear equations which has
to be solved in every time step T is an approximation of the
solution of [1][3] at T.
Spatial Discretization
In every element the solution is replaced by a polynomial which
interpolates the solution at the global nodes of the element. The
parametric representation is the polynomial interpolation of the
assigned geometrical nodes.
Time Discretization
From earlier time steps the solution at the current time is
computed by a difference formula of the order of consistency P and
step size H. P and H are controlled by an estimation of
the error on the level of the
equation so that the error of the solution is lower
than the given tolerance TOL. If the current time reaches the given
time mark T, the integration stops to save the solution at time T.
Solution at a Time Step
In every time step, the discretizations in space and time direction
produce a system of nonlinear equations. It is solved by the
Newton method. The iteration stops and the solution is accepted if
the space discretization error is reached or if it converges nearly
quadratically and the correction of the next iteration will be
smaller than TOL. It may be necessary to use underrelaxation. For
good convergence the simplified Newton method can be used.
Solution of a Newton Step
Systems of linear equations have to be solved which are in general
very large and extremely sparse. The linear equations solver
lsolpp in vemp02 uses iterative methods.
Error Estimation
The discretization error in space direction is estimated by
inserting other test functions into the functional equation as they
are used for the computation of the solution. The functions are
generated by halving the mesh size and the order. The
discretization error in time direction is estimated by the
difference of the derivative of the interpolation polynomials of the
current and of the next higher order. To compute the error function
EEST the global matrix of the last Newton step is used.
Accuracy
The accuracy is determined by the error tolerance TOL which
controls the error of the solution of the discretized functional
equation. Additionally the accuracy of the numerical solution
depends on the mesh width and the order of the proposal functions,
which is controlled by the error estimation, the order of
the geometrical approximation and the used order of integration formulas,
which is not considered in the estimation.
RESTRICTIONS

The existence and uniqueness of the solution of the problem have to be ensured.

lsolpp uses iterative methods. In
some cases problems with the convergence can be avoided if
the coefficients of the linear form are ordered in a way that
the derivatives of the COMPVth component of the solution occur
in the coefficients of the derivatives of the COMPVth component
from the test functions and/or the COMPVth component of the
solution occurs in the coefficients of the COMPVth component from
the test functions.

The NewtonRaphson iteration or the lsolpp iteration may be
divergent. In general this occurs if the coefficients for the derivatives of
the test functions get a dominant term of the solution relative
to terms of the derivative of the solution or the coefficients
for the test functions get a dominant term of the derivative of
the solution relative to terms of the solution. If it is
possible, scale the coefficients in the linear form so that
they are in the same order of magnitude. In case of divergence you
should check the Frechet derivatives and the mesh. Sometimes
an increase of the integration order will produce convergence.

An initial solution which equals zero may cause problems during
the iteration. In the beginning, the computation of the
tolerance on the level of equation may be unsafe because of the
small value of the norm of the solution. In most cases the
solution and the error estimation become more accurate when
more integration steps have been executed.

If the error estimation is active, the order of the used proposal
functions for a component must be greater than one, if the
derivative of this component of the test function or the solution
is used in the functional equation. In the other case, order one
is sufficient.
REFERENCES
[FAQ], [THEOMAN], [DATAMAN], [DATAMAN2], [LINSOL], [P_MPI],
SEE ALSO
VECFEM, vemcompile, vemrun, vemhint,
mesh, vemexamples, vemdis, lsolpp,
xvem, equation, userb, userf,
usrfu, vembuild, vemdis, vemfre,
vemge2, vemgen(later), vemopt(later), vemu06.
COPYRIGHTS
Program by L. Grosz, C. Roll, P. Sternecker, 19891996.
Copyrights by Universitaet Karlsruhe 19891996.
Copyrights by Lutz Grosz 1996.
All rights reserved. More details see VECFEM.
by L. Grosz, Auckland , 6. June, 2000.
