Merge branch 'feature-mhd' into fix_intel

examples/MHD/input_mhdchan.i3d

@@ -121,8 +121,8 @@ nclzSn = 0
 &MHDParam
 !====================
 
-mhd_equation = 'potential'
-hartmann = 6.d0
+mhd_equation = 'potential' ! MHD type of equation (potential or induction)
+hartmann = 6.d0            ! Hartmann number (Ha) ratio of electromagnetic force to the viscous force (BLsqrt(sigma/nu)) 
 
 /End
 

examples/MHD/input_mhdchan_test.i3d

@@ -121,8 +121,8 @@ nclzSn = 0
 &MHDParam
 !====================
 
-mhd_equation = 'potential'
-hartmann = 6.d0
+mhd_equation = 'potential' ! MHD type of equation (potential or induction)
+hartmann = 6.d0            ! Hartmann number (Ha) ratio of electromagnetic force to the viscous force (BLsqrt(sigma/nu)) 
 
 /End
 

examples/MHD/input_otv.i3d

@@ -8,8 +8,8 @@
 itype = 2
 
 ! Domain decomposition
-p_row=2                ! Row partition
-p_col=2                ! Column partition
+p_row=0                ! Row partition
+p_col=0                ! Column partition
 
 ! Mesh 
 nx=128                 ! X-direction nodes
@@ -97,9 +97,9 @@ nstat = 1             ! Size arrays for statistic collection
 &MHDParam
 !====================
 
-mhd_equation = 'induction'
-Stuart = 17.54386 
-Rem    = 17.54386
+mhd_equation = 'induction' ! MHD type of equation (potential or induction)
+Stuart = 17.54386          ! Stuart number ratio of electromagnetic to inertial forces (Ha^2/Rem)
+Rem    = 17.54386          ! Magnetic Reynolds number (UL/eta)
 
 /End
 

examples/MHD/input_otv_test.i3d

@@ -97,9 +97,9 @@ nstat = 1             ! Size arrays for statistic collection
 &MHDParam
 !====================
 
-mhd_equation = 'induction'
-Stuart = 50.d0
-Rem    = 50.d0
+mhd_equation = 'induction' ! MHD type of equation (potential or induction)
+Stuart = 50.0              ! Stuart number ratio of electromagnetic to inertial forces (Ha^2/Rem)
+Rem    = 50.0              ! Magnetic Reynolds number (UL/eta)
 
 /End
 

src/Case-ABL.f90

@@ -119,21 +119,6 @@ subroutine init_abl(ux1,uy1,uz1,ep1,phi1)
           endif
         endif
       enddo
-
-      !! Add random noise
-      !do j=1,xsize(2)
-      !  if (istret==0) y=real(j + xstart(2)-1-1,mytype)*dy
-      !  if (istret/=0) y=yp(j+xstart(2)-1)
-      !  !if (y.lt.50) then 
-      !  !  do k=1,xsize(3)
-      !  !  do i=1,xsize(1)
-      !  !    call random_number(phinoise)
-      !  !    phinoise=0.1*(phinoise*2.-1.)
-      !  !    phi1(i,j,k,1)=phi1(i,j,k,1)+phinoise
-      !  !  enddo
-      !  !  enddo
-      !  !endif
-      !enddo
     endif
 
     return

src/Case-Channel.f90

@@ -29,7 +29,7 @@ subroutine init_channel (ux1,uy1,uz1,ep1,phi1)
     use variables
     use param
     use MPI
-    use mhd, only : mhd_active, mhd_equation,Bm,Bmean
+    use mhd, only : mhd_equation,Bm,Bmean
 
     implicit none
 
@@ -93,9 +93,9 @@ subroutine init_channel (ux1,uy1,uz1,ep1,phi1)
                 if (idir_stream == 1) then
                    ux1(i,j,k)=one-y*y
                    uy1(i,j,k)=zero
