Version seems to be working but untested

gafusion · Sep 12, 2024 · c2d5e58 · c2d5e58
1 parent 8ab540a
commit c2d5e58
Show file tree

Hide file tree

Showing 6 changed files with 32 additions and 60 deletions.
diff --git a/cgyro/src/cgyro_init_arrays.F90 b/cgyro/src/cgyro_init_arrays.F90
@@ -462,7 +462,7 @@ subroutine cgyro_init_arrays
            sm = sdlnndr(is)+sdlntdr(is)*(energy(ie)-1.5)
 
            ! generalized beta/drift shear (acts on H)
-           sb = beta_star(0)*sbeta_star(is)*energy(ie)/bmag(it)**2
+           sb = sbeta_star(is)*energy(ie)/bmag(it)**2
 
            arg = -k_theta_base*itor*length/(2*pi)
 

diff --git a/cgyro/src/cgyro_make_profiles.F90 b/cgyro/src/cgyro_make_profiles.F90
@@ -122,13 +122,14 @@ subroutine cgyro_make_profiles
      z_eff   = z_eff_loc
      b_unit  = b_unit_loc
 
-     dens(1:n_species) = dens_loc(1:n_species)     
-     temp(1:n_species) = temp_loc(1:n_species)     
-     dlnndr(1:n_species) = dlnndr_loc(1:n_species)     
-     dlntdr(1:n_species) = dlntdr_loc(1:n_species)
+     dens(1:n_species)    = dens_loc(1:n_species)     
+     temp(1:n_species)    = temp_loc(1:n_species)     
+     dlnndr(1:n_species)  = dlnndr_loc(1:n_species)     
+     dlntdr(1:n_species)  = dlntdr_loc(1:n_species)
      sdlnndr(1:n_species) = sdlnndr_loc(1:n_species)     
      sdlntdr(1:n_species) = sdlntdr_loc(1:n_species)     
-
+     sbeta_star(1:n_species) = sbeta_loc(1:n_species)     
+
      if (ae_flag == 1) then
         is_ele = n_species+1
      else
@@ -252,13 +253,13 @@ subroutine cgyro_make_profiles
         dlntdr_ele = dlntdr_ae
      else
         is_ele = -1
-        do is=1, n_species
-           if(z(is) < 0.0) then
+        do is=1,n_species
+           if (z(is) < 0.0) then
               is_ele = is
               exit
            endif
         enddo
-        if(is_ele == -1) then
+        if (is_ele == -1) then
            call cgyro_error('No electron species specified')
            return
         endif
@@ -284,8 +285,8 @@ subroutine cgyro_make_profiles
      ! Always compute beta_* consistently with parameters in input.cgyro and then re-scale
      ! note: beta_star(0) will be over-written with sonic rotation
      call set_betastar
-     beta_star(0) = beta_star(0)*beta_star_scale
-     beta_star_fac  = beta_star_fac * beta_star_scale
+     beta_star(0)  = beta_star(0)*beta_star_scale
+     beta_star_fac = beta_star_fac*beta_star_scale
 
   endif
 

diff --git a/cgyro/src/cgyro_write_initdata.f90 b/cgyro/src/cgyro_write_initdata.f90
@@ -132,12 +132,12 @@ subroutine cgyro_write_initdata
      ! where x = r/rho. The half-domain has -L/4 < r < L/4.
      if (profile_shear_flag == 1) then
         write(io,*)
-        write(io,'(a)') ' i  s(a/Ln)  (a/Ln)_L  (a/Ln)_R  |  s(a/Lt)  (a/Lt)_L  (a/Lt)_R ' 
+        write(io,'(a)') ' i  s(a/Ln)  (a/Ln)_L  (a/Ln)_R  |  s(a/Lt)  (a/Lt)_L  (a/Lt)_R  |  sbeta_*' 
         do is=1,n_species
            dn = sdlnndr(is)*length/rho/4
            dt = sdlntdr(is)*length/rho/4
