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

Go to the documentation of this file.

module myFunctions

implicit none



contains



!#################################################################################
function distance(p1,p2)
implicit none

real distance

real(8) p1(:),p2(:)


distance=dsqrt((p1(1)-p2(1))**2.d0+(p1(2)-p2(2))**2.d0+(p1(3)-p2(3))**2.d0)




end function distance

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


function get_ptr_lb(array,lb1) result(ptr)

integer,intent(in):: lb1

integer,dimension(lb1:),intent(in),target:: array

integer,dimension(:),pointer:: ptr
ptr => array


!! usage
!integer,pointer :: aa(:)


!integer,target:: bb(2)


!bb(1)=11
!bb(2)=12

!aa => get_ptr_lb(bb,2)

!print*,aa(3)


end function

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



function getElementNodes(myNx,myNy,myNz) result(info)

implicit none

integer myNx,myNy,myNz

integer info(2)

integer numberOfElementalNodes

integer axis

        select case(myNz)

                case(1)

                        select case(myNy)

                                case(1)

                                        select case (myNx)

                                                case(1)

                                                        print*,'no grid point !'
                                                        
                                                        print*,'demona stop'

                                                        stop

                                                        axis=0

                                                case default
                                                
                                                        numberOfElementalNodes=2
                                                        axis=1

                                        end select


                                case default

                                        numberOfElementalNodes=4
                                        axis=2

                        end select


                case default
        
                        numberOfElementalNodes=8

                        axis=3

        
        end select


info(1)=numberOfElementalNodes
info(2)=axis





end function getElementNodes






!#################################################################################
FUNCTION reallocate_vector(p, n)               ! reallocate REAL
        REAL(8), POINTER, DIMENSION(:) :: p, reallocate_vector
        INTEGER, intent(in) :: n
        INTEGER :: nold, ierr

        ALLOCATE(reallocate_vector(1:n), STAT=ierr)

        REALLOCATE_vector=0.d0

        IF(ierr /= 0) STOP "allocate error"
        IF(.NOT. ASSOCIATED(p)) RETURN

        nold = MIN(SIZE(p), n)
        reallocate_vector(1:nold) = p(1:nold)

        DEALLOCATE(p) 

END FUNCTION REALLOCATE_vector
!#################################################################################
FUNCTION reallocate_vector_integer(p, n)               ! reallocate REAL
        INTEGER, POINTER, DIMENSION(:) :: p, reallocate_vector_integer
        INTEGER, intent(in) :: n
        INTEGER :: nold, ierr

        ALLOCATE(reallocate_vector_integer(1:n), STAT=ierr)

        reallocate_vector_integer=0

        IF(ierr /= 0) STOP "allocate error"
        IF(.NOT. ASSOCIATED(p)) RETURN

        nold = MIN(SIZE(p), n)
        reallocate_vector_integer(1:nold) = p(1:nold)

        DEALLOCATE(p) 

END FUNCTION REALLOCATE_vector_integer
!#################################################################################
FUNCTION reallocate_matrix(p,n,m)               ! reallocate REAL
        REAL(8), POINTER, DIMENSION(:,:) :: p, reallocate_matrix
        INTEGER, intent(in) :: n,m
        INTEGER :: nold,mold,ierr

        ALLOCATE(reallocate_matrix(n,m), STAT=ierr)

        reallocate_matrix=0.d0

        IF(ierr /= 0) STOP "allocate error"
        IF(.NOT. ASSOCIATED(p)) RETURN

        nold = SIZE(p,DIM=1)
        mold = SIZE(p,DIM=2)
        
        reallocate_matrix(1:nold,1:mold) = p(1:nold,1:mold)


        DEALLOCATE(p) 

END FUNCTION REALLOCATE_matrix
!#################################################################################
FUNCTION reallocate_matrix_integer(p,n,m)               ! reallocate REAL
        INTEGER, POINTER, DIMENSION(:,:) :: p, reallocate_matrix_integer
        INTEGER, intent(in) :: n,m
        INTEGER :: nold,mold,ierr

        ALLOCATE(reallocate_matrix_integer(n,m), STAT=ierr)

        reallocate_matrix_integer=0

        IF(ierr /= 0) STOP "allocate error"
        IF(.NOT. ASSOCIATED(p)) RETURN

        nold = SIZE(p,DIM=1)
        mold = SIZE(p,DIM=2)
        
        reallocate_matrix_integer(1:nold,1:mold) = p(1:nold,1:mold)


        DEALLOCATE(p) 

END FUNCTION REALLOCATE_matrix_integer
!#################################################################################
function multiplier(R,sint,cost,j)
        implicit none
        integer i,j
        real(8) sint,cost
        real(8),dimension(:,:) ::R
        real(8) multiplier(size(R,DIM=1),size(R,DIM=2))
        j=size(R,DIM=2)
        multiplier=R
        multiplier(j,j)=cost*R(j,j)+sint*R(j+1,j)
        multiplier(j+1,j)=-sint*R(j,j)+cost*R(j+1,j)
