TEHD Implementation in MagIC

Add the dielectrophoretic force in MagIC using a
temperature invarient electric field. The mode
currently availavble is  mode 1 with enables
convection only. The electric Rayleigh number RaE
(rae) controls the amount of electric buoyancy.

Author: ygaillard

src/Namelists.f90

@@ -77,7 +77,7 @@ subroutine readNamelists(tscheme)
       &    mpi_transp,l_adv_curl,mpi_packing
 
       namelist/phys_param/                                        &
-      &    ra,raxi,pr,sc,prmag,ek,epsc0,epscxi0,radratio,Bn,      &
+      &    ra,rae,raxi,pr,sc,prmag,ek,epsc0,epscxi0,radratio,Bn,      &
       &    ktops,kbots,ktopv,kbotv,ktopb,kbotb,kbotxi,ktopxi,     &
       &    s_top,s_bot,impS,sCMB,xi_top,xi_bot,impXi,xiCMB,       &
       &    nVarCond,con_DecRate,con_RadRatio,con_LambdaMatch,     &
@@ -323,6 +323,7 @@ subroutine readNamelists(tscheme)
       l_AB1    =.false.
       l_bridge_step=.true.
       l_onset  =.false.
+      l_ehd_dep=.true.
 
       if ( mode == 1 ) then
          !-- Only convection:
@@ -427,7 +428,7 @@ subroutine readNamelists(tscheme)
          l_chemical_conv = .false.
       end if
 
-      if ( ra == 0.0_cp ) l_heat=.false.
+      if ( ra == 0.0_cp .and. rae == 0.0_cp ) l_heat=.false.
 
       if ( ek < 0.0_cp ) l_non_rot= .true.
       if ( l_non_rot ) then
@@ -928,6 +929,7 @@ subroutine writeNamelists(n_out)
 
       write(n_out,*) "&phys_param"
       write(n_out,'(''  ra              ='',ES14.6,'','')') ra
+      write(n_out,'(''  rae             ='',ES14.6,'','')') rae
       write(n_out,'(''  raxi            ='',ES14.6,'','')') raxi
       write(n_out,'(''  pr              ='',ES14.6,'','')') pr
       write(n_out,'(''  sc              ='',ES14.6,'','')') sc

src/get_nl.f90

@@ -22,15 +22,15 @@ module grid_space_arrays_mod
    use truncation, only: n_phi_max, nlat_padded
    use radial_functions, only: or2, orho1, beta, otemp1, visc, r, or3, &
        &                       lambda, or4, or1
-   use physical_parameters, only: LFfac, n_r_LCR, prec_angle, ViscHeatFac,    &
-        &                         oek, po, dilution_fac, ra, opr, OhmLossFac, &
+   use physical_parameters, only: radratio, LFfac, n_r_LCR, prec_angle, ViscHeatFac,    &
+        &                         oek, po, dilution_fac, ra, rae, opr, OhmLossFac, &
         &                         epsPhase, phaseDiffFac, penaltyFac, tmelt
    use horizontal_data, only: sinTheta, cosTheta, phi, O_sin_theta_E2, &
        &                      cosn_theta_E2, O_sin_theta
    use parallel_mod, only: get_openmp_blocks
    use constants, only: two, third, one
    use logic, only: l_conv_nl, l_heat_nl, l_mag_nl, l_anel, l_mag_LF, l_adv_curl, &
-       &            l_chemical_conv, l_precession, l_centrifuge, l_phase_field
+       &            l_chemical_conv, l_precession, l_centrifuge, l_phase_field, l_ehd_dep
 
    implicit none
 
@@ -46,6 +46,7 @@ module grid_space_arrays_mod
       real(cp), allocatable :: VSr(:,:), VSt(:,:), VSp(:,:)
       real(cp), allocatable :: VXir(:,:), VXit(:,:), VXip(:,:)
       real(cp), allocatable :: heatTerms(:,:), phiTerms(:,:)
+      real(cp), allocatable :: DEPFr(:,:)
 
       !----- Fields calculated from these help arrays by legtf:
       real(cp), allocatable :: vrc(:,:), vtc(:,:), vpc(:,:)
@@ -182,6 +183,12 @@ subroutine initialize(this)
          bytes_allocated=bytes_allocated+2*n_phi_max*nlat_padded*SIZEOF_DEF_REAL
       end if
 