-           write(io,'(t1,i2,3(1x,1pe9.2),2x,3(1x,1pe9.2))') &
-             is,sdlnndr(is),dlnndr(is)-dn,dlnndr(is)+dn,sdlntdr(is),dlntdr(is)-dt,dlntdr(is)+dt
+           write(io,'(t1,i2,3(1x,1pe9.2),2x,3(1x,1pe9.2),2x,1(1x,1pe9.2))') &
+             is,sdlnndr(is),dlnndr(is)-dn,dlnndr(is)+dn,sdlntdr(is),dlntdr(is)-dt,dlntdr(is)+dt,sbeta_star(is)
         enddo
      endif
 

diff --git a/f2py/expro/expro.f90 b/f2py/expro/expro.f90
@@ -121,7 +121,6 @@ module expro
        expro_dlntedr,&
        expro_sdlnnedr,&
        expro_sdlntedr,&
-       expro_sbetae,&
        expro_dlnptotdr,&
        expro_w0p,&
        expro_surf,&
@@ -161,8 +160,7 @@ module expro
        expro_dlnnidr,&
        expro_dlntidr,&
        expro_sdlnnidr,&
-       expro_sdlntidr,&
-       expro_sbetai
+       expro_sdlntidr
   !-----------------------------------------------------------
 
   ! input.gacode.geo dimension and arrays
@@ -289,7 +287,6 @@ subroutine expro_init(flag)
        allocate(expro_dlntedr(nexp))      ; expro_dlntedr = 0.0
        allocate(expro_sdlnnedr(nexp))     ; expro_sdlnnedr = 0.0
        allocate(expro_sdlntedr(nexp))     ; expro_sdlntedr = 0.0
-       allocate(expro_sbetae(nexp))       ; expro_sbetae = 0.0
        allocate(expro_dlnptotdr(nexp))    ; expro_dlnptotdr = 0.0
        allocate(expro_w0p(nexp))          ; expro_w0p = 0.0
        allocate(expro_surf(nexp))         ; expro_surf = 0.0
@@ -331,7 +328,6 @@ subroutine expro_init(flag)
        allocate(expro_dlntidr(nion,nexp))  ; expro_dlntidr = 0.0
        allocate(expro_sdlnnidr(nion,nexp)) ; expro_sdlnnidr = 0.0
        allocate(expro_sdlntidr(nion,nexp)) ; expro_sdlntidr = 0.0
-       allocate(expro_sbetai(nion,nexp))   ; expro_sbetai = 0.0
 
     else
 
@@ -414,7 +410,6 @@ subroutine expro_init(flag)
        deallocate(expro_dlntedr)      
        deallocate(expro_sdlnnedr)      
        deallocate(expro_sdlntedr)      
-       deallocate(expro_sbetae)      
        deallocate(expro_dlnptotdr)    
        deallocate(expro_w0p)          
        deallocate(expro_surf)          
@@ -456,7 +451,6 @@ subroutine expro_init(flag)
        deallocate(expro_dlntidr)  
        deallocate(expro_sdlnnidr) 
        deallocate(expro_sdlntidr)  
-       deallocate(expro_sbetai)  
 
     endif
 

diff --git a/f2py/expro/expro_locsim.f90 b/f2py/expro/expro_locsim.f90
@@ -73,6 +73,7 @@ module expro_locsim_interface
   double precision :: beta_star_loc
 
   double precision, dimension(9) :: mass_loc
+
   double precision, dimension(9) :: z_loc
   double precision, dimension(9) :: dens_loc
   double precision, dimension(9) :: temp_loc
