[Nek5000-users] Gradients in an object : Few queries for parallel setting.

[nek5000-users at lists.mcs.anl.gov]

Shriram,

I think the following code should work in parallel, though I've
not tested it.

Note that I was unsure of your desired logic wrt how often you
wanted to write the data etc, so you'll need to modify accordingly.

Paul



       common /scrns/         pm1(lx1,ly1,lz1,lelv)
      $,vxx(lx1,ly1,lz1,lelv),vxy(lx1,ly1,lz1,lelv),vxz(lx1,ly1,lz1,lelv)
      $,vyx(lx1,ly1,lz1,lelv),vyy(lx1,ly1,lz1,lelv),vyz(lx1,ly1,lz1,lelv)
       common /scruz/
      $ vzx(lx1,ly1,lz1,lelv),vzy(lx1,ly1,lz1,lelv),vzz(lx1,ly1,lz1,lelv)
      $,temp(4,ly1**2),work(4,ly1**2)

       if (istep.eq.0) call set_obj  ! objects for surface integrals


       if (nid.eq.0.and.istep.eq.0) then   ! ONLY NODE 0 opens file
          open(unit=103,file='coord.txt') ! few files for post-processing
          open(unit=104,file='elenos.txt')
       endif

       call gradm1 (vxx,vxy,vxz,vx)   ! MORE EFFICIENT than dudxyz


       DO 100 II=1,NHIS

          IF (HCODE(10,II).NE.'I') GOTO 100
          IOBJ   = LOCHIS(1,II)
          MEMTOT = NMEMBER(IOBJ)

          IF (HCODE(1,II).NE.' ' .OR. HCODE(2,II).NE.' ' .OR.
     $        HCODE(3,II).NE.' ' ) THEN
            IFIELD = 1

            DO 50 MEM=1,MEMTOT
               ISK   = 0
               IEG   = OBJECT(IOBJ,MEM,1)
               IFC   = OBJECT(IOBJ,MEM,2)
               if (nid.eq.0) write(104,*) ieg,gllel(ieg),gllnid(ieg)

               call rzero(temp,4*ly1*ly1)

               if (gllnid(ieg).eq.nid) then ! this processor has a contribution

                  ie = gllel(ieg)

                  CALL DSSET(NX1,NY1,NZ1)
                  IFACE  = EFACE1(IFC)
                  JS1    = SKPDAT(1,IFACE)
                  JF1    = SKPDAT(2,IFACE)
                  JSKIP1 = SKPDAT(3,IFACE)
                  JS2    = SKPDAT(4,IFACE)
                  JF2    = SKPDAT(5,IFACE)
                  JSKIP2 = SKPDAT(6,IFACE)

                  I = 0
                  DO 110 J2=JS2,JF2,JSKIP2
                  DO 110 J1=JS1,JF1,JSKIP1
                     I = I+1
                     temp(1,i) = xm1(j1,j2,1,ie)
                     temp(2,i) = xm1(j1,j2,1,ie)
                     temp(3,i) = xm1(j1,j2,1,ie)
                     temp(4,i) = vxy(j1,j2,1,ie)
  110             CONTINUE

               endif
               call gop(temp,work,'+  ',4*ly1*ly1)  ! gather results onto rank 0

               if (nid.eq.0) write(103,18) ((temp(k,i),k=1,4),i=1,ly1*ly1)
   18          format(1p4e15.7)

   50       CONTINUE
          ENDIF
       ENDDO




On Thu, 10 Jun 2010, nek5000-users at lists.mcs.anl.gov wrote:

