C******************************************************************* C** C** v e m e 0 2 e x m 1 0 : C** C** System of three nonlinear patrial differential equations of C** the Navier-Stokes type on a 3-dimensional domain. The examples C** demonstrates the usage of restarts. The mesh is read from a C** VECFEM input file. C** C** by L. Grosz Karlsruhe, Mar. 1995 C** C** C******************************************************************* C** C** The problem is a system of three partial differential C** equations, which are a simplification of the Navier-Stokes C** equations (it is assummed that the pressure is constant). C** The domain is the 3-dimensional cube with length pi. An all C** boundaries Neuman boundary conditions are prescribed and for C** all components a Dirichlet condition at the origin is set. C** Using the notations in equation the problem is given by C** the functional equation: C** C** Dirichlet conditions: C** u1=0 C** u2=0 C** u3=0 C** C** linear functional equation: F{u}(v)=0 C** C** with F{u}(v)= C** volume{v1x1*u1x1+v1x2*u1x2+v1x3*u1x3+ C** Re*v1*(u1*u1x1+u2*u1x2+u3*u1x3-f1) C** + v2x1*u2x1+v2x2*u2x2+v2x3*u2x3+ C** Re*v2*(u1*u2x1+u2*u2x2+u3*u2x3-f2) C** + v3x1*u3x1+v3x2*u3x2+v3x3*u3x3+ C** Re*v3*(u1*u3x1+u2*u3x2+u3*u3x3-f3)} C** + area{v1*g1+v2*g2+v3*g3} C** C** The functions f and g are selected so that C** C** u1=cos(x1)*sin(x2)*sin(x3) C** u2=sin(x1)*cos(x2)*sin(x3) C** u3=sin(x1)*sin(x2)*cos(x3) C** C** is the exact solution of this problem. The real value Re C** controls numerical behaviour of the problem. For increasing C** values of Re the problem will get ill-posed and you have to C** expect that the calculation fails. C** C** The mesh is read from the vecfem input file wurf.vem. C** It uses hexahedron elements of order 2 for the subdivision C** of the cube and quadrilateral elements for its boundary. C** C**----------------------------------------------------------------- C** PROGRAM VEMEXM C** C**----------------------------------------------------------------- C** IMPLICIT NONE include 'bytes.h' C** C**----------------------------------------------------------------- C** C** some parameters which may be chanced: C** C** MAXNG = maximal number of groups C** INPUT = name of the vecfem input file C** STORE = total storage of process in Mbytes. C** INTEGER NPROC,MAXNG,STORE CHARACTER*80 INPUT PARAMETER (MAXNG=5, & INPUT='wurf.vem', & STORE=65) C** C**----------------------------------------------------------------- C** C** the parameters are explained in mesh. C** INTEGER NK,DIM,MESH,LOUT PARAMETER (NK=3, & DIM=3, & MESH=500, & LOUT=6) C** C**----------------------------------------------------------------- C** C** the length of the array for the mesh are set: C** INTEGER LU,LNODN,LNOD,LNOPRM,LNEK,LRPARM,LIPARM, & LDNOD,LIDPRM,LRDPRM,LIVEM,LRVEM,LLVEM,LBIG PARAMETER (LU =30 000, & LNODN =10 000, & LNOD =30 000, & LNOPRM=1, & LNEK=64 000, & LIPARM=6800, & LRPARM=1, & LDNOD =6, & LIDPRM=3, & LRDPRM=3, & LIVEM =MESH+800+LU+LDNOD/2, & LLVEM =500, & LRVEM =60+2*LU) C** C**----------------------------------------------------------------- C** C** RBIG should be as large as possible: the available C** storage STORE is reduced by all allocated array. C** the remaining storage is reserved for RBIG. C** PARAMETER ( LBIG=(STORE * 1 000 000)/IREAL & - (3*LU+LNOD+LNOPRM+LRPARM+LRDPRM) & - (LIVEM+LNODN+LNEK+LIPARM+LDNOD+LIDPRM)/RPI ) C** C**----------------------------------------------------------------- C** C** variables and arrays : C** -------------------- C** DOUBLE PRECISION T,NOD(LNOD),NOPARM(LNOPRM),RPARM(LRPARM), & RDPARM(LRDPRM),RBIG(LBIG),U(LU),RVEM(LRVEM), & EEST(LU),ERRG(LU),NRMERR(NK) INTEGER IVEM(LIVEM),NODNUM(LNODN),NEK(LNEK), & IPARM(LIPARM),DNOD(LDNOD),IDPARM(LIDPRM), & IBIG(RPI*LBIG) LOGICAL MASKL(NK,NK,MAXNG),MASKF(NK,MAXNG),LVEM(LLVEM) C*** INTEGER MYPROC,INFO,OUTFLG,GINFO,GINFO1,CLASS,NGROUP,G, & J,I,RUNIT,ERR CHARACTER*80 NAME,RFILE LOGICAL RSTRT DOUBLE PRECISION RE C*** EXTERNAL VEM630,VEM500 EXTERNAL DUMMY,USRFU,USERF,USERC,USERB C** C**----------------------------------------------------------------- C** C** The equivalence of RBIG and IBIG is very important : C** EQUIVALENCE (RBIG,IBIG) C** C**----------------------------------------------------------------- C** C** The common block PROP transports the real value RE to the C** subroutines: C** COMMON/PROP/RE RE=0.0 C** C**----------------------------------------------------------------- C** C** get task ids : C** NAME='a.out' CALL COMBGN(IVEM(200),MYPROC,LIVEM-203,IVEM(204),NAME,INFO) IF (INFO.NE.0) GOTO 9999 IVEM(201)=MYPROC IVEM(202)=0 IVEM(203)=IVEM(204) C** C**----------------------------------------------------------------- C** C** a protocol is printed only on process 1 : C** IF (MYPROC.EQ.1) THEN OUTFLG=1 ELSE OUTFLG=0 ENDIF C** C**----------------------------------------------------------------- C** C**** read the restart file : C** --------------------- C** WRITE(RFILE,9000) MYPROC RUNIT=99 9000 FORMAT('restart.p',I3.3) C** C**** If the restart file 'restart.p<NPROC>' exists, the file C** is read an the calculation is continued with veme02: C** INQUIRE(FILE=RFILE,EXIST=RSTRT) IF (RSTRT) THEN OPEN(UNIT=RUNIT,FILE=RFILE,STATUS= 'UNKNOWN', & FORM='UNFORMATTED') CALL VEM680 (RUNIT,LU,U,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, & LOUT,ERR) CLOSE(RUNIT) IF (ERR.GT.0) GOTO 9999 GOTO 5000 ENDIF C** C**----------------------------------------------------------------- C** C**** the parameters are copied into IVEM : C** ----------------------------------- C** IVEM(1)=MESH IVEM(MESH+ 2)=NK IVEM(MESH+ 3)=DIM C** C**----------------------------------------------------------------- C** C**** read the mesh from a universal file : C** ----------------------------------- C** IVEM(27)=LOUT IVEM(28)=OUTFLG IVEM(29)=9 IF (MYPROC.EQ.1) OPEN(9,FILE=INPUT,STATUS= 'UNKNOWN', & FORM='FORMATTED') CALL VEMU02 (LIVEM,IVEM,LNEK,NEK,LRPARM,RPARM,LIPARM,IPARM, & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN, & NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG) CLOSE(9) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** print mesh on processor 1 C** ------------------------- C** IVEM(20)=LOUT IVEM(21)=0000*OUTFLG IVEM(22)=2 CALL VEMU01(LIVEM,IVEM,LNEK,NEK,LRPARM,RPARM,LIPARM,IPARM, & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM, & LNODN,NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** distribute mesh : C** ---------------- C** IVEM(80)=LOUT IVEM(81)=OUTFLG IVEM(51)=2 CALL VEMDIS (LIVEM,IVEM,LNEK,NEK,LRPARM,RPARM,LIPARM,IPARM , & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM, & LNODN,NODNUM,LNOD,NOD,LNOPRM,NOPARM, & LBIG,RBIG,IBIG) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** set masks : C** --------- C** NGROUP=IVEM(IVEM(1)+4) GINFO=IVEM(IVEM(1)+21)+IVEM(1) GINFO1=IVEM(IVEM(1)+22) DO 400 G=1,NGROUP DO 410 I=1,NK MASKF(I,G)=.FALSE. DO 410 J=1,NK 410 MASKL(I,J,G)=.FALSE. CLASS=IVEM(GINFO+GINFO1*(G-1)+3) IF (CLASS.EQ.DIM) THEN DO 420 I=1,NK MASKF(I,G)=.TRUE. DO 420 J=1,NK 420 MASKL(I,J,G)=.TRUE. ELSEIF (CLASS.EQ.DIM-1) THEN DO 430 I=1,NK 430 MASKF(I,G)=.TRUE. ENDIF 400 CONTINUE C** C**----------------------------------------------------------------- C** C**** call of VECFEM : C** -------------- C** 5000 CONTINUE OPEN(10,FORM='UNFORMATTED',STATUS='SCRATCH') OPEN(11,FORM='UNFORMATTED',STATUS='SCRATCH') OPEN(12,FORM='UNFORMATTED',STATUS='SCRATCH') LVEM(1)=.FALSE. LVEM(4)=.FALSE. LVEM(5)=.FALSE. LVEM(6)=.TRUE. LVEM(7)=.TRUE. LVEM(8)=.TRUE. LVEM(9)=.FALSE. LVEM(10)=.TRUE. LVEM(11)=.FALSE. RVEM(1)=10. RVEM(3)=1.D-2 RVEM(10)=1.D-8 IVEM(3)=0 IVEM(10)=10 IVEM(11)=11 IVEM(12)=12 IVEM(40)=LOUT IVEM(41)=50*OUTFLG IVEM(45)=500 IVEM(46)=0 IVEM(60)=0 IVEM(70)=10 IVEM(71)=11 IVEM(72)=10 000 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,VEM500,VEM630) IF (IVEM(2).GT.1) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** save the mash files for restart : C** ------------------------------- C** IF (IVEM(2).EQ.1) THEN OPEN(UNIT=RUNIT,FILE=RFILE,FORM='UNFORMATTED') CALL VEM681 (RUNIT,LU,U,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, & LOUT,ERR) CLOSE(RUNIT) GOTO 9999 ENDIF C** C**----------------------------------------------------------------- C** C**** compute the error on the geometrical nodes : C** ------------------------------------------ C** IVEM(4)=30 IVEM(30)=LOUT IVEM(31)=OUTFLG*0 IVEM(32)=LNODN IVEM(33)=NK CALL VEMU05 (T,LU,ERRG,LU,U,LIVEM,IVEM, & LNEK, NEK ,LRPARM ,RPARM ,LIPARM ,IPARM , & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN, & NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG, & USERC) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** print the error and its norm : C** ---------------------------- C** IVEM(23)=LOUT IVEM(24)=OUTFLG IVEM(25)=IVEM(32) IVEM(26)=IVEM(33) CALL VEMU13 (LU,ERRG,NRMERR,LIVEM,IVEM, & LNEK, NEK ,LRPARM ,RPARM ,LIPARM ,IPARM , & LDNOD,DNOD,LRDPRM,RDPARM,LIDPRM,IDPARM,LNODN, & NODNUM,LNOD,NOD,LNOPRM,NOPARM,LBIG,RBIG,IBIG) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** 9999 CALL COMEND(IVEM(200),INFO) E N D SUBROUTINE USERB(T,COMPU,RHS, & NRSDP,RSDPRM,NRVDP,RVDP1,RVDPRM, & NISDP,ISDPRM,NIVDP,IVDP1,IVDPRM, & NDC,DIM,X,NOP,NOPARM,B) C** C******************************************************************* C** C** the routine USERB sets the Dirichlet boundary conditions, C** see userb. here the solution is set to zero, which is C** the exact solution at the origin. C** C******************************************************************* C** INTEGER COMPU,RHS,NRSDP,NRVDP,RVDP1,NISDP,NIVDP,IVDP1, & NDC,DIM,NOP DOUBLE PRECISION T,RSDPRM(NRSDP),RVDPRM(RVDP1,NRVDP), & X(NDC,DIM),NOPARM(NDC,NOP),B(NDC) INTEGER ISDPRM(NISDP),IVDPRM(IVDP1,NIVDP) C** C**----------------------------------------------------------------- C** INTEGER Z C** C**----------------------------------------------------------------- C** C**** start of calculation : C** -------------------- C** DO 10 Z=1,NDC B(Z) = 0 10 CONTINUE C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USERB-------------------------------------------------- E N D SUBROUTINE USRFU(T,GROUP,CLASS,COMPV,COMPU,LAST, & NELIS,L,DIM,X,TAU,NK,U,DUDX, & LT,UT,DUTDX,NOP,NOPARM,DNOPDX, & NRSP,RSPARM,NRVP,RVP1,RVPARM, & NISP,ISPARM,NIVP,IVP1,IVPARM, & F1UX,F1U,F0UX,F0U) C** C******************************************************************* C** C** the routine USRFU sets the Frechet derivative of the linear C** form F, see usrfu: C** C******************************************************************* C** INTEGER GROUP,CLASS,COMPV,COMPU,LAST,NELIS,L,LT, & DIM,NK,NOP,NRSP,RVP1,NRVP,NISP,IVP1,NIVP DOUBLE PRECISION T,X(L,DIM),TAU(L,DIM,CLASS),U(L,NK),UT(LT,NK), & DUDX(L,NK,CLASS),DUTDX(LT,NK,CLASS), & NOPARM(L,NOP),DNOPDX(L,NOP,CLASS), & RSPARM(NRSP),RVPARM(RVP1,NRVP), & F1UX(L,CLASS,CLASS),F1U(L,CLASS),F0UX(L,CLASS), & F0U(L) INTEGER ISPARM(NISP),IVPARM(IVP1,NIVP) C** C**----------------------------------------------------------------- C** INTEGER Z DOUBLE PRECISION RE COMMON/PROP/RE C** C**----------------------------------------------------------------- C** C**** start of calculation : C** --------------------- C** IF ((COMPV.EQ.1).AND.(COMPU.EQ.1).AND.(CLASS.EQ.3)) THEN DO 4 Z=1,NELIS F0U(Z)=RE*DUDX(Z,1,1) F0UX(Z,1)=RE*U(Z,1) F0UX(Z,2)=RE*U(Z,2) F0UX(Z,3)=RE*U(Z,3) F1UX(Z,1,1)=1 F1UX(Z,2,2)=1 F1UX(Z,3,3)=1 4 CONTINUE ENDIF IF ((COMPV.EQ.1).AND.(COMPU.EQ.2).AND.(CLASS.EQ.3)) THEN DO 8 Z=1,NELIS F0U(Z)=RE*DUDX(Z,1,2) 8 CONTINUE ENDIF IF ((COMPV.EQ.1).AND.(COMPU.EQ.3).AND.(CLASS.EQ.3)) THEN DO 12 Z=1,NELIS F0U(Z)=RE*DUDX(Z,1,3) 12 CONTINUE ENDIF IF ((COMPV.EQ.2).AND.(COMPU.EQ.1).AND.(CLASS.EQ.3)) THEN DO 16 Z=1,NELIS F0U(Z)=RE*DUDX(Z,2,1) 16 CONTINUE ENDIF IF ((COMPV.EQ.2).AND.(COMPU.EQ.2).AND.(CLASS.EQ.3)) THEN DO 20 Z=1,NELIS F0U(Z)=RE*DUDX(Z,2,2) F0UX(Z,1)=RE*U(Z,1) F0UX(Z,2)=RE*U(Z,2) F0UX(Z,3)=RE*U(Z,3) F1UX(Z,1,1)=1 F1UX(Z,2,2)=1 F1UX(Z,3,3)=1 20 CONTINUE ENDIF IF ((COMPV.EQ.2).AND.(COMPU.EQ.3).AND.(CLASS.EQ.3)) THEN DO 24 Z=1,NELIS F0U(Z)=RE*DUDX(Z,2,3) 24 CONTINUE ENDIF IF ((COMPV.EQ.3).AND.(COMPU.EQ.1).AND.(CLASS.EQ.3)) THEN DO 28 Z=1,NELIS F0U(Z)=RE*DUDX(Z,3,1) 28 CONTINUE ENDIF IF ((COMPV.EQ.3).AND.(COMPU.EQ.2).AND.(CLASS.EQ.3)) THEN DO 32 Z=1,NELIS F0U(Z)=RE*DUDX(Z,3,2) 32 CONTINUE ENDIF IF ((COMPV.EQ.3).AND.(COMPU.EQ.3).AND.(CLASS.EQ.3)) THEN DO 36 Z=1,NELIS F0U(Z)=RE*DUDX(Z,3,3) F0UX(Z,1)=RE*U(Z,1) F0UX(Z,2)=RE*U(Z,2) F0UX(Z,3)=RE*U(Z,3) F1UX(Z,1,1)=1 F1UX(Z,2,2)=1 F1UX(Z,3,3)=1 36 CONTINUE ENDIF C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USRFU-------------------------------------------------- E N D SUBROUTINE USERF (T,GROUP,CLASS,COMPV,RHS,LAST, & NELIS,L,DIM,X,TAU,NK,U,DUDX, & LT,UT,DUTDX,NOP,NOPARM,DNOPDX, & NRSP,RSPARM,NRVP,RVP1,RVPARM, & NISP,ISPARM,NIVP,IVP1,IVPARM, & F1,F0) C** C******************************************************************* C** C** the routine USERF sets the coefficients of the linear form F, C** see userf. C** C******************************************************************* C** INTEGER GROUP,CLASS,COMPV,RHS,LAST,NELIS,L,LT,DIM,NK,NOP, & NRSP,RVP1,NRVP,NISP,IVP1,NIVP DOUBLE PRECISION T,X(L,DIM),TAU(L,DIM,CLASS),U(L,NK),UT(LT,NK), & DUDX(L,NK,CLASS),DUTDX(LT,NK,CLASS), & NOPARM(L,NOP),DNOPDX(L,NOP,CLASS), & RSPARM(NRSP),RVPARM(RVP1,NRVP), & F1(L,CLASS),F0(L) INTEGER ISPARM(NISP),IVPARM(IVP1,NIVP) C** C**----------------------------------------------------------------- C** INTEGER Z DOUBLE PRECISION RE COMMON/PROP/RE C** C**----------------------------------------------------------------- C** C**** start of calculation : C** -------------------- C** C** the coefficients for the volume integration : C** IF (CLASS.EQ.3) THEN IF (COMPV.EQ.1) THEN DO 4 Z=1,NELIS F1(Z,1)=DUDX(Z,1,1) F1(Z,2)=DUDX(Z,1,2) F1(Z,3)=DUDX(Z,1,3) F0(Z)=RE*(U(Z,1)*DUDX(Z,1,1)+ & U(Z,2)*DUDX(Z,1,2)+U(Z,3)*DUDX(Z,1,3)) & +RE*COS(X(Z,1))*SIN(X(Z,1)) & -2*RE*COS(X(Z,1))*SIN(X(Z,1))*COS(X(Z,3))**2 & -2*RE*COS(X(Z,1))*SIN(X(Z,1))*COS(X(Z,2))**2 & +3*RE*COS(X(Z,1))*SIN(X(Z,1))*COS(X(Z,2))**2*COS(X(Z,3))**2 & -3*COS(X(Z,1))*SIN(X(Z,2))*SIN(X(Z,3)) 4 CONTINUE ENDIF IF (COMPV.EQ.2) THEN DO 8 Z=1,NELIS F1(Z,1)=DUDX(Z,2,1) F1(Z,2)=DUDX(Z,2,2) F1(Z,3)=DUDX(Z,2,3) F0(Z)=RE*(U(Z,1)*DUDX(Z,2,1)+ & U(Z,2)*DUDX(Z,2,2)+U(Z,3)*DUDX(Z,2,3)) & +RE*COS(X(Z,2))*SIN(X(Z,2)) & -2*RE*COS(X(Z,1))**2*SIN(X(Z,2))*COS(X(Z,2)) & -2*RE*COS(X(Z,2))*SIN(X(Z,2))*COS(X(Z,3))**2 & +3*RE*COS(X(Z,1))**2*SIN(X(Z,2))*COS(X(Z,2))*COS(X(Z,3))**2 & -3*SIN(X(Z,1))*COS(X(Z,2))*SIN(X(Z,3)) 8 CONTINUE ENDIF IF (COMPV.EQ.3) THEN DO 12 Z=1,NELIS F1(Z,1)=DUDX(Z,3,1) F1(Z,2)=DUDX(Z,3,2) F1(Z,3)=DUDX(Z,3,3) F0(Z)=RE*(U(Z,1)*DUDX(Z,3,1)+ & U(Z,2)*DUDX(Z,3,2)+U(Z,3)*DUDX(Z,3,3)) & +RE*COS(X(Z,3))*SIN(X(Z,3)) & -2*RE*COS(X(Z,1))**2*SIN(X(Z,3))*COS(X(Z,3)) & -2*RE*COS(X(Z,2))**2*SIN(X(Z,3))*COS(X(Z,3)) & +3*RE*COS(X(Z,1))**2*SIN(X(Z,3))*COS(X(Z,3))*COS(X(Z,2))**2 & -3*SIN(X(Z,1))*SIN(X(Z,2))*COS(X(Z,3)) 12 CONTINUE ENDIF ENDIF C** C**----------------------------------------------------------------- C** C** the coefficients for the area integration : C** IF (CLASS.EQ.2) THEN IF (COMPV.EQ.1) THEN DO 3 Z=1,NELIS IF (X(Z,2).LT.1.D-2) F0(Z)=COS(X(Z,1))*SIN(X(Z,3)) IF (X(Z,2).GT.3.14) F0(Z)=COS(X(Z,1))*SIN(X(Z,3)) IF (X(Z,3).LT.1.D-2) F0(Z)=COS(X(Z,1))*SIN(X(Z,2)) IF (X(Z,3).GT.3.14) F0(Z)=COS(X(Z,1))*SIN(X(Z,2)) 3 CONTINUE ENDIF IF (COMPV.EQ.2) THEN DO 7 Z=1,NELIS IF (X(Z,1).LT.1.D-2) F0(Z)=COS(X(Z,2))*SIN(X(Z,3)) IF (X(Z,1).GT.3.14) F0(Z)=COS(X(Z,2))*SIN(X(Z,3)) IF (X(Z,3).LT.1.D-2) F0(Z)=SIN(X(Z,1))*COS(X(Z,2)) IF (X(Z,3).GT.3.14) F0(Z)=SIN(X(Z,1))*COS(X(Z,2)) 7 CONTINUE ENDIF IF (COMPV.EQ.3) THEN DO 11 Z=1,NELIS IF (X(Z,1).LT.1.D-2) F0(Z)=SIN(X(Z,2))*COS(X(Z,3)) IF (X(Z,1).GT.3.14) F0(Z)=SIN(X(Z,2))*COS(X(Z,3)) IF (X(Z,2).LT.1.D-2) F0(Z)=SIN(X(Z,1))*COS(X(Z,3)) IF (X(Z,2).GT.3.14) F0(Z)=SIN(X(Z,1))*COS(X(Z,3)) 11 CONTINUE ENDIF ENDIF C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USERF------------------------------------------------- E N D SUBROUTINE USERC(T,GROUP,LAST,NELIS, & NRSP,RSPARM,NRVP,RVP1,RVPARM, & NISP,ISPARM,NIVP,IVP1,IVPARM, & L,DIM,X,NK,U,DUDX,NOP,NOPARM,DNOPDX,N,CU) C** C******************************************************************* C** C** the routine USERC computes in this case the error of the C** computed solution, see userc. C** C******************************************************************* C** INTEGER GROUP,LAST,NELIS,L,DIM,NK,N, & NRSP,RVP1,NRVP,NISP,IVP1,NIVP,NOP DOUBLE PRECISION T,X(L,DIM),U(L,NK),DUDX(L,NK,DIM), & RSPARM(NRSP),RVPARM(RVP1,NRVP), & NOPARM(L,NOP),DNOPDX(L,NOP,DIM),CU(L,N) INTEGER ISPARM(NISP),IVPARM(IVP1,NIVP) C** C**----------------------------------------------------------------- C** INTEGER Z C** C**----------------------------------------------------------------- C** C**** start of calculation : C** -------------------- C** DO 10 Z=1,NELIS CU(Z,1) = ABS( U(Z,1) - COS(X(Z,1))*SIN(X(Z,2))*SIN(X(Z,3))) CU(Z,2) = ABS( U(Z,2) - SIN(X(Z,1))*COS(X(Z,2))*SIN(X(Z,3))) CU(Z,3) = ABS( U(Z,3) - SIN(X(Z,1))*SIN(X(Z,2))*COS(X(Z,3))) 10 CONTINUE C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USERC-------------------------------------------------- E N D