+      if ( l_ehd_dep ) then
+         allocate( this%DEPFr(nlat_padded,n_phi_max) )
+         this%DEPFr(:,:)=0.0_cp
+         bytes_allocated=bytes_allocated + 1*n_phi_max*nlat_padded*SIZEOF_DEF_REAL
+         end if
+         
    end subroutine initialize
 !----------------------------------------------------------------------------
    subroutine finalize(this)
@@ -195,6 +202,7 @@ subroutine finalize(this)
       deallocate( this%VxBr, this%VxBt, this%VxBp, this%VSr, this%VSt, this%VSp )
       if ( l_chemical_conv ) deallocate( this%VXir, this%VXit, this%VXip )
       if ( l_precession ) deallocate( this%PCr, this%PCt, this%PCp )
+      if ( l_ehd_dep ) deallocate( this%DEPFr )
       if ( l_centrifuge ) deallocate( this%CAr, this%CAt )
       if ( l_adv_curl ) deallocate( this%cvtc, this%cvpc )
       if ( l_phase_field ) deallocate( this%phic, this%phiTerms )
@@ -257,6 +265,12 @@ subroutine get_nl(this, time, nR, nBc, lRmsCalc)
             &        this%cbtc(:,nPhi)*this%brc(:,nPhi) )
          end if      ! Lorentz force required ?
 
+         if ( l_ehd_dep .and. (nBc == 0 .or. lRmsCalc) .and. nR>n_r_LCR ) then
+            !------ Get the dielectrophoretic force:
+            !---- RaE * eta**2 / (1-eta)**4 * Sc / r**5
+            this%DEPFr(:,nPhi)=  rae * opr *  radratio**2/(1.0D0-radratio)**4 * this%sc(:,nPhi) * or3(nR)
+         end if      ! DEP force required ?
+
          if ( l_conv_nl .and. (nBc == 0 .or. lRmsCalc) ) then
 
             if ( l_adv_curl ) then ! Advection is \curl{u} \times u

src/logic.f90

@@ -72,6 +72,7 @@ module logic
    logical :: l_perpPar      ! Switch for calculation of of kinetic energy perpendicular+parallel to the rotation axis
    logical :: l_LCR          ! Switch for zero electrical conductivity beyond r_LCR
    logical :: lVerbose       ! Switch for detailed information about run progress
+   logical :: l_ehd_dep      ! Switch for dilectrophoretic force
 
    logical :: l_PressGraph   ! Store pressure in graphic files
 

src/phys_param.f90

@@ -37,6 +37,7 @@ module physical_parameters
    !-- Dimensionless parameters:
    real(cp) :: radratio       ! aspect ratio
    real(cp) :: ra             ! Rayleigh number
+   real(cp) :: rae            ! Electric Rayleigh number
    real(cp) :: raxi           ! Chemical composition-based Rayleigh number
    real(cp) :: sc             ! Schmidt number (i.e. chemical Prandtl number)
    real(cp) :: ek             ! Ekman number

src/rIter.f90

@@ -21,7 +21,7 @@ module rIter_mod
        &            l_precession, l_centrifuge, l_adv_curl,          &
        &            l_double_curl, l_parallel_solve, l_single_matrix,&
        &            l_temperature_diff, l_RMS, l_phase_field,        &
-       &            l_onset, l_DTrMagSpec
+       &            l_onset, l_DTrMagSpec, l_ehd_dep
    use radial_data, only: n_r_cmb, n_r_icb, nRstart, nRstop, nRstartMag, &
        &                  nRstopMag
    use radial_functions, only: or2, orho1, l_R
@@ -667,6 +667,10 @@ subroutine transform_to_lm_space(this, nR, lRmsCalc, dVSrLM, dVXirLM, dphidt)
                this%gsa%Advr(:, nPhi)=this%gsa%Advr(:,nPhi) + this%gsa%CAr(:,nPhi)
                this%gsa%Advt(:, nPhi)=this%gsa%Advt(:,nPhi) + this%gsa%CAt(:,nPhi)
             end if
+
+            if ( l_ehd_dep ) then
+               this%gsa%Advr(:, nPhi)=this%gsa%Advr(:,nPhi) + this%gsa%DEPFr(:,nPhi)
+            end if   
          end do
          !$omp end parallel