D:/demona/2.0/manual/doc/2-31set/problems/test/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) Tc,Tw,Te,Ts,Tn

real(8) T1,T2,T3,T4

real(8) d2TdX2,d2TdY2



real(8) dtau

real(8) upu,upv

real(8) myU,myV



real(8) sourceT
real(8) sourceU

real(8) dx,dy



!< 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






!stop




discrete=-1.d0


x=>discreteY


! 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
modelDomain=1


sourceT=1.d0






dx=1.d0/(myModel%physics(modelPhysics)%domain(modelDomain)%maxI-myModel%physics(modelPhysics)%domain(modelDomain)%minI)
dy=1.d0/(myModel%physics(modelPhysics)%domain(modelDomain)%maxJ-myModel%physics(modelPhysics)%domain(modelDomain)%minJ)




!< 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





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


dof=>myModel%physics(modelPhysics)%dof

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

IM2=>myModel%physics(2)%domain(1)%dofMatrix


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

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


do myP=1,myModel%endIndex(modelPhysics)


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


!boundary

!---------------------------------------------------------------------------------------------
!bottom
!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
        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
                        TC=x(IM(i,j,1,absolute))

                        !print*,TC,i,j

                        RHS(t)=TC-0.d0

                end do



                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
!---------------------------------------------------------------------------------------------
!top
!bc=3
myBC=3 

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
        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
                        TC=x(IM(i,j,1,absolute))

                        !print*,TC,i,j


                        ! here the coupling is defined explicitly
                        ! need and interpolation interface
                        ! such as:
                        !
                        ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
                

                        T1=x(IM2(1,j,1,absolute))
                        T2=x(IM2(1+1,j,1,absolute))
                        T3=x(IM2(1+1,j+1,1,absolute))
                        T4=x(IM2(1,j+1,1,absolute))

                        !interpolation
                        RHS(t)=TC-(T1+T2+T3+T4)/4.d0

                end do


                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
        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)


                        TC=x(IM(i,j,1,absolute))
                        TW=x(IM(i-1,j,1,absolute))
                        TE=x(IM(i+1,j,1,absolute))
                        TS=x(IM(i,j-1,1,absolute))
                        TN=x(IM(i,j+1,1,absolute))

                        d2TdX2=(TE+TW-2.d0*TC)/(dx**2.d0)

                        d2TdY2=(TN+TS-2.d0*TC)/(dy**2.d0)

                        RHS(t)=d2TdX2+d2TdY2+sourceT

                        !RHS(t)=TC

                        !print*,TC,i,j

                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
!---------------------------------------------------------------------------------------------
!stop

end do ! domain 1 loop

end do ! physics loop

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






deallocate(RHS)




!x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
!x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
!x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x
!x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x-x

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


modelPhysics=2


modelDomain=1



dof=>myModel%physics(modelPhysics)%dof

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

IM2=>myModel%physics(1)%domain(1)%dofMatrix

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


! these declarations are needed for physics one where the location
! of other physics are @ i=1, here on the other hand, the locatiob
! is nx and nx-1 so nx will be called.
nx=>myModel%physics(1)%domain(1)%nx
ny=>myModel%physics(1)%domain(1)%ny


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

do myP=1,myModel%endIndex(modelPhysics)


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

!boundary

!---------------------------------------------------------------------------------------------
!bottom
!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
        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
                        TC=x(IM(i,j,1,absolute))

                        RHS(t)=TC-0.d0

                end do


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

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


                end do

        end do


end do
!---------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------
!top
!bc=3
myBC=3 

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
        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
                        TC=x(IM(i,j,1,absolute))


                        ! here the coupling is defined explicitly
                        ! need and interpolation interface
                        ! such as:
                        !
                        ! TC - interpolate(TC_related_info_i,j,t,physics,domain)
                

                        T1=x(IM2(nx-1,j-1,1,absolute))
                        T2=x(IM2(nx-1+1,j-1,1,absolute))
                        T3=x(IM2(nx-1+1,j,1,absolute))
                        T4=x(IM2(nx-1,j,1,absolute))

                        !interpolation
                        RHS(t)=TC-(T1+T2+T3+T4)/4.d0

                end do


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

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


                end do

        end do


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

!---------------------------------------------------------------------------------------------
! 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
        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)
                        
                        TC=x(IM(i,j,1,absolute))
                        TW=x(IM(i-1,j,1,absolute))
                        TE=x(IM(i+1,j,1,absolute))
                        TS=x(IM(i,j-1,1,absolute))
                        TN=x(IM(i,j+1,1,absolute))

                        d2TdX2=(TE+TW-2.d0*TC)/(dx**2.d0)

                        d2TdY2=(TN+TS-2.d0*TC)/(dy**2.d0)

                        RHS(t)=d2TdX2+d2TdY2+sourceT

                        !RHS(t)=TC


                end do



                do myT=1,maxT

                        t=tVec(myT)

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

                end do

        end do


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

end do


end do


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













deallocate(RHS)



discreteY=>y

end function discrete












end module model

---