D:/demona/2.0/manual/doc/2-31set/problems/bratu/1/model.f90

Go to the documentation of this file.

module model
!module modelSet  ! use to differentiate between set and no set models.


use global

use myFunctions

implicit none


contains 


function discrete(a)


real(8) discrete(globalCurrentSetLength)

integer i,j,t

integer a

integer myI,myJ,myT,myBC,myUnk

integer myP,myD

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

integer,pointer :: IM2(:,:,:,:)

integer,pointer :: maxI,maxJ,maxT

integer,pointer :: iVec(:),jVec(:),bcVec(:),tVec(:)

integer,pointer :: endIndex(:)

integer modelPhysics
integer modelDomain

real(8), allocatable :: RHS(:)

integer,pointer :: dof

integer,pointer :: nx,ny

!integer,pointer :: dof


! model specific definitions

real(8) STFc,STFw,STFe,STFs,STFn,stfcpto
real(8) VORc,VORw,VORe,VORs,VORn,vorcpto

real(8) VORCOREAL

real(8) dSTFdX,dSTFdY
real(8) dVORdX,dVORdY
real(8) d2STFdX2,d2STFdY2
real(8) d2VORdX2,d2VORdY2

real(8) dvordtau,dstfdtau


real(8) UdVORdX,VdVORdY


!real(8) dtau

real(8) upu,upv

real(8) myU,myV



real(8) sourceT
real(8) sourceU

real(8),pointer :: lambda

real(8) dx,dy

real(8),pointer :: myLx,myLy

real(8),pointer :: Re


!< if x+eps*v is to be evaluated, then tempy is updated before calling discrete(1)
!< such that tempy(setVec)=tempy(setVec)+lambda*deltax
!< if discrete(0) is called, original y is evaluated
do i=1,a
        discreteY=>tempy
end do








discrete=0.d0


x=>discreteY

xPTO=>yPTO


discrete=x(setVEc)



! In this problem, discretization for 1 and 2 is the same although they are located in different
! part of the domains of different physics
! they might be coded in a loop so double explicit definition for P1 and P2 is avoided
! this is not the case for BC 3 since it is defined explicitly
! however, if and interpolation interface is used, than the bourden is passed to the interface routine
! so that the discretizations looks the same. After then, P1 and P2 can be defined in a loop.








!!!!!!!!!!!!!!
!< assume that  physics-1 starts here
!< also that the domain-1 of physics-1 starts here, as well
!#############################################################################################
!#############################################################################################
!#############################################################################################
!
!############                                   ###
!#############                                  ###
!###            ###                                     ###
!###            ###                                     ###
!###            ###                                     ###
!###            ###                                     ###
!#############          ####            ###
!############           ####            ###
!###                                                    ###
!###                                                    ###
!###                                                    ###
!###                                                    ###
!###                                                    ###
!###                                                    ###     
!
!#############################################################################################
!#############################################################################################
!#############################################################################################


!modelPhysics=currentPhysics
!modelDomain=currentDomain

modelPhysics=1


!< the idea is  mymodel%discrete%set(1)
!<                              mymodel%discrete%set(2)
!<                              mymodel%discrete%set(3)
! with set vec pointing to the set of the physics(..)%domain%(..)%set


lambda=>myModel%physics(modelPhysics)%physicalParameters(1)



!Re=100.d0
!Re=1000.d0
!Re=2000.d0
!Re=3000.d0
!Re=3500.d0
!Re=5000.d0
!Re=10000.d0
!stop



!print*,'dx,dy=',dx,dy
!stop


dof=>myModel%physics(modelPhysics)%dof




allocate(RHS(dof))
RHS=0.d0

!#####################################################################
!user supplied variables
!#####################################################################