-                   uz1(i,j,k)=zero !sin(real(i-1,mytype)*dx)+cos(real(k-1,mytype)*dz)
+                   uz1(i,j,k)=sin(real(i-1,mytype)*dx)+cos(real(k+xstart(3)-2,mytype)*dz)
                 else
-                        print *,'test'
+                   ! TODO: check if ux1 needs sin and cos 
                    uz1(i,j,k)=one-y*y
                    uy1(i,j,k)=zero
                    ux1(i,j,k)=zero
@@ -155,7 +155,7 @@ subroutine boundary_conditions_channel (ux,uy,uz,phi)
     use param
     use var, only : di2
     use variables
-    use mhd, only : Bm, mhd_active, mhd_equation 
+    use mhd, only : Bm, mhd_equation 
 
     implicit none
 
@@ -276,100 +276,13 @@ subroutine postprocess_channel(ux1,uy1,uz1,pp3,phi1,ep1)
     real(mytype), intent(in), dimension(xsize(1),xsize(2),xsize(3),numscalar) :: phi1
     real(mytype), intent(in), dimension(ph1%zst(1):ph1%zen(1),ph1%zst(2):ph1%zen(2),nzmsize,npress) :: pp3
 
-    real(mytype) :: ufriction,umax,resdux1,resduy1,resduz1,temp1,temp2
-    integer :: nxc, nyc, nzc, xsize1, xsize2, xsize3
-    integer :: i,j,k
-    integer :: code
-
-    if (nclx1==1 .and. xend(1)==nx) then
-       xsize1=xsize(1)-1
-    else
-       xsize1=xsize(1)
-    endif
-    if (ncly1==1 .and. xend(2)==ny) then
-       xsize2=xsize(2)-1
-    else
-       xsize2=xsize(2)
-    endif
-    if (nclz1==1 .and. xend(3)==nz) then
-       xsize3=xsize(3)-1
-    else
-       xsize3=xsize(3)
-    endif
-    if (nclx1==1) then
-       nxc=nxm
-    else
-       nxc=nx
-    endif
-    if (ncly1==1) then
-       nyc=nym
-    else
-       nyc=ny
-    endif
-    if (nclz1==1) then
-       nzc=nzm
-    else
-       nzc=nz
-    endif
-
-    !we only collect statistics every 10 time steps to save computational time
-    if (mod(itime, 10) == 0) then
-
-       temp1=zero
-
-       if (ncly1==2 .and. xstart(2)==1) then
-
-         ! bottom wall, one node away from wall
-         j=2
-         do k=1,xsize(3)
-             do i=1,xsize(1)
-              temp1=temp1+ux1(i,j,k)*ppy(j-1)
-             enddo
-         enddo
-
-       endif
-       !
-       if (nclyn==2 .and. xend(2)==ny) then
-
-         ! upper wall, one node away from wall
-         j=xsize(2)-1
-         do k=1,xsize(3)
-             do i=1,xsize(1)
-              temp1=temp1-ux1(i,j,k)*ppy(j+1)
-             enddo
-         enddo
-
-       endif
-       !
-       temp2=zero
-       do k=1,xsize3
-          do j=1,xsize2
-             do i=1,xsize1
-              temp2=max(temp2,ux1(i,j,k))
-             enddo
-          enddo
-       enddo
-
-       call MPI_ALLREDUCE(temp1,ufriction,1,real_type,MPI_SUM,MPI_COMM_WORLD,code)
-       call MPI_ALLREDUCE(temp2,     umax,1,real_type,MPI_MAX,MPI_COMM_WORLD,code)
-
-       ufriction=xnu*1.5_mytype*ufriction/(nxc*nzc) ! ub=2/3
-
-       if (nrank==0) then
-          write(52,'(3(E20.12))') (itime-1)*dt,ufriction,umax
-          flush(52)
-       endif
-
-    endif
-
-
   end subroutine postprocess_channel
 
   subroutine visu_channel_init(visu_initialised)
 
     use decomp_2d_io, only : decomp_2d_register_variable
     use visu, only : io_name, output2D
