|
| C**************************************************************** C** C** v e m e 0 2 e x m 0 9 : C** C** The Navier-Stokes equations on 3-dimensional domain. The mesh is C** generated distributed on the processors. C** C** by L. Grosz Canberra, Nov. 1997 C** C**************************************************************** C** C** The problem is a system of four partial differential C** equations four the three velocity components u1,u2 and u3 and C** the pressure u4. C** The domain is the 3-dimensional unit cube [0,1]^3. C** Using the notations in equation the problem is given by C** the functional equation: C** C** Dirichlet conditions: C** u1=b1 at x1=0 and x1=1 C** u2=b2 at x2=0 and x2=1 C** u3=b3 at x3=0 and x3=1 C** u4=b4 at (x1,x2,x3)=(1/2,1/2,1/2) C** C** nonlinear functional equation: F{u}(v)=0 C** C** with F{u}(v)= C** volume{v1x1*(u1x1+u4)+v1x2*u1x2+v1x3*u1x3+ C** Re*v1*(u1*u1x1+u2*u1x2+u3*u1x3+f1) C** + v2x1*u2x1+v2x2*(u2x2+u4)+v2x3*u2x3+ C** Re*v2*(u1*u2x1+u2*u2x2+u3*u2x3+f2) C** + v3x1*u3x1+v3x2*u3x2+v3x3*(u3x3+u4)+ C** Re*v3*(u1*u3x1+u2*u3x2+u3*u3x3+f3) C** + v4*(u1x1+u2x2+u3x3)} C** + area{v1*g1+v2*g2+v3*g3} C** C** The functions b, f and g are selected so that C** C** u1=x1*(x3-x2) C** u2=x2*(x1-x3) C** u3=x3*(x2-x1) C** u4=(x1+x2+x3)/3 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 of ELEM1 X ELEM2 X ELEM3 elements is generated distributed C** on the processors. The mixed FEM method with order two elements C** for the velocity components (u1,u2,u3) and order one elements for C** the pressure component u4 are used. C** C**----------------------------------------------------------------- C** PROGRAM VEMEXM C** C**----------------------------------------------------------------- C** IMPLICIT NONE include 'bytes.h' 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** ELEM2 = number of elements in x2 direction on process C** ELEM3 = number of elements in x3 direction on process C** STORE = total storage of process in Mbytes. C** INTEGER NPROC,ELEM1,ELEM2,ELEM3 REAL STORE PARAMETER (NPROC=1, & STORE=15, & ELEM1=5, & ELEM2=5, & ELEM3=5) C** C**----------------------------------------------------------------- C** C** N1,N2,N3 = number of nodes in x1,x2,x3-direction on the C** process C** E3DIM = maximal number of elements in x3-direction generated C** on a processor C** N3DIM = maximal number of nodes in x3-direction generated C** on a processor C** C** other parameters are explained in mesh. C** INTEGER NK,DIM,N1,N2,N3DIM,E3DIM,LOUT,NGROUP,N3 PARAMETER (N1=2*ELEM1+1, & N2=2*ELEM2+1, & N3=2*ELEM3+1, & E3DIM=(ELEM3/NPROC)+1, & N3DIM=2*E3DIM+1, & NK=4, & NGROUP=2, & DIM=3, & 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*N3DIM*1.5, & LU =NN*NK, & LNODN =NN, & LNOD =NN*DIM, & LNOPRM=1, & LNEK=(2*NK*20*(ELEM1*ELEM2*E3DIM)+ & 2*NK*8*2*(ELEM1*ELEM2+E3DIM*ELEM1+E3DIM*ELEM2)), & LIPARM=(ELEM1*ELEM2*E3DIM+ & 2*(ELEM1*ELEM2+E3DIM*ELEM1+E3DIM*ELEM2)), & LRPARM=1, & LDNOD =2*NK*2*(N3DIM*N2+N3DIM*N1+N2*N1), & LIDPRM=1, & LRDPRM=1, & LIVEM =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),PROPO(NK,2) LOGICAL MASKL(NK,NK,2),MASKF(NK,2),LVEM(LLVEM) C INTEGER MYPROC,INFO,OUTFLG,GINFO,GINFO1,ELEM3L,ELEM30, & N30,N3L,MESH,DINFO,DINFO1,NDNUM0,ELNUM0, & Z3,Z2,Z1,ADGEO1,ADIVP1,NE1,S,HERE,ELMID, & ADGEO2,ADIVP2,NE2,NE0,ADDC,NDC,DEF CHARACTER*80 NAME DOUBLE PRECISION RE, X30 C EXTERNAL VEM630,VEM500 EXTERNAL DUMMY,USRFU,USERF,USERC,USERB,USERU0 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=10.0 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 (IVEM(200).NE.NPROC) then IF (MYPROC.EQ.0) then PRINT*,'Number of processors does not fit !' ENDIF 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** ELEM3L is the number of elements in x3-direction C** generated on this processor. ELEM30/N30 is the first element/node C** that is generated on this processor. C** ELEM3L=ELEM3/NPROC DEF=ELEM3-ELEM3L*NPROC IF (MYPROC.LE.DEF) THEN ELEM3L=ELEM3L+1 ELEM30=ELEM3L*(MYPROC-1) ELSE ELEM30=DEF*(ELEM3L+1)+ELEM3L*(MYPROC-DEF-1) ENDIF N30=2*ELEM30+1 N3L=2*ELEM3L+1 C** C**----------------------------------------------------------------- C** C** This process generates the nodes in the subdomain C** [0,1] x [0,1] x [ X30,X30+1/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+1/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=DBLE(N30-1)/DBLE(N3-1) NDNUM0=N1*N2*(N30-1)+1 ELNUM0=ELEM1*ELEM2*ELEM30+ & 2*(ELEM1*ELEM2+ELEM30*ELEM2+ELEM1*ELEM30)+1 C** C**----------------------------------------------------------------- C** C**** the parameters are copied into IVEM : C** ----------------------------------- C** MESH =210+NPROC GINFO =30 GINFO1=23+2*NK DINFO =GINFO+GINFO1*NGROUP DINFO1=17 IVEM(1)=MESH IVEM(MESH+ 1)=N1*N2*N3L 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**** 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,N3L 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)= & DBLE(Z3-1)/DBLE(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*ELEM3L DO 20 Z3=1,ELEM3L 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)*ELEM3L IF (MYPROC.EQ.1) NE2=NE2+ELEM1*ELEM2 IF (MYPROC.EQ.NPROC) 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*ELEM3L+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.NPROC) 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*ELEM3L+NE0+ELNUM0 S=2*(Z1-1)+2*(Z2-1)*N1+2*N1*N2*ELEM3L+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,ELEM3L DO 33 Z2=1,ELEM2 HERE=Z2+ELEM2*(Z3-1)+NE0+ADGEO2-1 ELMID=Z2+ELEM2*(Z3-1)+ELEM1*ELEM2*ELEM3L+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+ELEM3L*ELEM2 C** C**** elements on boundary 3 (x1=1): C** DO 34 Z3=1,ELEM3L DO 34 Z2=1,ELEM2 HERE=Z2+ELEM2*(Z3-1)+NE0+ADGEO2-1 ELMID=Z2+ELEM2*(Z3-1)+ELEM1*ELEM2*ELEM3L+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+ELEM3L*ELEM2 C** C**** elements on boundary 2 (x2=0): C** DO 35 Z3=1,ELEM3L DO 35 Z1=1,ELEM1 HERE=Z1+ELEM1*(Z3-1)+NE0+ADGEO2-1 ELMID=Z3+ELEM3L*(Z1-1)+ELEM1*ELEM2*ELEM3L+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+ELEM3L*ELEM1 C** C**** elements on boundary 4 (x2=1): C** DO 36 Z3=1,ELEM3L DO 36 Z1=1,ELEM1 HERE=Z1+ELEM1*(Z3-1)+NE0+ADGEO2-1 ELMID=Z3+ELEM3L*(Z1-1)+ELEM1*ELEM2*ELEM3L+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** ADDC=0 C** C** component 1: C** NDC=0 DO 41 Z3=1,N3L DO 41 Z2=1,N2 NDC=NDC+1 DNOD(ADDC+NDC)=1+N1*(Z2-1)+N1*N2*(NDNUM0-1+Z3-1) NDC=NDC+1 DNOD(ADDC+NDC)=N1+N1*(Z2-1)+N1*N2*(NDNUM0-1+Z3-1) 41 CONTINUE IVEM(MESH+DINFO ) = NDC IVEM(MESH+DINFO+ 2) = ADDC+1 IVEM(MESH+DINFO+ 4) = 0 IVEM(MESH+DINFO+ 7) = 0 IVEM(MESH+DINFO+ 9) = 0 IVEM(MESH+DINFO+12) = 0 ADDC=ADDC+NDC C** C** component 2: C** NDC=0 DO 42 Z3=1,N3L DO 42 Z1=1,N1 NDC=NDC+1 DNOD(ADDC+NDC)=Z1+N1*N2*(NDNUM0-1+Z3-1) NDC=NDC+1 DNOD(ADDC+NDC)=Z1+N1*(N2-1)+N1*N2*(NDNUM0-1+Z3-1) 42 CONTINUE IVEM(MESH+DINFO+DINFO1 ) = NDC IVEM(MESH+DINFO+DINFO1+ 2) = ADDC+1 IVEM(MESH+DINFO+DINFO1+ 4) = 0 IVEM(MESH+DINFO+DINFO1+ 7) = 0 IVEM(MESH+DINFO+DINFO1+ 9) = 0 IVEM(MESH+DINFO+DINFO1+12) = 0 ADDC=ADDC+NDC C** C** component 3: C** NDC=0 IF (MYPROC.EQ.1) THEN DO 43 Z1=1,N1 DO 43 Z2=1,N2 NDC=NDC+1 DNOD(ADDC+NDC)=Z1+N1*(Z2-1) 43 CONTINUE ENDIF IF (MYPROC.EQ.NPROC) THEN DO 44 Z1=1,N1 DO 44 Z2=1,N2 NDC=NDC+1 DNOD(ADDC+NDC)=Z1+N1*(Z2-1)+N1*N2*(N3-1) 44 CONTINUE ENDIF IVEM(MESH+DINFO+2*DINFO1 ) = NDC IVEM(MESH+DINFO+2*DINFO1+ 2) = ADDC+1 IVEM(MESH+DINFO+2*DINFO1+ 4) = 0 IVEM(MESH+DINFO+2*DINFO1+ 7) = 0 IVEM(MESH+DINFO+2*DINFO1+ 9) = 0 IVEM(MESH+DINFO+2*DINFO1+12) = 0 ADDC=ADDC+NDC C** C** C** component 4: C** DNOD(ADDC+1)=N1*N2*(((ELEM3+1)/2)*2+1) IVEM(MESH+DINFO+3*DINFO1 ) = 1 IVEM(MESH+DINFO+3*DINFO1+ 2) = ADDC+1 IVEM(MESH+DINFO+3*DINFO1+ 4) = 0 IVEM(MESH+DINFO+3*DINFO1+ 7) = 0 IVEM(MESH+DINFO+3*DINFO1+ 9) = 0 IVEM(MESH+DINFO+3*DINFO1+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**** generate mixed mesh : C** ------------------- C** PROPO(1,1)=2 PROPO(2,1)=2 PROPO(3,1)=2 PROPO(4,1)=1 PROPO(1,2)=2 PROPO(2,2)=2 PROPO(3,2)=2 PROPO(4,2)=1 IVEM(101)=LOUT IVEM(102)=OUTFLG CALL VEMGE2 (PROPO,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)=3 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 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**** set masks : C** --------- C** MASKF(1,1)=.TRUE. MASKF(2,1)=.TRUE. MASKF(3,1)=.TRUE. MASKF(4,1)=.TRUE. MASKF(1,2)=.TRUE. MASKF(2,2)=.TRUE. MASKF(3,2)=.TRUE. MASKF(4,2)=.FALSE. MASKL(1,1,1)=.TRUE. MASKL(1,2,1)=.TRUE. MASKL(1,3,1)=.TRUE. MASKL(1,4,1)=.TRUE. MASKL(2,1,1)=.TRUE. MASKL(2,2,1)=.TRUE. MASKL(2,3,1)=.TRUE. MASKL(2,4,1)=.TRUE. MASKL(3,1,1)=.TRUE. MASKL(3,2,1)=.TRUE. MASKL(3,3,1)=.TRUE. MASKL(3,4,1)=.TRUE. MASKL(4,1,1)=.TRUE. MASKL(4,2,1)=.TRUE. MASKL(4,3,1)=.TRUE. MASKL(4,4,1)=.FALSE. MASKL(1,1,2)=.FALSE. MASKL(1,2,2)=.FALSE. MASKL(1,3,2)=.FALSE. MASKL(1,4,2)=.FALSE. MASKL(2,1,2)=.FALSE. MASKL(2,2,2)=.FALSE. MASKL(2,3,2)=.FALSE. MASKL(2,4,2)=.FALSE. MASKL(3,1,2)=.FALSE. MASKL(3,2,2)=.FALSE. MASKL(3,3,2)=.FALSE. MASKL(3,4,2)=.FALSE. MASKL(4,1,2)=.FALSE. MASKL(4,2,2)=.FALSE. MASKL(4,3,2)=.FALSE. MASKL(4,4,2)=.FALSE. 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(5)=.FALSE. LVEM(6)=.FALSE. LVEM(7)=.TRUE. LVEM(8)=.TRUE. LVEM(9)=.FALSE. LVEM(10)=.TRUE. LVEM(11)=.FALSE. RVEM(1)=0. RVEM(3)=1.D-2 RVEM(10)=1.D-4 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**** 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 USERU0(T,NE,L,DIM,X,NOP,NOPARM,COMPU,U0) C** C******************************************************************* C** C** USERU0 sets the initial guess for veme02 or the initial C** solution for vemp02. C** C******************************************************************* C** IMPLICIT DOUBLE PRECISION (S,T) C** C**----------------------------------------------------------------- C** C** Formal Parameters : see man vemu08 C** ----------------- C** INTEGER NE,L,DIM,NOP,COMPU DOUBLE PRECISION T,X(L,DIM),NOPARM(L,NOP),U0(NE) C** C**----------------------------------------------------------------- C** INTEGER Z Z=0 C** C**----------------------------------------------------------------- C** C**** Start of Calculation : C** -------------------- C** IF (COMPU.EQ.1) THEN DO 10 Z=1,NE U0(Z) = X(Z,1)*(X(Z,3)-X(Z,2)) 10 CONTINUE ENDIF IF (COMPU.EQ.2) THEN DO 20 Z=1,NE U0(Z) = X(Z,2)*(X(Z,1)-X(Z,3)) 20 CONTINUE ENDIF IF (COMPU.EQ.3) THEN DO 30 Z=1,NE U0(Z) = X(Z,3)*(X(Z,2)-X(Z,1)) 30 CONTINUE ENDIF IF (COMPU.EQ.4) THEN DO 40 Z=1,NE U0(Z) = (X(Z,1)+X(Z,2)+X(Z,3))/3.D0 40 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 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.D0 F1UX(Z,2,2)=1.D0 F1UX(Z,3,3)=1.D0 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.1).AND.(COMPU.EQ.4).AND.(CLASS.EQ.3)) THEN DO 16 Z=1,NELIS F1U(Z,1)=1.D0 16 CONTINUE ENDIF IF ((COMPV.EQ.2).AND.(COMPU.EQ.1).AND.(CLASS.EQ.3)) THEN DO 20 Z=1,NELIS F0U(Z)=RE*DUDX(Z,2,1) 20 CONTINUE ENDIF IF ((COMPV.EQ.2).AND.(COMPU.EQ.2).AND.(CLASS.EQ.3)) THEN DO 24 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.D0 F1UX(Z,2,2)=1.D0 F1UX(Z,3,3)=1.D0 24 CONTINUE ENDIF IF ((COMPV.EQ.2).AND.(COMPU.EQ.3).AND.(CLASS.EQ.3)) THEN DO 28 Z=1,NELIS F0U(Z)=RE*DUDX(Z,2,3) 28 CONTINUE ENDIF IF ((COMPV.EQ.2).AND.(COMPU.EQ.4).AND.(CLASS.EQ.3)) THEN DO 32 Z=1,NELIS F1U(Z,2)=1.D0 32 CONTINUE ENDIF IF ((COMPV.EQ.3).AND.(COMPU.EQ.1).AND.(CLASS.EQ.3)) THEN DO 36 Z=1,NELIS F0U(Z)=RE*DUDX(Z,3,1) 36 CONTINUE ENDIF IF ((COMPV.EQ.3).AND.(COMPU.EQ.2).AND.(CLASS.EQ.3)) THEN DO 40 Z=1,NELIS F0U(Z)=RE*DUDX(Z,3,2) 40 CONTINUE ENDIF IF ((COMPV.EQ.3).AND.(COMPU.EQ.3).AND.(CLASS.EQ.3)) THEN DO 44 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.D0 F1UX(Z,2,2)=1.D0 F1UX(Z,3,3)=1.D0 44 CONTINUE ENDIF IF ((COMPV.EQ.3).AND.(COMPU.EQ.4).AND.(CLASS.EQ.3)) THEN DO 48 Z=1,NELIS F1U(Z,3)=1.D0 48 CONTINUE ENDIF IF ((COMPV.EQ.4).AND.(COMPU.EQ.1).AND.(CLASS.EQ.3)) THEN DO 52 Z=1,NELIS F0UX(Z,1)=1.D0 52 CONTINUE ENDIF IF ((COMPV.EQ.4).AND.(COMPU.EQ.2).AND.(CLASS.EQ.3)) THEN DO 56 Z=1,NELIS F0UX(Z,2)=1.D0 56 CONTINUE ENDIF IF ((COMPV.EQ.4).AND.(COMPU.EQ.3).AND.(CLASS.EQ.3)) THEN DO 60 Z=1,NELIS F0UX(Z,3)=1.D0 60 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)+U(Z,4) 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)) & +1.D0/3.D0-RE*(X(Z,1)*(X(Z,3)**2+X(Z,2)**2)- & X(Z,1)**2*(X(Z,2)+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)+U(Z,4) 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)) & +1.D0/3.D0-RE*(X(Z,2)*(X(Z,1)**2+X(Z,3)**2)- & X(Z,2)**2*(X(Z,3)+X(Z,1))) 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)+U(Z,4) F0(Z)=RE*(U(Z,1)*DUDX(Z,3,1)+ & U(Z,2)*DUDX(Z,3,2)+U(Z,3)*DUDX(Z,3,3)) & +1.D0/3.D0-RE*(X(Z,3)*(X(Z,1)**2+X(Z,2)**2)- & X(Z,3)**2*(X(Z,1)+X(Z,2))) 12 CONTINUE ENDIF IF (COMPV.EQ.4) THEN DO 13 Z=1,NELIS F0(Z)=DUDX(Z,1,1)+DUDX(Z,2,2)+DUDX(Z,3,3) 13 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.0.0001D0) F0(Z)=-X(Z,1) IF (X(Z,2).GT.0.9999D0) F0(Z)= X(Z,1) IF (X(Z,3).LT.0.0001D0) F0(Z)= X(Z,1) IF (X(Z,3).GT.0.9999D0) F0(Z)=-X(Z,1) 3 CONTINUE ENDIF IF (COMPV.EQ.2) THEN DO 7 Z=1,NELIS IF (X(Z,1).LT.0.0001D0) F0(Z)= X(Z,2) IF (X(Z,1).GT.0.9999D0) F0(Z)=-X(Z,2) IF (X(Z,3).LT.0.0001D0) F0(Z)=-X(Z,2) IF (X(Z,3).GT.0.9999D0) F0(Z)= X(Z,2) 7 CONTINUE ENDIF IF (COMPV.EQ.3) THEN DO 11 Z=1,NELIS IF (X(Z,1).LT.0.0001D0) F0(Z)=-X(Z,3) IF (X(Z,1).GT.0.9999D0) F0(Z)= X(Z,3) IF (X(Z,2).LT.0.0001D0) F0(Z)= X(Z,3) IF (X(Z,2).GT.0.9999D0) F0(Z)=-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 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** IF (COMPU.EQ.1) THEN DO 10 Z=1,NDC B(Z) = X(Z,1)*(X(Z,3)-X(Z,2)) 10 CONTINUE ENDIF IF (COMPU.EQ.2) THEN DO 20 Z=1,NDC B(Z) = X(Z,2)*(X(Z,1)-X(Z,3)) 20 CONTINUE ENDIF IF (COMPU.EQ.3) THEN DO 30 Z=1,NDC B(Z) = X(Z,3)*(X(Z,2)-X(Z,1)) 30 CONTINUE ENDIF IF (COMPU.EQ.4) THEN DO 40 Z=1,NDC B(Z) = (X(Z,1)+X(Z,2)+X(Z,3))/3.D0 40 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 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,1)*(X(Z,3)-X(Z,2))) CU(Z,2) = ABS( U(Z,2) - X(Z,2)*(X(Z,1)-X(Z,3))) CU(Z,3) = ABS( U(Z,3) - X(Z,3)*(X(Z,2)-X(Z,1))) CU(Z,4) = ABS( U(Z,4) - (X(Z,1)+X(Z,2)+X(Z,3))/3.D0) 10 CONTINUE C** C**----------------------------------------------------------------- C** C**** end of calculation C** ------------------ C** R E T U R N C**---end of USERC-------------------------------------------------- E N D |