main.f90 4.53 KB
``````!
! -------------------------------------------------------------------
!                              main.f90
! This code do:
!       1) generate à initial function psi
!       2) Find the solution of psi with dirichlet condition
!       3) Find the solution of V with neuman condition
!
!       In the code all equation refer of psi, RHS as
!                   Laplacien(psi) = rhs
!
!       zeta, psi, v are the field to be compute
!       stage 1,2,3 refer as computing psi,zeta,v
!
!       psi = cos(2*pi*x) * sin(2*pi*y)
!           with periodic boundary at xa and xb and
!           diriclet bndy at yc,yd -> psi = 0
!
!       V =   cos(2*pi*x) *  cos(pi*y)
!           with periodic boundary at xa and xb and
!           neuman bndy at yc,yd -> dV/dx  = 0
!           V(yc) = -tau/po and V(yd) = 0
!
!
!          * * * * * * * * * * * *  y(nx)=yd = 0
!          *       kbdy=4        *
!          *                     *
!          *                     *
!          *                     *
!          * kbdy=1       kbdy=2 *
!          *                     *
!          *                     *
!          *                     *
!          *       kbdy=3        *
!          * * * * * * * * * * * *  y(1)=yc = Ly
!
!          x(0)=xa= 0            x(nx)=xb = Lx
!

!
! -------------------------------------------------------------------

module gridInfo
! Declare constante
real, parameter :: Pi = 3.1415927, nu = 0.016, f = 1.E-4, tau= 1E-3, Vmean = 0.0031

! LOW RES TEST
integer, parameter :: nx=97, ny=49
real, parameter :: Lx=1075, Ly=195, dx = (Lx-0) / float(nx-1), dy = (Ly-0)/float(ny-1)

! HIGH RES 10KM
!integer, parameter :: nx=3585, ny=65
!real, parameter :: Lx=10755, Ly=195, dx = (Lx-0) / float(nx-1), dy = (Ly-0)/float(ny-1)

end module gridInfo

program main
use gridInfo
implicit none

! Time parameter
integer niter, it
parameter (niter =3)

! Field
real zeta(nx+1,ny), psi(nx+1,ny), v(nx+1,ny), rhs(nx,ny), laplace(nx+1,ny), jacoZ(nx+1,ny), jacoV(nx+1,ny), realV(nx+1,ny)
real realU(nx+1,ny)
real A, epsi, sumErr, psi0(nx+1,ny)

! MudPack Multigrid parameter
integer iixp,jjyq,iiex,jjey,nnx,nny,isx,jsy,llwork
parameter (iixp = 3 , jjyq = 3 , iiex = 6, jjey = 5)         ! LOW RES
!parameter (iixp = 112 , jjyq = 2 , iiex = 6, jjey = 6)      ! HIGH RES D parameter (nnx=iixp*2**(iiex-1)+1, nny=jjyq*2**(jjey-1)+1)

include 'mudPackVar.f90'      ! Variable needed for mudPack package
include 'initialize.f90'      ! Initialising field and discretization of mud2
!include 'rhsRealV.f90'        ! Create a V field as calcul in MIT

write(*,*) '------ Starting Main.f90 ------'

! Set V ini as realV
v = realV

epsi=0.0    ! Time dependante contante Y = A + epsi(it) *  B

!--------------- ITERATIVE SOLVER --------------------
do it=1,niter

!if (it.gt.1) then         ! Set Guessing matrice to TRUE
iguessD   = 1
iguessN   = 1
!endif

!if (it>5) then            ! Set epsilon as function of time
!epsi = float((it-5))/float(niter-5)
!end if

psi0=psi

print*, 'Itertation:' , it , epsi

! Executing solver for a diri bndy
print*, 'zeta'
include 'rhsZeta.f90'
if (it.eq.1) then
call writeField('zetaini.bin',zeta)
endif
! set initialisation and iguess to 0
call mud2(iprmD,fprm,workD,cof,bndc,rhs,zeta(1:nx,:),mgopt,ierror)
call checkIerror(ierror, 'zeta')

! Executing sovler for a dirichlet bndy
print*, 'phi'
include 'rhsPsi.f90'
call mud2(iprmD,fprm,workD,cof,bndc,rhs,psi(1:nx,:),mgopt,ierror)
call checkIerror(ierror, 'psi')

! Executing solver for a neumman bndy
!print*, 'v'
!include 'rhsV.f90'
!call mud2(iprmN,fprm,workN,cof,bndc,rhs,v(1:nx,:),mgopt,ierror)
!call checkIerror(ierror, 'v')
! If we wante to force V
!print*, 'V sumdiff',sum(realV-v)
v = realV

print*, 'Sum Err',sum(psi-psi0)
! print*, 'Mean V', sum(v)/(nx*ny), 'Mean Jacov', sum(jacoV)/(nx*ny)
!print*, 'Mean Z', sum(v)/(nx*ny), 'Mean Jacoz', sum(jacoZ)/(nx*ny)

end do

! Wrinting field
call writeField("psi.bin", psi)
call writeField("zeta.bin",zeta)
call writeField("v.bin",v)
!call writeField("jacoV.bin",jacoV)
!call writeField("jacoZ.bin",jacoZ)
!call writeField("laplace.bin",laplace)

!call writeField("v.bin",v)

write(*,*) '---- End of Main.f90 ----'
end program main

include 'subMud.f90'

``````