do myP=1,myModel%endIndex(modelPhysics)

        !#############################################################################################
        !
        !########                       ##
        !##             ##                      ##
        !##             ##                      ##
        !##             ##      ###             ##
        !##             ##                      ##
        !##             ##                      ##
        !########                       ##
        !
        !#############################################################################################
        
        modelDomain=1

        IM=>myModel%physics(modelPhysics)%domain(modelDomain)%dofMatrix

        myLx=>myModel%physics(modelPhysics)%domain(modelDomain)%Lx
        myLy=>myModel%physics(modelPhysics)%domain(modelDomain)%Ly

        dx=myLx/(myModel%physics(modelPhysics)%domain(modelDomain)%maxI-myModel%physics(modelPhysics)%domain(modelDomain)%minI)
        dy=myLy/(myModel%physics(modelPhysics)%domain(modelDomain)%maxJ-myModel%physics(modelPhysics)%domain(modelDomain)%minJ)
        ! normally, especially in finite volume coordinates should be read so that face lengths and center-to-center distances
        ! are calculated. for this uniform finite difference grid dx,dy is enough to proceed

        ! caution, how is dx dy in finite volume!!! was lx for the whole domain with ghostcells 
        ! or just the real interior domain w/o ghost cells!!! - if so number of ghostcells connection!!

        !print*,dx,dy



        do myD=1,myModel%physics(modelPhysics)%endIndex(modelDomain)


        !boundary

        !---------------------------------------------------------------------------------------------
        !left
        !bc=2
        myBC=2   

        maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
        jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec


        do myJ=1,maxJ
                
                j=jVec(myJ)
                maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
                iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec

                RHS=0.d0
                do myI=1,maxI

                        i=iVec(myI)
                        
                        maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
                        tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
                        endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex

                        !unk1
                        t=1
                        do myT=1,endIndex(t)

                                ! center values
                                STFC=x(IM(i,j,t,absolute)) !< or TC=x(IM(i,j,*1*,absolute))

                                !print*,TC,i,j

                                RHS(t)=STFC-0.d0

                        end do

                        !unk1


                        do myT=1,maxT
                                
                                t=tVec(myT)

                                !print*,IM(i,j,t,relative)

                                discrete(IM(i,j,t,relative))=rhs(t)

                        end do

                end do



        end do
        !---------------------------------------------------------------------------------------------
        !stop
        !---------------------------------------------------------------------------------------------
        !interior
        !bc=1
        myBC=1

        maxJ=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%maxJ
        jVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%jVec

        do myJ=1,maxJ
                
                j=jVec(myJ)
                maxI=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%maxI
                iVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%iVec

                RHS=0.d0
                do myI=1,maxI

                        i=iVec(myI)
                        
                        maxT=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%maxT
                        tVec=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%tVec
                        endIndex=>myModel%physics(modelPhysics)%domain(modelDomain)%set%bc(myBC)%j(myJ)%i(myI)%endIndex

                        !unk1
                        t=1
                        do myT=1,endIndex(t)


                                STFC=x(IM(i,j,t,absolute))
                                STFW=x(IM(i-1,j,t,absolute))
                                STFE=x(IM(i+1,j,t,absolute))
                                STFS=x(IM(i,j-1,t,absolute))
                                STFN=x(IM(i,j+1,t,absolute))

                                d2STFdX2=(STFW+STFE-2.d0*STFC)/(dx**2.d0)

                                d2STFdY2=(STFN+STFS-2.d0*STFC)/(dy**2.d0)

                                RHS(t)=d2STFdX2+d2STFdY2+lambda*dexp(STFC)





                        end do




                        do myT=1,maxT

                                t=tVec(myT)

                                !print*,IM(i,j,t,relative),rhs(t)

                                discrete(IM(i,j,t,relative))=rhs(t)

                        end do

                end do


        end do
        !---------------------------------------------------------------------------------------------

        end do ! domain 1 loop




end do ! physics loop

!#####################################################################


deallocate(RHS)






!===================================================================
discreteY=>y

end function discrete












end module model

---