> Hello,
>
> I am trying to find gradients and the coordinates of an object in nek. I
> used set_obj routine to fix my walls as object and also wrote a routine that
> is
> similar to the drag_calc routine to dump out the co-ordinates and the
> gradients on the object. This code works fine for a serial setting while for
> parallel setting,
> only a part of the object coordinates are dumped out. I suspect that I have
> not included a few routines that would gather information from all the
> processors.
> It would be great if you could throw in some suggestions that would set it
> to work in parallel setting. I have pasted a part of the usrchk routine
> herewith :
>
> <code>
>
>       DIMENSION TEMP1(lz1**2) ! Stores the du/dy of all the GLL points in
> one face (in the object)
>
>       parameter(lt=lx1*ly1*lz1)
>       common /scrns/         pm1(lx1,ly1,lz1,lelv)
>     $,vxx(lx1,ly1,lz1,lelv),vxy(lx1,ly1,lz1,lelv),vxz(lx1,ly1,lz1,lelv)
>     $,vyx(lx1,ly1,lz1,lelv),vyy(lx1,ly1,lz1,lelv),vyz(lx1,ly1,lz1,lelv)
>       common /scruz/
>     $ vzx(lx1,ly1,lz1,lelv),vzy(lx1,ly1,lz1,lelv),vzz(lx1,ly1,lz1,lelv)
>     $,one(lx1,ly1,lz1,lelv)
>
>      if (istep.eq.0) then
>         call set_obj                   ! objects for surface integrals
>      endif
>
>       call izero(TEMP1,lz1**2) ! Initializing to zeros.
>
>        open(unit=103,file='coord.txt') ! few files for post-processing
>        open(unit=104,file='elenos.txt')
>
>       CALL DUDXYZ (vxx,vx,RXM1,SXM1,TXM1,JACM1,1,1)
>       CALL DUDXYZ (vxy,vx,RYM1,SYM1,TYM1,JACM1,1,1)
>       CALL DUDXYZ (vxz,vx,RZM1,SZM1,TZM1,JACM1,1,1)
>
>       ntot1 = nx1*ny1*nz1*nelv
>
>       call rone (one,ntot1)
> c------------------------------------------------------------------
>
>      DO 100 II=1,NHIS
>
>         IF (HCODE(10,II).NE.'I') GOTO 100
>         IOBJ   = LOCHIS(1,II)
>         MEMTOT = NMEMBER(IOBJ)
>
>         IF (HCODE(1,II).NE.' ' .OR. HCODE(2,II).NE.' ' .OR.
>     $       HCODE(3,II).NE.' ' ) THEN
>            IFIELD = 1
>
>            DO 50 MEM=1,MEMTOT
>               ISK   = 0
>               IEG   = OBJECT(IOBJ,MEM,1)
>               IFC   = OBJECT(IOBJ,MEM,2)
> c               IF (GLLNID(IEG).EQ.NID) THEN
> C                 This processor has a contribution
>                  IE = GLLEL(IEG)
>                   write(104,*)ieg,ie,gllnid(ieg),nid
>
>       call rone (one,ntot1)
>       call dotabf(TEMP1,vxy,one,ifc,ie) ! subroutine that multiplies vxy by
> one and stores in temp for that element
>
>      CALL DSSET(NX1,NY1,NZ1)
>      IFACE  = EFACE1(IFC)
>      JS1    = SKPDAT(1,IFACE)
>      JF1    = SKPDAT(2,IFACE)
>      JSKIP1 = SKPDAT(3,IFACE)
>      JS2    = SKPDAT(4,IFACE)
>      JF2    = SKPDAT(5,IFACE)
>      JSKIP2 = SKPDAT(6,IFACE)
>
>      I = 0
>      DO 110 J2=JS2,JF2,JSKIP2
>      DO 110 J1=JS1,JF1,JSKIP1
>         I = I+1
>         if(istep.eq.2) then
>         write(103,18)xm1(J1,J2,1,ie),ym1(J1,J2,1,ie),zm1(J1,J2,1,ie)
>     $ ,TEMP1(I)
>   18  format(1x,4(E14.7,1X))
>         endif
>  110 CONTINUE
>
> c            ENDIF
>   50       CONTINUE
>
>
> </code>
>
> If the rea/usr file is needed, please let me know. Thanks.
>
> Thanks
> Shriram
>