-    use mhd, only : mhd_active
+    use param, only : mhd_active
 
     implicit none
 
@@ -403,7 +316,8 @@ subroutine visu_channel(ux1, uy1, uz1, pp3, phi1, ep1, num)
     use var, only : ta2,tb2,tc2,td2,te2,tf2,di2,ta3,tb3,tc3,td3,te3,tf3,di3
     use var, ONLY : nzmsize
     use visu, only : write_field
-    use mhd, only : mhd_active,mhd_equation,Je, Bm
+    use param, only : mhd_active
+    use mhd, only : mhd_equation,Je, Bm
 
     use ibm_param, only : ubcx,ubcy,ubcz
 

src/Case-TGV.f90

@@ -28,7 +28,7 @@ subroutine init_tgv (ux1,uy1,uz1,ep1,phi1)
     use variables
     use param
     use MPI
-    use mhd,  only: mhd_active,Bm,Bmean
+    use mhd,  only: Bm,Bmean
     use visu, only: write_snapshot, end_snapshot
 
     implicit none
@@ -196,7 +196,7 @@ subroutine postprocess_tgv(ux1,uy1,uz1,phi1,ep1)
     USE var, only : ta1,tb1,tc1,td1,te1,tf1,tg1,th1,ti1,di1
     USE var, only : ta2,tb2,tc2,td2,te2,tf2,di2,ta3,tb3,tc3,td3,te3,tf3,di3
     USE ibm_param
-    use mhd, only : mhd_active,Bm,Je,Rem
+    use mhd, only : Bm,Je,Rem
 
     real(mytype),intent(in),dimension(xsize(1),xsize(2),xsize(3)) :: ux1, uy1, uz1,ep1
     real(mytype),intent(in),dimension(xsize(1),xsize(2),xsize(3),numscalar) :: phi1
@@ -485,7 +485,6 @@ subroutine visu_tgv_init (visu_initialised)
 
     use decomp_2d_io, only : decomp_2d_register_variable
     use visu, only : io_name, output2D
-    use mhd, only : mhd_active
 
     implicit none
 
@@ -521,7 +520,7 @@ subroutine visu_tgv(ux1, uy1, uz1, num)
     USE var, only : ta2,tb2,tc2,td2,te2,tf2,di2,ta3,tb3,tc3,td3,te3,tf3,di3
     use var, ONLY : nxmsize, nymsize, nzmsize
     use visu, only : write_field
-    use mhd, only : mhd_active,Bm,Je,del_cross_prod,rem
+    use mhd, only : Bm,Je,del_cross_prod,rem
     use ibm_param, only : ubcx,ubcy,ubcz
 
     implicit none

src/case.f90

@@ -493,7 +493,7 @@ end subroutine visu_case
   !##################################################################
   subroutine momentum_forcing(dux1, duy1, duz1, rho1, ux1, uy1, uz1, phi1)
 
-    use mhd, only: mhd_active,momentum_forcing_mhd
+    use mhd, only: momentum_forcing_mhd
 
     implicit none
 

src/mhd.f90

@@ -11,11 +11,9 @@ module mhd
   !
   implicit none
   !
-  logical :: mhd_active
   character(len=9) :: mhd_equation
   real(8) :: hartmann,stuart,rem
   !+------------+------------------------------------------------------+
-  !|  mhd_active| the swith to activate the mhd module.                |
   !|    hartmann| hartmann number, the ratio of lorentz force to       |
   !|            | viscous force                                        |
   !|      stuart| Stuart number, magnetic interaction parameter, ratio |
@@ -48,7 +46,9 @@ subroutine mhd_init
     !
     use param, only: re,ntime,itype,itype_channel,itype_tgv,zlz,dz
     use param, only: zero, one
+    use decomp_2d_mpi, only : decomp_2d_abort
     !
