C******************************************************************* C** C** v e m p 0 2 e x m 0 3 : C** C** velocity driven, time-dependent diffusion on the 3-dimensional C** unit cube. For the solution the nonsteady solver is used. C** The isoparametrical mesh is generated distributed on the C** processors. C** C** by L. Grosz Karlsruhe, Jan. 1995 C** C** C******************************************************************* C** C** B4 C** (0,1,2) / (1,1,2) C** *----------* C** / : B6 / | C** / : / | C** B5 *----------* B3 | C** \ | *----|-----* (1,1,0) C** | / | / C** | / B2 | / C** *----------* C** x3 ^ x2 \ (1,0,0) C** | / B1 C** -->x1 C** C** The problem is the velocity driven diffusion problem, which C** is solved by the nonsteady solver vemp02. The domain is the C** 3-dimensional [0,1] x [0,1] x [0,2] cube. An all C** boundaries Neuman boundary conditions are prescribed and on C** one point a Dirichlet condition is set. C** Using the notations in equation the problem is given by C** the functional equation: C** C** Dirichlet conditions: 0<t<=4 C** u1=b C** C** Initial solution: at t=0 C** u1=u01 C** C** functional equation: F{t,ut,u}(v)=0 on 0<t<=4 with C** C** F{t,ut,u}(v):= volume{v1x1 * u1x1 + v1x2 * u1x2 + v1x3 * u1x3 C** + v1 * ( ut1 + w1 * u1x1 + w2 * u1x2 + w3 * u1x3 + f)} C** + area{v1 * g} C** C** The functions b, u01, f and g are selected so that C** u1=x3*cos(t^2) is the exact solution of this problem. We set C** w1=w2=0 and w3=16*x1*x2*(1-x1)*(1-x2). C** C** The domain is subdivided into hexahedron elements of C** order 2. Therefore the boundary is subdivided into C** quadrilateral elements of order 2. The mesh is generated C** distributed onto the processes. The error of the computed C** solution approximation is calculated. C** PROGRAM VEMEXM C** C**----------------------------------------------------------------- C** IMPLICIT NONE include 'bytes.h' C** C**----------------------------------------------------------------- C** C** some parameters which may be chanced: C** C** NPROC = number of processors C** ELEM1 = number of elements in x1 direction, C** in x2 direction also ELEM1 elements will be C** generated, but only about ELEM1/NPROC on this C** process. C** STORE = total storage of process in Mbytes. C** INTEGER NPROC,ELEM1,STORE PARAMETER (NPROC=1, & ELEM1=5, & STORE=25) C** C**----------------------------------------------------------------- C** C** ELEM2 = number of elements in x2 direction on process C** ELEM3 = number of elements in x3 direction on process C** =2*ELEM2, since the channel has the length 2. C** N1,N2,N3 = number of nodes in x1,x2,x3-direction on the C** process C** C** other parameters are explained in mesh. C** INTEGER NK,NGROUP,DIM,MESH,GINFO,GINFO1,DINFO,DINFO1, & N1,N2,N3,ELEM2,ELEM3,LOUT PARAMETER (ELEM2=ELEM1,ELEM3=(2*ELEM1+NPROC-1)/NPROC, & N1=2*ELEM1+1, & N2=2*ELEM2+1, & N3=2*ELEM3+1, & NK=1, & NGROUP=2, & DIM=3, & MESH =210+NPROC, & GINFO =30, & GINFO1=23+2*NK, & DINFO =GINFO+GINFO1*NGROUP, & DINFO1=17, & LOUT=6) C** C**----------------------------------------------------------------- C** C** the length of the array for the mesh are set: C** they are a little bit greater than actual used in the C** mesh generation. this is necessary for the mesh distribution. C** INTEGER NN,LU,LNODN,LNOD,LNOPRM,LNEK,LRPARM,LIPARM, & LDNOD,LIDPRM,LRDPRM,LIVEM,LRVEM,LLVEM,LBIG PARAMETER (NN=N1*N2*N3*1.5, & LU =NN*NK, & LNODN =NN, & LNOD =NN*DIM, & LNOPRM=1, & LNEK=(40*(ELEM1*ELEM2*ELEM3)+ & 32*(ELEM1*ELEM2+ELEM3*ELEM1+ELEM3*ELEM2))*1.5, & LIPARM=(ELEM1*ELEM2*ELEM3+ & 2*(ELEM1*ELEM2+ELEM3*ELEM1+ELEM3*ELEM2))*1.5, & LRPARM=1, & LDNOD =2*NK*1.5, & LIDPRM=NK*1.5, & LRDPRM=1, & LIVEM =MESH+DINFO+DINFO1*NK+600+LU+LDNOD/2, & LLVEM =500, & LRVEM =60+15*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,NGROUP),MASKF(NK,NGROUP),LVEM(LLVEM) C*** INTEGER MYPROC,INFO,OUTFLG,NDNUM0,HERE,S,NE1,ADGEO1, & NE2,ADGEO2,ADIVP2,ADIVP1,NE0,NDC1,NDC2, & ADDCG1,ADDCG2,ELNUM0,ELMID,SPACE,LSPACE,STEP INTEGER Z1,Z2,Z3 DOUBLE PRECISION X30 CHARACTER*80 NAME C*** EXTERNAL VEM630,VEM500 EXTERNAL DUMMY,USRFU,USERF,USERC,USERB,USRFUT,USERU0 C** C**----------------------------------------------------------------- C** C** The equivalence of RBIG and IBIG is very important : C** EQUIVALENCE (RBIG,IBIG) C** C**----------------------------------------------------------------- C** C** get task ids : C** NAME='a.out' IVEM(200)=NPROC 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) IF (NPROC.NE.IVEM(200)) THEN PRINT*,'Set NPROC=',IVEM(200) GOTO 9999 ENDIF 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**** the parameters are copied into IVEM : C** ----------------------------------- C** IVEM(1)=MESH IVEM(MESH+ 1)=N1*N2*N3 IVEM(MESH+ 2)=NK IVEM(MESH+ 3)=DIM IVEM(MESH+ 4)=NGROUP IVEM(MESH+ 5)=NN IVEM(MESH+13)=NN IVEM(MESH+14)=0 IVEM(MESH+15)=0 IVEM(MESH+18)=0 IVEM(MESH+21)=GINFO IVEM(MESH+22)=GINFO1 IVEM(MESH+23)=DINFO IVEM(MESH+24)=DINFO1 C** C**----------------------------------------------------------------- C** C** This process generates the nodes in the subdomain C** [0,1] x [0,1] x [ X30,X30+2/NPROC] starting with the node id C** number NDNUM0. The nodes with x3=X30 are also generated on C** process MYPROC-1 and the nodes with x3=X30+2/NPROC are also C** generated on process MYPROC+1. The first element generated C** on the process gets the element id number ELNUM0. C** X30=2*DBLE(MYPROC-1)/DBLE(IVEM(200)) NDNUM0=(MYPROC-1)*N1*N2*(N3-1)+1 ELNUM0=(MYPROC-1)*(ELEM1*ELEM2*ELEM3+ & 2*(ELEM1*ELEM2+ELEM3*ELEM2+ELEM1*ELEM3))+1 C** C**----------------------------------------------------------------- C** C**** the generation of the geometrical nodes : C** --------------------------------------- C** C** the grid is regular with N1 points in x1- and N2 points in C** x2 direction. C** DO 10 Z3=1,N3 DO 10 Z2=1,N2 DO 10 Z1=1,N1 NOD(Z1+N1*(Z2-1)+N1*N2*(Z3-1) )=DBLE(Z1-1)/DBLE(N1-1) NOD(Z1+N1*(Z2-1)+N1*N2*(Z3-1)+ NN)=DBLE(Z2-1)/DBLE(N2-1) NOD(Z1+N1*(Z2-1)+N1*N2*(Z3-1)+2*NN)= & 2*DBLE(Z3-1)/DBLE(IVEM(200)*(N3-1))+X30 NODNUM(Z1+N1*(Z2-1)+N1*N2*(Z3-1))=Z1+N1*(Z2-1)+N1*N2*(Z3-1) & +NDNUM0-1 10 CONTINUE C** C**----------------------------------------------------------------- C** C**** the generation of the elements : C** ------------------------------- C** C** The domain is covered by hexahedron elements of order 2 C** and consequently the boundaries are described by C** quadrilateral elements of order 2. The succession of the C** nodes in the element is defined in vemu02 and vembuild(3). C** The lowest node id in an element is S. C** C** ADGEO1 defines the start address of the hexahedrons C** elements in NEK and ADIVP1 defines the start address of C** the element id number assigned to the elements. The element C** id number is unique over all processes. NE1 is the number of C** hexahedrons elements generated on the process. HERE gives C** the address of the element in NEK, which lowest vertex has C** the node id S over all processes. C** ADGEO1=1 ADIVP1=1 NE1=ELEM1*ELEM2*ELEM3 DO 20 Z3=1,ELEM3 DO 20 Z2=1,ELEM2 DO 20 Z1=1,ELEM1 S=2*(Z1-1)+2*(Z2-1)*N1+2*N1*N2*(Z3-1)+NDNUM0 HERE=Z1+ELEM1*(Z2-1)+ELEM1*ELEM2*(Z3-1)+ADGEO1-1 ELMID=Z1+ELEM1*(Z2-1)+ELEM1*ELEM2*(Z3-1)+ELNUM0 IPARM(ADIVP1-1+Z1+ELEM1*(Z2-1)+ELEM1*ELEM2*(Z3-1))=ELMID NEK(HERE )=S NEK(HERE+ NE1)=S+2 NEK(HERE+ 2*NE1)=S+2*N1+2 NEK(HERE+ 3*NE1)=S+2*N1 NEK(HERE+ 4*NE1)=S+2*N1*N2 NEK(HERE+ 5*NE1)=S+2*N1*N2+2 NEK(HERE+ 6*NE1)=S+2*N1*N2+2*N1+2 NEK(HERE+ 7*NE1)=S+2*N1*N2+2*N1 NEK(HERE+ 8*NE1)=S+1 NEK(HERE+ 9*NE1)=S+N1+2 NEK(HERE+10*NE1)=S+2*N1+1 NEK(HERE+11*NE1)=S+N1 NEK(HERE+12*NE1)=S+N1*N2 NEK(HERE+13*NE1)=S+N1*N2+2 NEK(HERE+14*NE1)=S+N1*N2+2*N1+2 NEK(HERE+15*NE1)=S+N1*N2+2*N1 NEK(HERE+16*NE1)=S+2*N1*N2+1 NEK(HERE+17*NE1)=S+2*N1*N2+N1+2 NEK(HERE+18*NE1)=S+2*N1*N2+2*N1+1 NEK(HERE+19*NE1)=S+2*N1*N2+N1 20 CONTINUE C** C** ADGEO2 defines the start address of the line elements C** in NEK and ADIVP2 defines the start address of the C** element id number assigned to the elements. The entries 1 to C** 20*NE1 in NEK and 1 to NE1 in IPARM are already used by C** the elements in group 1. NE2 is the number of C** quadrilateral elements generated on the process, where the C** elements on boundary B1/B6 are only generated on process 1 C** or NPROC. HERE gives the address of the element in NEK, which C** is a boundary element of the hexahedrons element with lowest C** node id S. C** ADGEO2=ADGEO1+20*NE1 ADIVP2=ADIVP1+NE1 NE2=2*(ELEM1+ELEM2)*ELEM3 IF (MYPROC.EQ.1) NE2=NE2+ELEM1*ELEM2 IF (MYPROC.EQ.IVEM(200)) NE2=NE2+ELEM1*ELEM2 NE0=0 C** C** these are the quadrilateral elements on boundary 1 (x3=0): C** only on process 1. NE0 counts the already generated line C** elements C** C**** elements on boundary 1 (x3=0): (only on process 1) C** IF (MYPROC.EQ.1) THEN DO 31 Z2=1,ELEM2 DO 31 Z1=1,ELEM1 HERE=Z1+ELEM1*(Z2-1)+NE0+ADGEO2-1 ELMID=Z1+ELEM1*(Z2-1)+ELEM1*ELEM2*ELEM3+NE0+ELNUM0 S=2*(Z1-1)+2*(Z2-1)*N1+NDNUM0 IPARM(ADIVP2-1+NE0+Z1+ELEM1*(Z2-1))=ELMID NEK(HERE )= S NEK(HERE+ NE2)= S+2*N1 NEK(HERE+2*NE2)= S+2*N1+2 NEK(HERE+3*NE2)= S+2 NEK(HERE+4*NE2)= S+N1 NEK(HERE+5*NE2)= S+2*N1+1 NEK(HERE+6*NE2)= S+N1+2 NEK(HERE+7*NE2)= S+1 31 CONTINUE NE0=NE0+ELEM1*ELEM2 ENDIF C** C**** elements on boundary 6 (x3=2): (only on process NPROC) C** IF (MYPROC.EQ.IVEM(200)) THEN DO 32 Z2=1,ELEM2 DO 32 Z1=1,ELEM1 HERE=Z1+ELEM1*(Z2-1)+NE0+ADGEO2-1 ELMID=Z1+ELEM1*(Z2-1)+ELEM1*ELEM2*ELEM3+NE0+ELNUM0 S=2*(Z1-1)+2*(Z2-1)*N1+2*N1*N2*ELEM3+NDNUM0 IPARM(ADIVP2-1+NE0+Z1+ELEM1*(Z2-1))=ELMID NEK(HERE )= S NEK(HERE+ NE2)= S+2 NEK(HERE+2*NE2)= S+2*N1+2 NEK(HERE+3*NE2)= S+2*N1 NEK(HERE+4*NE2)= S+1 NEK(HERE+5*NE2)= S+N1+2 NEK(HERE+6*NE2)= S+2*N1+1 NEK(HERE+7*NE2)= S+N1 32 CONTINUE NE0=NE0+ELEM1*ELEM2 ENDIF C** C**** elements on boundary 5 (x1=0): C** DO 33 Z3=1,ELEM3 DO 33 Z2=1,ELEM2 HERE=Z2+ELEM2*(Z3-1)+NE0+ADGEO2-1 ELMID=Z2+ELEM2*(Z3-1)+ELEM1*ELEM2*ELEM3+NE0+ELNUM0 S=2*(Z2-1)*N1+2*N1*N2*(Z3-1)+NDNUM0 IPARM(ADIVP2-1+NE0+Z2+ELEM2*(Z3-1))=ELMID NEK(HERE )= S NEK(HERE+ NE2)= S+2*N1*N2 NEK(HERE+2*NE2)= S+2*N1*N2+2*N1 NEK(HERE+3*NE2)= S+2*N1 NEK(HERE+4*NE2)= S+N1*N2 NEK(HERE+5*NE2)= S+2*N1*N2+N1 NEK(HERE+6*NE2)= S+N1*N2+2*N1 NEK(HERE+7*NE2)= S+N1 33 CONTINUE NE0=NE0+ELEM3*ELEM2 C** C**** elements on boundary 3 (x1=1): C** DO 34 Z3=1,ELEM3 DO 34 Z2=1,ELEM2 HERE=Z2+ELEM2*(Z3-1)+NE0+ADGEO2-1 ELMID=Z2+ELEM2*(Z3-1)+ELEM1*ELEM2*ELEM3+NE0+ELNUM0 S=2*ELEM1+2*(Z2-1)*N1+2*N1*N2*(Z3-1)+NDNUM0 IPARM(ADIVP2-1+NE0+Z2+ELEM2*(Z3-1))=ELMID NEK(HERE )= S NEK(HERE+ NE2)= S+2*N1 NEK(HERE+2*NE2)= S+2*N1*N2+2*N1 NEK(HERE+3*NE2)= S+2*N1*N2 NEK(HERE+4*NE2)= S+N1 NEK(HERE+5*NE2)= S+N1*N2+2*N1 NEK(HERE+6*NE2)= S+2*N1*N2+N1 NEK(HERE+7*NE2)= S+N1*N2 34 CONTINUE NE0=NE0+ELEM3*ELEM2 C** C**** elements on boundary 2 (x2=0): C** DO 35 Z3=1,ELEM3 DO 35 Z1=1,ELEM1 HERE=Z1+ELEM1*(Z3-1)+NE0+ADGEO2-1 ELMID=Z3+ELEM3*(Z1-1)+ELEM1*ELEM2*ELEM3+NE0+ELNUM0 S=2*(Z1-1)+2*N1*N2*(Z3-1)+NDNUM0 IPARM(ADIVP2-1+NE0+Z1+ELEM1*(Z3-1))=ELMID NEK(HERE )= S NEK(HERE+ NE2)= S+2 NEK(HERE+2*NE2)= S+2*N2*N1+2 NEK(HERE+3*NE2)= S+2*N1*N2 NEK(HERE+4*NE2)= S+1 NEK(HERE+5*NE2)= S+N1*N2+2 NEK(HERE+6*NE2)= S+2*N2*N1+1 NEK(HERE+7*NE2)= S+N1*N2 35 CONTINUE NE0=NE0+ELEM3*ELEM1 C** C**** elements on boundary 4 (x2=1): C** DO 36 Z3=1,ELEM3 DO 36 Z1=1,ELEM1 HERE=Z1+ELEM1*(Z3-1)+NE0+ADGEO2-1 ELMID=Z3+ELEM3*(Z1-1)+ELEM1*ELEM2*ELEM3+NE0+ELNUM0 S=2*(Z1-1)+2*ELEM2*N1+2*N1*N2*(Z3-1)+NDNUM0 IPARM(ADIVP2-1+NE0+Z1+ELEM1*(Z3-1))=ELMID NEK(HERE )= S NEK(HERE+ NE2)= S+2*N1*N2 NEK(HERE+2*NE2)= S+2*N1*N2+2 NEK(HERE+3*NE2)= S+2 NEK(HERE+4*NE2)= S+N1*N2 NEK(HERE+5*NE2)= S+2*N2*N1+1 NEK(HERE+6*NE2)= S+N1*N2+2 NEK(HERE+7*NE2)= S+1 36 CONTINUE C** C** C**----------------------------------------------------------------- C** C** the start addresses, etc are written to IVEM: C** C** group 1: hexahedrons elements C** IVEM(MESH+GINFO ) = NE1 IVEM(MESH+GINFO+ 2) = 8 IVEM(MESH+GINFO+ 3) = 3 IVEM(MESH+GINFO+ 8) = 0 IVEM(MESH+GINFO+11) = 0 IVEM(MESH+GINFO+13) = 0 IVEM(MESH+GINFO+14) = ADIVP1 IVEM(MESH+GINFO+15) = NE1 IVEM(MESH+GINFO+16) = 1 IVEM(MESH+GINFO+20) = ADGEO1 IVEM(MESH+GINFO+21) = NE1 IVEM(MESH+GINFO+23) = 20 C** C** group 2: quadrilateral elements C** IVEM(MESH+GINFO+GINFO1 ) = NE2 IVEM(MESH+GINFO+GINFO1+ 2) = 4 IVEM(MESH+GINFO+GINFO1+ 3) = 2 IVEM(MESH+GINFO+GINFO1+ 8) = 0 IVEM(MESH+GINFO+GINFO1+11) = 0 IVEM(MESH+GINFO+GINFO1+13) = 0 IVEM(MESH+GINFO+GINFO1+14) = ADIVP2 IVEM(MESH+GINFO+GINFO1+15) = NE2 IVEM(MESH+GINFO+GINFO1+16) = 1 IVEM(MESH+GINFO+GINFO1+20) = ADGEO2 IVEM(MESH+GINFO+GINFO1+21) = NE2 IVEM(MESH+GINFO+GINFO1+23) = 8 C** C**----------------------------------------------------------------- C** C** generation of the nodes with Dirichlet conditions : C** ------------------------------------------------- C** C** The node with node id number 1 gets a Dirichlet condition: C** (only on processor 1) C** IF (MYPROC.EQ.1) THEN DNOD(1)=1 NDC1=1 ELSE NDC1=0 ENDIF C** C**----------------------------------------------------------------- C** C** the start addresses, etc are written to IVEM: C** C** component 1: C** IVEM(MESH+DINFO ) = NDC1 IVEM(MESH+DINFO+ 2) = 1 IVEM(MESH+DINFO+ 4) = 0 IVEM(MESH+DINFO+ 7) = 0 IVEM(MESH+DINFO+ 9) = 0 IVEM(MESH+DINFO+12) = 0 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 the initial solution : C** ------------------------ C** IVEM(30)=LOUT IVEM(31)=OUTFLG*0 T=0. CALL VEMU08(T,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,USERU0) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C**** call of VECFEM : C** -------------- C** 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(6)=.TRUE. LVEM(7)=.TRUE. LVEM(8)=.TRUE. LVEM(9)=.FALSE. LVEM(10)=.TRUE. LVEM(11)=.FALSE. LVEM(15)=.FALSE. LVEM(16)=.TRUE. RVEM(3)=1.D-2 RVEM(10)=1.D-3 RVEM(11)=T RVEM(12)=4. RVEM(13)=0.1 RVEM(14)=1.D-8 RVEM(15)=0 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)=10 IVEM(70)=10 IVEM(71)=1 IVEM(72)=5000 MASKL(1,1,1)=.TRUE. MASKL(1,1,2)=.FALSE. MASKF(1,1)=.TRUE. MASKF(1,2)=.TRUE. C** C** the solution is returned after T+.25 is reached : C** 9998 T=T+.25 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, & VEM500,VEM630) STEP=IVEM(3) IF (IVEM(2).NE.0) GOTO 9999 C** C** C**----------------------------------------------------------------- C** C**** compute the error on the geometrical nodes : C** ------------------------------------------ C** SPACE=IVEM(8) LSPACE=IVEM(9) IVEM(4)=30 IVEM(30)=LOUT IVEM(31)=OUTFLG*0 IVEM(32)=NN 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,LSPACE, & RBIG(SPACE),IBIG(RPI*(SPACE-1)+1),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)=1 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, & LSPACE,RBIG(SPACE),IBIG(RPI*(SPACE-1)+1)) IF (IVEM(2).NE.0) GOTO 9999 C** C**----------------------------------------------------------------- C** C** If step is equal to 1 the T-integration is continued : C** IF (STEP.EQ.1) GOTO 9998 C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** 9999 CALL COMEND(IVEM(200),INFO) E N D SUBROUTINE USERU0(T,NE,L,DIM,X,NOP,NOPARM,COMPU,U0) C** C******************************************************************* C** C** the routine USERU0 sets the initial solution, see vemu08. C** C******************************************************************* C** INTEGER NE,L,DIM,NOP,COMPU DOUBLE PRECISION T,X(L,DIM),NOPARM(L,NOP),U0(NE) C** C**----------------------------------------------------------------- C** INTEGER Z C** C**----------------------------------------------------------------- C** C**** start of calculation : C** -------------------- C** IF (COMPU.EQ.1) THEN DO 10 Z=1,NE U0(Z) = X(Z,3) * COS(T**2) 10 CONTINUE ENDIF C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USERU0------------------------------------------------- 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 exact solution x3 is prescribed. 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** IF (COMPU.EQ.1) THEN DO 10 Z=1,NDC B(Z) = X(Z,3) * COS(T**2) 10 CONTINUE ENDIF C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USERB-------------------------------------------------- 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** It is f=-16*x1*x2*(1-x1)*(1-x2) * u1x3 C** g=-(n1*u1x1+n2*u1x1+n3*u1x3)=-n3 C** C** computed by the exact solution u1=x3. (n1,n2,n3) is the outer C** normal field, which is computed by the tangential field: C** n3=(tau11*tau23-tau21*tau12)/norm , norm is the normalization. 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 NORM 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 12 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)=16*X(Z,1)*X(Z,2)*(1.-X(Z,1))*(1.-X(Z,2))* & (DUDX(Z,1,3)-COS(T**2))+UT(Z,1)+2*T*SIN(T**2)*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 11 Z=1,NELIS NORM = (TAU(Z,2,1)*TAU(Z,3,2)-TAU(Z,3,1)*TAU(Z,2,2))**2 & + (TAU(Z,1,1)*TAU(Z,3,2)-TAU(Z,3,1)*TAU(Z,1,2))**2 & + (TAU(Z,1,1)*TAU(Z,2,2)-TAU(Z,2,1)*TAU(Z,1,2))**2 F0(Z)=-(TAU(Z,1,1)*TAU(Z,2,2)-TAU(Z,2,1)*TAU(Z,1,2)) & /SQRT(NORM)*COS(T**2) 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 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 C** C**----------------------------------------------------------------- C** C**** start of calculation : C** --------------------- C** C** the coefficients for the area integration : C** IF (CLASS.EQ.3) THEN IF ((COMPV.EQ.1).AND.(COMPU.EQ.1)) THEN DO 112 Z=1,NELIS F1UX(Z,1,1)=1. F1UX(Z,2,2)=1. F1UX(Z,3,3)=1. F0UX(Z,3)=16*X(Z,1)*X(Z,2)*(1.-X(Z,1))*(1.-X(Z,2)) 112 CONTINUE ENDIF ENDIF C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USRFU-------------------------------------------------- E N D SUBROUTINE USRFUT(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, & F1UTX,F1UT,F0UTX,F0UT) C** C******************************************************************* C** C** the routine USRFU sets the Frechet derivative of F with C** respect of ut, 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), & F1UTX(L,CLASS,CLASS),F1UT(L,CLASS), & F0UTX(L,CLASS),F0UT(L) INTEGER ISPARM(NISP),IVPARM(IVP1,NIVP) C** C**----------------------------------------------------------------- C** INTEGER Z C** C**----------------------------------------------------------------- C** C**** start of calculation : C** --------------------- C** C** the coefficients for the area integration : C** IF (CLASS.EQ.3) THEN IF ((COMPV.EQ.1).AND.(COMPU.EQ.1)) THEN DO 112 Z=1,NELIS F0UT(Z)=1. 112 CONTINUE ENDIF ENDIF C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USRFUT------------------------------------------------- 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) - X(Z,3) * COS(T*T) ) 10 CONTINUE C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USERC-------------------------------------------------- E N D