end function multiplier
!#################################################################################
function previous(H,sint,cost)
        implicit none
        integer i,j
        real(8),dimension(:) :: sint,cost
        real(8),dimension(:) ::H
        real(8) :: previous(size(H))
        real(8) aux
        j=size(sint)
        previous=H
        ! ilk bakista burada aux gerekmiyor gibi duruyor cunku sol va sag taraflar farklýsol:previous,sag H
        ! ama buradaki transformationlar arka arkaya uygulanmali, bu yüzden degisen bir deger tekrar degismeli
        do i=1,j-1
                aux=previous(i)
                previous(i)=cost(i)*previous(i)+sint(i)*previous(i+1)
                previous(i+1)=-sint(i)*aux+cost(i)*previous(i+1)
        end do
end function previous
!#################################################################################
function backwardsub(U,b)
        implicit none
        real(8)  U(:,:),b(:)
        real(8)  backwardsub(size(b))
        real(8)  c
        integer i,j,n
        backwardsub=0.d0
        n=size(b)
        backwardsub(n)=b(n)/U(n,n)
        do i=n-1,1,-1
                c=0.d0
                do j=i+1,n
                        c=backwardsub(j)*U(i,j)+c
                end do
                backwardsub(i)=(b(i)-c)/U(i,i)
        end do
end function backwardsub
!#################################################################################
function scaler(a,b)
implicit none
integer i,n
real(8) a(:),b(:)
real(8) scaler(size(a))
n=size(a)
scaler=0.d0
!$OMP PARALLEL DO
do i=1,n        
        scaler(i)=a(i)*b(i)
end do
!$OMP END PARALLEL DO

end function scaler


!#################################################################################
function invscaler(a,b)
implicit none
integer i,n
real(8) a(:),b(:)
real(8) invscaler(size(a))
n=size(a)
invscaler=0.d0
!$OMP PARALLEL DO
do i=1,n        
        invscaler(i)=a(i)/b(i)
end do
!$OMP END PARALLEL DO

end function invscaler


!#################################################################################
function signum(x)
implicit none
real(8) signum,eps,x


eps=epsilon(0.d0)

signum=x/dabs(x+eps)


end function signum
!#################################################################################
function wNorm(myVec,myWVec)

implicit none 

real(8) myVec(:),myWVec(:)

real(8) wNorm

wNorm=0.d0

myVec=invscaler(myVec,myWvec)

wNorm=dsqrt(dot_product(myVec,myVec))

end function wNorm
!#################################################################################


subroutine fileNamer(phyNum,domNum,file,timeNum)

! also consider time dependence.

implicit none
character(len=1) :: phy
character(len=1) :: dom
character(len=36) :: file
character(len=4) :: ext
character(len=2) :: phyStr
character(len=2) :: domStr
character(len=15) :: timeStr
character(len=8) :: folderStr



integer phyNum
integer domNum
integer timeNum

WRITE(phyStr,'(I2)') phyNum
WRITE(domStr,'(I2)') domNum
WRITE(timeStr,'(I15)') timeNum

folderStr='results/'
phy='P'
dom='D'
ext='.dat'
file=folderStr//phy//phyStr//' - '//dom//domStr//timeStr//ext

end subroutine fileNamer



!#################################################################################
function cubic(df0,f0,fprev1,fprev2,lambdaprev1,lambdaprev2) !<cubic backtracking
implicit none
real(8) df0,f0,fprev1,fprev2,lambdaprev1,lambdaprev2
real(8) cubic,lambda
real(8) constant,BS(2,2),C(2,1),D(2,1)   ! 'A' is changed to 'constant','B' is changed to 'BS'
real(8) a,b
cubic=0.d0
lambda=0.d0
constant=1.d0/(lambdaprev1-lambdaprev2)
BS(1,1)= 1.d0/lambdaprev1**2.d0
BS(1,2)=-1.d0/lambdaprev2**2.d0
BS(2,1)=-lambdaprev2/lambdaprev1**2.d0
BS(2,2)= lambdaprev1/lambdaprev2**2.d0
C(1,1)=fprev1-f0-df0*lambdaprev1
C(2,1)=fprev2-f0-df0*lambdaprev2
D=constant*matmul(BS,C)
a=D(1,1)
b=D(2,1)
lambda=(-b+dsqrt(b**2.d0-3.d0*a*df0))/(3.d0*a)
cubic=lambda
end function
!#################################################################################
function machineps()
!<this function is equivalent to epsilon(0.d0) epsilon being a built-in function
implicit none
real(8) machineps
machineps=1.d0
do 
        machineps=machineps/2.d0
        if ((1.d0+machineps).EQ.1.d0) exit
end do
machineps=machineps*2.d0
end function machineps
!#################################################################################
function quadratic(df0,f0,f1) !<quadratic backtracking
implicit none
real(8) df0,f0,f1
real(8) quadratic,lambda
quadratic=0.d0
lambda=0.d0
lambda=-0.5d0*df0/(f1-f0-df0)
quadratic=lambda
end function
!#################################################################################
function ser(normOld,normNew,oldDt,serpower)
implicit none

real(8) ser
real(8) serTest
real(8) serpower
real(8) normOld,normNew,oldDt



ser=oldDt

!if ((100.d0*dabs(normOld-normNew)/normOld).gt.2.d0) return


!print*,'old=',normOld
!print*,'new=',normnew
!print*,'old dt=',oldDt


serTest=((normOld/normNew)**serpower)*oldDt

!print*,'sertest=',sertest


serTest=min(serTest,2.d0*oldDt)
serTest=max(serTest,0.1d0*OldDt)





ser=serTest

!print*,'ser=',ser
!pause


end function ser













end module myFunctions

---