+    ! TODO fix these if
     ! stuart=hartmann**2/re
     if(stuart<=1.d-15) then
       stuart=hartmann**2/re
@@ -78,27 +78,24 @@ subroutine mhd_init
     !
     allocate(Bmean(xsize(1),xsize(2),xsize(3),1:3))
     !
-    if(nrank==0) print*,'** MHD fields allocated'
-    !
     if(mhd_equation == 'induction') then
-
+      ! TODO Check what is going to happen for other flow types  
       if(itype.eq.itype_tgv) then
         Bmean(:,:,:,1)=zero
         Bmean(:,:,:,2)=zero
         Bmean(:,:,:,3)=zero
       endif
     elseif(mhd_equation == 'potential') then
+      ! TODO Check what is going to happen for other flow types  
       if(itype.eq.itype_channel) then
         Bmean=zero
         !
         Bm(:,:,:,1)=zero
         Bm(:,:,:,2)=one
         Bm(:,:,:,3)=zero
       endif
-      !
-      if(nrank==0) print*,'** MHD fields initlised'
     else
-      stop ' mhd_equation not induction or potential'
+      call decomp_2d_abort(1, "Error in setup mhd equation induction or potential not selected")
     endif
     !
   end subroutine mhd_init
@@ -183,6 +180,10 @@ end function vortcal
   !+-------------------------------------------------------------------+
   subroutine momentum_forcing_mhd(dux1,duy1,duz1,ux1,uy1,uz1)
     !
+    USE decomp_2d_mpi, only : decomp_2d_abort
+
+    implicit none
+
     ! arguments
     real(mytype),intent(in),dimension(xsize(1),xsize(2),xsize(3)) ::   &
                                                 ux1,uy1,uz1
@@ -198,7 +199,7 @@ subroutine momentum_forcing_mhd(dux1,duy1,duz1,ux1,uy1,uz1)
     elseif(mhd_equation == 'potential') then
       Je=solve_mhd_potential_poisson(ux1,uy1,uz1)
     else
-      stop ' mhd_equation error '
+      call decomp_2d_abort(1, "Error in setup mhd equation induction or potential not selected")
     endif
     !
     do k = 1, xsize(3)
@@ -431,7 +432,8 @@ subroutine mhd_rhs_eq(dB,B,B0,ux1,uy1,uz1)
     integer :: i,j,k,is
     !
 
-    ! local 
+    ! local
+    ! TODO Need to check is all these local variable are necessary -> potential issue with INTEL  
     real(mytype) :: rrem
     real(mytype), save, allocatable, dimension(:,:,:) :: tx1,ty1,tz1,tx2,ty2,tz2, &
                                                          tx3,ty3,tz3,bx2,by2,bz2, &
@@ -696,6 +698,8 @@ subroutine mhd_rhs_eq(dB,B,B0,ux1,uy1,uz1)
 
   end subroutine mhd_rhs_eq
   !
+  ! TODO Check if the normal divergence can be used  
+  ! TODO Check if already allocated arrays can be re-used
   subroutine solve_poisson_mhd
     !
     use decomp_2d, only : zsize, ph1
@@ -749,6 +753,7 @@ function divergence_scalar(vec,nlock) result(pp3)
 
     nvect3=(ph1%zen(1)-ph1%zst(1)+1)*(ph1%zen(2)-ph1%zst(2)+1)*nzmsize
 
+    ! TODO If we properly use tmp var this should not be used 
     ta1(:,:,:) = vec(:,:,:,1)
     tb1(:,:,:) = vec(:,:,:,2)
     tc1(:,:,:) = vec(:,:,:,3)

src/module_param.f90

@@ -327,6 +327,7 @@ module param
   real(mytype) :: C_filter
   character(len=512) :: inflowpath
   logical :: validation_restart
+  logical :: mhd_active
 
   ! Logical, true when synchronization is needed
   logical, save :: sync_vel_needed = .true.