@@ -283,7 +284,6 @@ subroutine expro_locsim_profiles(&
   dens_exp(n_species_exp,:)    = expro_ne(:)
   dlnndr_exp(n_species_exp,:)  = expro_dlnnedr(:)*a_meters 
   sdlnndr_exp(n_species_exp,:) = expro_sdlnnedr(:)*a_meters
-  sbeta_exp(n_species_exp,:)   = expro_sbetae(:)*a_meters
 
   mass_loc(n_species_exp) = expro_masse
   z_loc(n_species_exp) = -1d0
@@ -309,8 +309,6 @@ subroutine expro_locsim_profiles(&
         sdlnndr_exp(i_ion,:) = expro_sdlnnidr(i_ion,:)*a_meters
      endif
 
-     sbeta_exp(i_ion,:) = expro_sbetai(i_ion,:)*a_meters 
-
   enddo
 
   ! Rotation
@@ -377,7 +375,6 @@ subroutine expro_locsim_profiles(&
      call cub_spline1(rmin_exp,dlnndr_exp(i,:),expro_n_exp,rmin,dlnndr_loc(i))
      call cub_spline1(rmin_exp,sdlntdr_exp(i,:),expro_n_exp,rmin,sdlntdr_loc(i))
      call cub_spline1(rmin_exp,sdlnndr_exp(i,:),expro_n_exp,rmin,sdlnndr_loc(i))
-     call cub_spline1(rmin_exp,sbeta_exp(i,:),expro_n_exp,rmin,sbeta_loc(i))
      beta_star_loc = beta_star_loc+dens_loc(i)*temp_loc(i)*(dlnndr_loc(i)+dlntdr_loc(i))
   enddo
 
@@ -392,6 +389,11 @@ subroutine expro_locsim_profiles(&
   betae_loc = 4.027e-3*dens_loc(n_species_exp)*temp_loc(n_species_exp)/b_unit_loc**2
   beta_star_loc = beta_star_loc*betae_loc/(dens_loc(n_species_exp)*temp_loc(n_species_exp))
 
+  ! Add extra effective curvature terms to omega_star shear
+  sdlnndr_loc = sdlnndr_loc+(1.5*dlntdr_loc**2+dlnndr_loc**2-dlnndr_loc*dlntdr_loc)*rhos_loc/a_meters
+  sdlntdr_loc = sdlntdr_loc+dlntdr_loc**2*rhos_loc/a_meters
+  sbeta_loc   = beta_star_loc*(sdlnndr_loc+sdlntdr_loc-2*dlntdr_loc*dlnndr_loc*rhos_loc/a_meters)/(dlnndr_loc+dlntdr_loc)
+
   if (numeq_flag == 1 .and. expro_nfourier > 0) then
 
      geo_ny_loc = expro_nfourier

diff --git a/f2py/expro/expro_util.f90 b/f2py/expro/expro_util.f90
@@ -22,8 +22,6 @@ subroutine expro_compute_derived
   double precision, dimension(:), allocatable :: temp
   double precision, dimension(:), allocatable :: cc
   double precision, dimension(:), allocatable :: loglam
-  double precision, dimension(:,:), allocatable :: sni,sti
-  double precision, dimension(:), allocatable  :: sne,ste
 
   double precision :: r_min
   double precision :: fa,fb
@@ -67,10 +65,6 @@ subroutine expro_compute_derived
   !
   allocate(torflux(expro_n_exp))
   allocate(temp(expro_n_exp))
-  allocate(sne(expro_n_exp))
-  allocate(ste(expro_n_exp))
-  allocate(sni(expro_n_ion,expro_n_exp))
-  allocate(sti(expro_n_ion,expro_n_exp))
 
   torflux(:) = expro_torfluxa*expro_rho(:)**2
 
@@ -145,17 +139,17 @@ subroutine expro_compute_derived
   ! NOTE: expro_sdln* will be renormalized after calculation of rhos later
 
   ! sne = -ne''/ne (1/m^2) [not fully normalized yet]
-  call bound_deriv(sne,expro_ne*expro_dlnnedr,expro_rmin,expro_n_exp)
-  sne = sne/expro_ne
+  call bound_deriv(expro_sdlnnedr,expro_ne*expro_dlnnedr,expro_rmin,expro_n_exp)
+  expro_sdlnnedr = expro_sdlnnedr/expro_ne
 
   ! sTe = -Te''/Te (1/m^2) [not fully normalized yet]
-  call bound_deriv(ste,expro_te*expro_dlntedr,expro_rmin,expro_n_exp)
-  ste = ste/expro_te
+  call bound_deriv(expro_sdlntedr,expro_te*expro_dlntedr,expro_rmin,expro_n_exp)
+  expro_sdlntedr = expro_sdlntedr/expro_te
 
   expro_dlnnidr = 0d0
   expro_dlntidr = 0d0
-  sni = 0d0
-  sti = 0d0
+  expro_sdlnnidr = 0d0
+  expro_sdlntidr = 0d0
 
   do is=1,expro_n_ion
      if (minval(expro_ni(is,:)) > 0d0 .and. minval(expro_ti(is,:)) > 0d0) then
@@ -166,37 +160,21 @@ subroutine expro_compute_derived
         call bound_deriv(expro_dlntidr(is,:),-log(expro_ti(is,:)),expro_rmin,expro_n_exp)
 
         ! sni = -ni''/ni (1/m^2) [not fully normalized yet]
-        call bound_deriv(sni(is,:),expro_ni(is,:)*expro_dlnnidr(is,:),expro_rmin,expro_n_exp)
-        sni(is,:) = sni(is,:)/expro_ni(is,:)
+        call bound_deriv(expro_sdlnnidr(is,:),expro_ni(is,:)*expro_dlnnidr(is,:),expro_rmin,expro_n_exp)
+        expro_sdlnnidr(is,:) = expro_sdlnnidr(is,:)/expro_ni(is,:)
 
         ! sTi = -Ti''/Ti (1/m^2) [not fully normalized yet]
         call bound_deriv(expro_sdlntidr(is,:),expro_ti(is,:)*expro_dlntidr(is,:),expro_rmin,expro_n_exp)
-        sti(is,:) = sti(is,:)/expro_ti(is,:)
+        expro_sdlntidr(is,:) = expro_sdlntidr(is,:)/expro_ti(is,:)
      endif
   enddo
 
-  ! omega_star shear
-  expro_sdlnnedr = sne+1.5*expro_dlntedr**2+expro_dlnnedr**2-expro_dlnnedr*expro_dlntedr
-  expro_sdlnnidr = sni-1.5*expro_dlntidr**2+expro_dlnnidr**2-expro_dlnnidr*expro_dlntidr
-
-  expro_sdlntedr = ste+expro_dlntedr**2
-  expro_sdlntidr = sti+expro_dlntidr**2
-
-  ! beta_star shear 
-  expro_sbetae = (sne+ste-2*expro_dlntedr*expro_dlnnedr)/(expro_dlnnedr+expro_dlntedr)
-  expro_sbetai = (sni+sti-2*expro_dlntidr*expro_dlnnidr)/(expro_dlnnidr+expro_dlntidr)
-
   ! 1/L_Ptot = -dln(Ptot)/dr (1/m)
   if (minval(expro_ptot) > 0d0) then
      call bound_deriv(expro_dlnptotdr,-log(expro_ptot),expro_rmin,expro_n_exp)
   else
      expro_dlnptotdr = 0d0
   endif
-
-  deallocate(sne)
-  deallocate(ste)
-  deallocate(sni)
-  deallocate(sti)
   !--------------------------------------------------------------------
 
   !-------------------------------------------------------------------
@@ -353,18 +331,15 @@ subroutine expro_compute_derived
   !-----------------------------------------------------------------
 
   !-----------------------------------------------------------------
-  ! Renormalize shearing parameters
+  ! Renormalize curvatures
   !
   ! sn = -n''/n*rhos (1/m) 
   ! sT = -T''/T*rhos (1/m)
-  ! sbeta = sbeta*rhos (1/m)
   expro_sdlnnedr = expro_sdlnnedr*expro_rhos(:)
   expro_sdlntedr = expro_sdlntedr*expro_rhos(:)
-  expro_sbetae = expro_sbetae*expro_rhos(:)
   do is=1,expro_n_ion
      expro_sdlnnidr(is,:) = expro_sdlnnidr(is,:)*expro_rhos(:)
      expro_sdlntidr(is,:) = expro_sdlntidr(is,:)*expro_rhos(:)
-     expro_sbetai(is,:) = expro_sbetai(is,:)*expro_rhos(:)
   enddo
   !-----------------------------------------------------------------