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
!


! see tmud2.m from the mudpack5.0.1 to have more information
!
! -------------------------------------------------------------------  

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 readV(nx,ny)
  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'