D:/demona/2.0/manual/doc/2-31set/gs.f90

Go to the documentation of this file.

subroutine gs

use global

use model


implicit none


integer myP,myD

integer,pointer :: IM(:,:,:,:)

integer,pointer :: IV(:,:)   ! absolute,i,j,t

integer,pointer :: gsSetVec(:)

integer maxNewtonBackup

integer i

integer gsIter

integer gsTol

real(8) normGS

integer myIndex

integer myI,myJ,myT,myN

integer mySize

integer,allocatable :: myPVec(:),myDvec(:)

integer gsCurrentPhysics,gsCurrentDomain


maxNewtonBackup=maxNewtonIter

maxNewtonIter=gsMaxNewton

gsTol=1e-5


pVec=>myModel%pVec


if (accumulating.eqv.yes) maxGsIter=1



!find a better way to avoid callin current set
! maybe two sets can be used. one globalset - which is permanent - like the domain
! and second set for just one unk



! zero all flag layers in the global set
! copy localsetvecs i.e. global set vec in a sense.
do myP=1,size(pVec)

        currentPhysics=pVec(myP)

        dVec=>myModel%physics(currentPhysics)%dVec

        do myD=1,size(dVec)

                currentDomain=dVec(myD)

                IM=>myModel%physics(currentPhysics)%domain(currentDomain)%dofMatrix
                
                IM(:,:,:,flag)=0

                mySize=size(myModel%physics(currentPhysics)%domain(currentDomain)%setVec)

                allocate(myModel%physics(currentPhysics)%domain(currentDomain)%gsSetVec(mySize))
                gsSetVec=>myModel%physics(currentPhysics)%domain(currentDomain)%gsSetVec
                gsSetVEc=myModel%physics(currentPhysics)%domain(currentDomain)%setVec

        end do

end do



allocate(myPvec(size(pVec)))
myPVec=pVEc

do myP=1,size(mypVec)

        gsCurrentPhysics=myPVec(myP)

        dVec=>myModel%physics(gsCurrentPhysics)%dVec

        allocate(myDvec(size(dVec)))
        myDVec=dVEc

        do myD=1,size(mydVec)

                !print*,myD

                gsCurrentDomain=myDVec(myD)

                IM=>myModel%physics(gsCurrentPhysics)%domain(gsCurrentDomain)%dofMatrix
                IV=>myModel%physics(gsCurrentPhysics)%domain(gsCurrentDomain)%dofVector

                relativeIndex=>myModel%physics(gsCurrentPhysics)%domain(gsCurrentDomain)%relativeIndex
        
                !allocate(gsSetVec(size(myModel%physics(currentPhysics)%domain(currentDomain)%setVec)))

                gsSetVec=>myModel%physics(gsCurrentPhysics)%domain(gsCurrentDomain)%gsSetVec


                do gsIter=1,maxGsIter
                        

                        if (mod(gsIter,10).eq.0) then

                        print*,'^^^^^^^^^^^^^^^^^^^^^'
                        print*,'gauss seidel iteration=',gsIter
                        print*,'^^^^^^^^^^^^^^^^^^^^^'

                        end if

                        call postprocessor

                        print*,'if block is on, then gsSetVec should be smaller'

                        IM(:,:,:,flag)=0

                        do i=1,size(gsSetVec)
                                
                                
                                myIndex=gsSetVec(i)

                                myIndex=myIndex-relativeIndex

                                myN=IV(myIndex,1)
                                myI=IV(myIndex,2)
                                myJ=IV(myIndex,3)
                                myT=IV(myIndex,4)

                                IM(myI,myJ,myT,flag)=1

                                if (blockGS.eqv.yes) then

                                        IM(myI,myJ,:,flag)=1

                                end if



                                call currentSet

                                !print*,setVec
                                !stop



                                call Newton


                                !print*,y(11)

                                !pause

                                
                                
                                
                                !IM(myI,myJ,myT,flag)=1 ! try also accumulating version?
                                if (accumulating.eqv.no) IM(myI,myJ,myT,flag)=0

                                if (blockGS.eqv.yes) then

                                        !IM(myI,myJ,:,flag)=1 ! try also accumulating version?
                                        if (accumulating.eqv.no) IM(myI,myJ,:,flag)=0

                                end if



                        end do

                        !stop
                        
                        !print*,'---------------------------------------------------------------------------'
                        !print*,y(9)
                        !print*,'---------------------------------------------------------------------------'

                        

                        normGS=dsqrt(dot_product(discrete(0),discrete(0)))



                        !if (normGS.le.gsTol) then

                                !print*,'gauss seidel converged at iter=',gsIter
                                !print*,'norm is',normGS
                                !maxNewtonIter=maxNewtonBackup
                                !return

                        !end if


                end do

                !stop


                
                !stop

                nullify(gsSetVec)

        end do

        deallocate(myDVec)


end do


! zero all flag layers in the global set
! copy localsetvecs i.e. global set vec in a sense.
do myP=1,size(pVec)

        currentPhysics=pVec(myP)

        dVec=>myModel%physics(currentPhysics)%dVec

        do myD=1,size(dVec)

                currentDomain=dVec(myD)

                IM=>myModel%physics(currentPhysics)%domain(currentDomain)%dofMatrix
                
                IM(:,:,:,flag)=1

                deallocate(myModel%physics(currentPhysics)%domain(currentDomain)%gsSetVec)


        end do

end do

call currentset

maxNewtonIter=maxNewtonBackup





deallocate(myPVec)



















end subroutine gs

---