From c2e0b2a98faba6403ae9382f2186b966cc587479 Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Wed, 2 Oct 2024 08:22:55 -0700 Subject: [PATCH] New selection logic for global mode --- cgyro/bin/cgyro_json | 44 ++++++++++++++++++++++++--- cgyro/bin/cgyro_parse.py | 4 +-- cgyro/src/cgyro_advect_wavenumber.F90 | 37 +++++++++++----------- cgyro/src/cgyro_globals.F90 | 4 +-- cgyro/src/cgyro_globalshear.F90 | 16 +++++----- cgyro/src/cgyro_init_arrays.F90 | 2 +- cgyro/src/cgyro_make_profiles.F90 | 12 +++++--- cgyro/src/cgyro_read_input.f90 | 4 +-- cgyro/src/cgyro_rhs.F90 | 12 +++++--- cgyro/src/cgyro_write_initdata.f90 | 4 +-- 10 files changed, 89 insertions(+), 50 deletions(-) diff --git a/cgyro/bin/cgyro_json b/cgyro/bin/cgyro_json index 1a01be4a6..2dfe845fe 100755 --- a/cgyro/bin/cgyro_json +++ b/cgyro/bin/cgyro_json @@ -1,5 +1,6 @@ #!/usr/bin/env python +import sys import json import numpy as np @@ -12,13 +13,26 @@ try: version = file.read().rstrip() except: version = 'unavailable' - -#-------------------------------------------------------------- +if sim.n_species == 1: + print('Conversion not implemented for adiabatic electrons.') + sys.exit() + +# Electron index +for i in range(sim.n_species): + if sim.z[i] < 0.0: + ielec = i + +#================================================================== # json.cgyro.localdump d = {} +# Code-specific parameters +d['N_THETA'] = sim.n_theta +d['N_RADIAL'] = sim.n_theta + +# Physics parameters d['RMIN'] = sim.rmin d['RMAJ'] = sim.rmaj d['Q'] = sim.q @@ -52,8 +66,11 @@ d['SHAPE_COS'] = sim.shape_cos[:].tolist() # JSON output with open('json.cgyro.localdump','w') as f: json.dump(d,f,indent=2,sort_keys=True) +#================================================================== +#================================================================== # standard localdump output +# nmax = 0 with open('out.cgyro.localdump','w') as f: for key in d: @@ -70,17 +87,24 @@ with open('out.cgyro.localdump','w') as f: else: f.write(key+'='+str(x)+'\n') -#-------------------------------------------------------------- +#================================================================== + +#================================================================== # json.cgyro.imas nmax = nmax+1 +# Reference dimensions in CGYRO units lref=sim.rmaj vthref = np.sqrt(2.0) bref = sim.b_gs2 +# Initialize dictionary d = {} +#-------------------------------------------------------------------- +# Species and profiles + d['name'] = 'CGYRO' d['version'] = version d['r_minor_norm'] = sim.rmin/lref @@ -92,15 +116,23 @@ d['density_norm'] = sim.dens[:].tolist() d['temperature_norm'] = sim.temp[:].tolist() d['density_log_gradient_norm'] = (lref*sim.dlnndr[:]).tolist() d['temperature_log_gradient_norm'] = (lref*sim.dlntdr[:]).tolist() +#-------------------------------------------------------------------- -# Rotation parameters +#-------------------------------------------------------------------- +# Rotation, collisons, beta d['velocity_tor_norm'] = sim.mach/vthref # -L/v 1/B (r/q Bu) (-dw/dr) = L/v Bu/B (r/q dw/dr) = L/v Bu/B (-GAMMA_E) d['shearing_rate_norm'] = -(lref/vthref)/bref*sim.gamma_e # -L^2/v dw/dr = -L/v (L dw/dr) = -L/v (-GAMMA_P) d['velocity_tor_gradient_norm'] = lref/vthref*sim.gamma_p - +# nu_ee +d['collisionality_norm'] = (lref/vthref)*sim.nu[ielec] +# electron beta d['beta_reference'] = sim.betae_unit/bref**2 +#----------------------------------------------------------------------- + +#----------------------------------------------------------------------- +# Geometry d['dgeometric_axis_r_dr_minor'] = sim.shift d['dgeometric_axis_z_dr_minor'] = sim.dzmag @@ -119,7 +151,9 @@ d['dc_dr_minor_norm'] = (sim.shape_s_cos[:nmax]*lref/sim.rmin).tolist() d['ds_dr_minor_norm'] = (sim.shape_s_sin[:nmax]*lref/sim.rmin).tolist() d['ds_dr_minor_norm'][1] = sim.s_delta*lref/sim.rmin/np.cos(np.arcsin(sim.delta)) d['ds_dr_minor_norm'][2] = -sim.s_zeta*lref/sim.rmin +#----------------------------------------------------------------------- with open('json.cgyro.imas','w') as f: json.dump(d,f,indent=2,sort_keys=True) +#================================================================== diff --git a/cgyro/bin/cgyro_parse.py b/cgyro/bin/cgyro_parse.py index 221e9016d..d2f613d11 100644 --- a/cgyro/bin/cgyro_parse.py +++ b/cgyro/bin/cgyro_parse.py @@ -79,9 +79,9 @@ x.add('HIPREC_FLAG','0') x.add('UDSYMMETRY_FLAG','0') x.add('SHEAR_METHOD','2') +x.add('GLOBAL_FLAG','0') x.add('N_GLOBAL','4') x.add('NU_GLOBAL','15.0') -x.add('PROFILE_SHEAR_FLAG','0') x.add('THETA_PLOT','1') x.add('GPU_BIGMEM_FLAG','1') x.add('UPWIND_SINGLE_FLAG','0') @@ -139,7 +139,7 @@ x.add('DLNTDR','1.0',n=n) x.add('SDLNNDR','0.0',n=n) x.add('SDLNTDR','0.0',n=n) -x.add('SBETA_STAR','0.0',n=n) +x.add('SBETA','0.0',n=n) x.add('DLNNDR_SCALE','1.0',n=n) x.add('DLNTDR_SCALE','1.0',n=n) diff --git a/cgyro/src/cgyro_advect_wavenumber.F90 b/cgyro/src/cgyro_advect_wavenumber.F90 index de6b83019..a53920fbd 100644 --- a/cgyro/src/cgyro_advect_wavenumber.F90 +++ b/cgyro/src/cgyro_advect_wavenumber.F90 @@ -16,10 +16,9 @@ subroutine cgyro_advect_wavenumber(ij) integer :: ir,l,ll,j,iccj,ivc,itor,llnt complex :: rl,he1,he2 - if (nonlinear_flag == 0) return + if (nonlinear_flag == 0 .or. source_flag == 0) return - if (source_flag == 1) then - call timer_lib_in('shear') + call timer_lib_in('shear') #if defined(OMPGPU) !$omp target teams distribute parallel do simd collapse(4) & @@ -31,14 +30,14 @@ subroutine cgyro_advect_wavenumber(ij) #else !$omp parallel do collapse(4) private(ivc,ir,l,iccj,j,ll,rl,llnt,he1,he2) #endif - do itor=nt1,nt2 - do ivc=1,nv_loc - do ir=1,n_radial - do j=1,n_theta - iccj = (ir-1)*n_theta+j + do itor=nt1,nt2 + do ivc=1,nv_loc + do ir=1,n_radial + do j=1,n_theta + iccj = (ir-1)*n_theta+j - ! Wavenumber advection ExB shear - if (shear_method == 2) then + ! Wavenumber advection ExB shear + if (shear_method == 2) then rl = 0.0 #if (!defined(OMPGPU)) && defined(_OPENACC) !$acc loop seq @@ -63,10 +62,10 @@ subroutine cgyro_advect_wavenumber(ij) rl = rl+c_wave(l)*(he1-he2) enddo rhs(iccj,ivc,itor,ij) = rhs(iccj,ivc,itor,ij) + omega_eb_base*itor*rl - endif + endif - ! Wavenumber advection profile shear - if (profile_shear_flag == 1) then + ! Wavenumber advection profile shear + if (global_flag == 1) then iccj = (ir-1)*n_theta+j rl = rhs(iccj,ivc,itor,ij) #if (!defined(OMPGPU)) && defined(_OPENACC) @@ -91,14 +90,12 @@ subroutine cgyro_advect_wavenumber(ij) rl = rl-c_wave(l)*(he1-he2) enddo rhs(iccj,ivc,itor,ij) = rl - endif - enddo - enddo - enddo + endif + enddo + enddo enddo + enddo - call timer_lib_out('shear') - - endif + call timer_lib_out('shear') end subroutine cgyro_advect_wavenumber diff --git a/cgyro/src/cgyro_globals.F90 b/cgyro/src/cgyro_globals.F90 index 67087a032..c91c3bdf0 100644 --- a/cgyro/src/cgyro_globals.F90 +++ b/cgyro/src/cgyro_globals.F90 @@ -91,10 +91,10 @@ module cgyro_globals integer :: hiprec_flag integer :: udsymmetry_flag integer :: shear_method + integer :: global_flag integer :: n_global real :: nu_global integer :: psym_flag - integer :: profile_shear_flag integer :: theta_plot integer :: gpu_bigmem_flag integer :: upwind_single_flag @@ -138,7 +138,7 @@ module cgyro_globals real, dimension(11) :: dlntdr real, dimension(11) :: sdlnndr real, dimension(11) :: sdlntdr - real, dimension(11) :: sbeta_star + real, dimension(11) :: sbeta integer :: subroutine_flag ! only used for cgyro_read_input diff --git a/cgyro/src/cgyro_globalshear.F90 b/cgyro/src/cgyro_globalshear.F90 index 9b4c34874..4ace2a0a1 100644 --- a/cgyro/src/cgyro_globalshear.F90 +++ b/cgyro/src/cgyro_globalshear.F90 @@ -16,7 +16,7 @@ subroutine cgyro_globalshear(ij) integer :: ir,l,ll,j,iccj,ivc,itor,llnt complex :: rl,h1,h2 - if (nonlinear_flag == 0) return + if (nonlinear_flag == 0 .or. source_flag == 0) return call timer_lib_in('shear') @@ -26,7 +26,7 @@ subroutine cgyro_globalshear(ij) #elif defined(_OPENACC) !$acc parallel loop collapse(4) gang vector & !$acc& private(ivc,ir,l,iccj,j,ll,rl,llnt,h1,h2) & -!$acc& present(rhs(:,:,:,ij),omega_ss,omega_sbeta,field,h_x,cap_h_c,c_wave) +!$acc& present(rhs(:,:,:,ij),omega_ss,omega_sbeta,field,cap_h_c,c_wave) #else !$omp parallel do collapse(4) private(ivc,ir,l,iccj,j,ll,rl,llnt,h1,h2) #endif @@ -48,22 +48,22 @@ subroutine cgyro_globalshear(ij) if ( (ir+ll) <= n_radial ) then ! ExB shear - h1 = omega_eb_base*itor*h_x(iccj+llnt,ivc,itor) - ! omega_star shear - h1 = h1-sum(omega_ss(:,iccj+llnt,ivc,itor)*field(:,iccj+llnt,itor)) + h1 = omega_eb_base*itor*cap_h_c(iccj+llnt,ivc,itor) ! beta_star shear h1 = h1-omega_sbeta(iccj+llnt,ivc,itor)*cap_h_c(iccj+llnt,ivc,itor) + ! omega_star shear + h1 = h1-sum(omega_ss(:,iccj+llnt,ivc,itor)*field(:,iccj+llnt,itor)) else h1 = 0.0 endif if ( (ir-ll) >= 1 ) then ! ExB shear - h2 = omega_eb_base*itor*h_x(iccj-llnt,ivc,itor) - ! omega_star shear - h2 = h2-sum(omega_ss(:,iccj-llnt,ivc,itor)*field(:,iccj-llnt,itor)) + h2 = omega_eb_base*itor*cap_h_c(iccj-llnt,ivc,itor) ! beta_star shear h2 = h2-omega_sbeta(iccj-llnt,ivc,itor)*cap_h_c(iccj-llnt,ivc,itor) + ! omega_star shear + h2 = h2-sum(omega_ss(:,iccj-llnt,ivc,itor)*field(:,iccj-llnt,itor)) else h2 = 0.0 endif diff --git a/cgyro/src/cgyro_init_arrays.F90 b/cgyro/src/cgyro_init_arrays.F90 index 9a3901fe7..9cc8e39eb 100644 --- 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 = sbeta_star(is)*energy(ie)/bmag(it)**2 + sb = sbeta(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 index 798fc0f15..bff4f817c 100644 --- a/cgyro/src/cgyro_make_profiles.F90 +++ b/cgyro/src/cgyro_make_profiles.F90 @@ -128,7 +128,7 @@ subroutine cgyro_make_profiles 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) + sbeta(1:n_species) = sbeta_loc(1:n_species) if (ae_flag == 1) then is_ele = n_species+1 @@ -418,7 +418,10 @@ subroutine cgyro_make_profiles case (1) call cgyro_info('ExB shear: Hammett discrete shift (do not use for production simulations)') case (2) - call cgyro_info('ExB shear: Wavenumber advection') + call cgyro_info('ExB shear: Wavenumber advection (h)') + source_flag = 1 + case (3) + call cgyro_info('ExB shear: Wavenumber advection (H)') source_flag = 1 case default call cgyro_error('Unknown ExB shear method') @@ -429,12 +432,13 @@ subroutine cgyro_make_profiles call cgyro_info('ExB shear: OFF') endif - if (profile_shear_flag == 1) then - call cgyro_info('Profile shear: Wavenumber advection') + if (global_flag == 1) then + call cgyro_info('Global profile variation: ON') source_flag = 1 else sdlnndr(1:n_species) = 0.0 sdlntdr(1:n_species) = 0.0 + sbeta(1:n_species) = 0.0 endif diff --git a/cgyro/src/cgyro_read_input.f90 b/cgyro/src/cgyro_read_input.f90 index c95124f36..be77de7dd 100644 --- a/cgyro/src/cgyro_read_input.f90 +++ b/cgyro/src/cgyro_read_input.f90 @@ -82,9 +82,9 @@ subroutine cgyro_read_input call cgyro_readbc_int(hiprec_flag,'HIPREC_FLAG') call cgyro_readbc_int(udsymmetry_flag,'UDSYMMETRY_FLAG') call cgyro_readbc_int(shear_method,'SHEAR_METHOD') + call cgyro_readbc_int(global_flag,'GLOBAL_FLAG') call cgyro_readbc_int(n_global,'N_GLOBAL') call cgyro_readbc_real(nu_global) - call cgyro_readbc_int(profile_shear_flag,'PROFILE_SHEAR_FLAG') call cgyro_readbc_int(theta_plot,'THETA_PLOT') call cgyro_readbc_int(gpu_bigmem_flag,'GPU_BIGMEM_FLAG') call cgyro_readbc_int(upwind_single_flag,'UPWIND_SINGLE_FLAG') @@ -129,7 +129,7 @@ subroutine cgyro_read_input call cgyro_readbc_realv(dlntdr,is) call cgyro_readbc_realv(sdlnndr,is) call cgyro_readbc_realv(sdlntdr,is) - call cgyro_readbc_realv(sbeta_star,is) + call cgyro_readbc_realv(sbeta,is) call cgyro_readbc_realv(dlnndr_scale,is) call cgyro_readbc_realv(dlntdr_scale,is) diff --git a/cgyro/src/cgyro_rhs.F90 b/cgyro/src/cgyro_rhs.F90 index c47a6038f..9eb18a87d 100644 --- a/cgyro/src/cgyro_rhs.F90 +++ b/cgyro/src/cgyro_rhs.F90 @@ -336,10 +336,14 @@ subroutine cgyro_rhs(ij,update_cap) call cgyro_rhs_comm_test ! Wavenumber advection shear terms - ! depends on h_x and field only, updates rhs - call cgyro_advect_wavenumber(ij) - !call cgyro_globalshear(ij) - + if (shear_method == 2) then + ! s* Phi + gamma_E h + call cgyro_advect_wavenumber(ij) + else if (shear_method == 3) then + ! s* Phi + (gamma_E + sbeta) H + call cgyro_globalshear(ij) + endif + call cgyro_rhs_comm_test ! Nonlinear evaluation [f,g] diff --git a/cgyro/src/cgyro_write_initdata.f90 b/cgyro/src/cgyro_write_initdata.f90 index 1c95d8192..bdb5e14f8 100644 --- a/cgyro/src/cgyro_write_initdata.f90 +++ b/cgyro/src/cgyro_write_initdata.f90 @@ -130,14 +130,14 @@ subroutine cgyro_write_initdata ! k(x) = k(0) + s * r/rho ! ! where x = r/rho. The half-domain has -L/4 < r < L/4. - if (profile_shear_flag == 1) then + if (global_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 | 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),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) + is,sdlnndr(is),dlnndr(is)-dn,dlnndr(is)+dn,sdlntdr(is),dlntdr(is)-dt,dlntdr(is)+dt,sbeta(is) enddo endif