
content="text/html; charset=iso88591">
content="C:\PROGRAM FILES\MICROSOFT OFFICE\OFFICE\html.dot">
VECFEM3 Reference Manual: veme02
Type: FORTRAN routine
NAME
veme02  solves a nonlinear functional equation by mixed
finite elements
SYNOPSIS
 CALL VEME02(
 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,
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, USRFU, USERF, VEM50X, VEM63X
PURPOSE
veme02 is a subroutine for the numerical solution of a nonlinear
steady functional equation on a space of smooth functions H. The
problem for the solution U in H has to be given in the form :
[1] functional equation
for all V in H with V(COMPV)D(COMPV)=0 for all COMPV=1,..,NK:
F{U}(V)=0
[2] Dirichlet conditions
for all COMPU=1,...,NK : U(COMPU)D(COMPU)=B(COMPU)
where F is a linear form depending on the searched solution U,
see userf. B must not depend on U, see userb.
Using the finite element method the linear form is discretized.
The resulting system of nonlinear equations is solved by the
NewtonRaphson method. An error indicator estimates the quality
of the computed solution. In every NewtonRaphson step and for
the calculation of the error indicator the program package lsolpp is used. The computed solution has
to be postprocessed by the routines vemu03,
vemu04 and vemu05
and can be handed over to a postprocessor program, see vemide or vempat.
ARGUMENTS
 T double precision, scalar, input, local
 Real number (e.g. the current time).
 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. U(i)
is the value at the global node i+PTRMBK(MYPROC),
see vemdis. If STARTU
is set, U gives an initial guess for the
NewtonRaphson iteration.
 EEST double precision, array: EEST(LU),
output, local
 The vector of the error indicator at the global nodes. EEST(i)
is the value at the global node i+PTRMBK(MYPROC),
see vemdis. For a mesh
adaptation the error indicator can be interpolated to the
centre 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+36*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 information 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. 
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

 Restart indicator.
0 
the first call of veme02. 
2 
the restart of veme02. 
veme02 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 veme02. If
the current value of SECOND is greater than JOBEND
on any processor, veme02 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 veme02exm10).
Before you restart veme02 you have to read the
mesh arrays, solution and information arrays (keep
the TIDS in IVEM!) (e.g.
see veme02exm10)
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.
 (10), input, local
 Unit for paging.
 (11), input, local
 Unit for paging.
 (12), input, local
 Unit for paging. If it is necessary, veme02
writes parts (pages) of RBIG/IBIG
to external data sets. For that purpose veme02
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
can not 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 lsolpp.
 (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 veme02 call.
 (204)=TIDS(1), input, global
 Begin of the list TIDS that defines
the mapping of the logical process ids to the
physical process ids. See combgn.
 (MESH), input, local
 Start of mesh information, see mesh.
 LRVEM integer, scalar, input, local
 Length of the real information vector, LRVEM>=40+
14*NK+ 2*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, veme02 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.
 LLVEM integer, scalar, input, local
 Length of the logical information vector, LLVEM>=20+2*NK+
4*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, 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.
 (5)=STARTU, input, global
 STARTU=true indicates that U
specified an initial guess for the NewtonRaphson
iteration. Use vemu08 or vemu06 to set an initial guess.
 (6)=ERRSTP, input, global
 If ERRSTP=true, the estimation of the
discretization error is considered in the stopping
criterion of the NewtonRaphson iteration, normally true.
The iteration ends if the estimated discretization is
reached, so that no computing time is wasted to reach a
too strict TOL.
 (7)=ERREST, input, global
 If ERREST=true, the error indicator 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 convergence
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 NewtonRaphson iteration
is controlled by the maximum error over all components,
else they are controlled by the error of each individual
component, normally false.
 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 veme02 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.
 IBIG integer, array: IBIG(*), work
array, local
 Integer work array, RBIG and IBIG
have to be defined by the EQUIVALENCE statement.
 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 the
linear form F depends on the COMPUth component of the
solution 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.
 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.
 VEM63X 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][2] is discretized by the finite
element method. The solution of the resulting system of linear
equations is an approximation of the solution of [1][2].
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.
Solution
The discretization produces a system of nonlinear equations.
It is solved by the NewtonRaphson 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 NewtonRaphson method can be used.
Solution
Systems of linear equations have to be solved, which are in
general very large and extremely sparse. The linear equations
solver lsolpp in veme02 uses iterative
methods.
Estimation
The discretization error is estimated by inserting other test
functions into the functional equation as they are used for the
computation of the solution. Halving the mesh size and the order
generates the functions. To compute the error function EEST
the global matrix of the last NewtonRaphson 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.
 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
VECFEM, vemcompile,
vemrun, vemhint,
mesh, vemexamples,
vemdis, lsolpp,
xvem, equation,
userb, userf, usrfu, vembuild,
vemfre, vemge2,
vemgen(later), vemopt(later), vemu06, vemu08.
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, 4 June, 2000.
