submodule(iol_source_m) iol_source_init_s !! Routines related with initialisation of IOL sources implicit none contains module subroutine init_iol(self, filename, equi, mesh, polars) character(len=*), intent(in) :: filename class(iol_source_t), intent(inout) :: self class(equilibrium_t), intent(inout) :: equi type(mesh_cart_t), intent(in) :: mesh type(polars_t), intent(in) :: polars ! Read parameters call self%read_params_iol(filename, equi) ! Create IOL grid call self%grid_init(equi, polars) ! Create IOL spline grid call self%spline_grid_init(equi, mesh, polars) ! Create Xplane call self%xplane_init(equi, mesh) ! Account for flipped equilibra as IOL can only handle lower single null select type(equi) type is(numerical_t) if (equi%flip_z) then call self%flip() endif end select ! Allocate fields allocate(self%dens(self%jmax, self%imax), source=0.0_GP) allocate(self%temp_ion(self%jmax, self%imax), source=0.0_GP) allocate(self%pot(self%jmax, self%imax), source=0.0_GP) allocate(self%pot_sp(self%nrho_sp), source=0.0_GP) ! Allocate iol source on inner points of GRILLIX grid allocate(self%src_iol(mesh%get_n_points_inner()), source=0.0_GP) ! Reset time averaging self%nt_av = 0 end subroutine module subroutine read_params_iol(self, filename, equi) character(len=*), intent(in) :: filename class(iol_source_t), intent(inout) :: self class(equilibrium_t), intent(inout) :: equi ! to be read from some parameterfile integer :: signBt !! Sign of toroidal field following AUG conventions integer :: signIp !! Sign of toroidal current following AUG conventions real(GP) :: xnrm !! Normalisation length within poloidal plane [m] real(GP) :: bnrm !! Normalisation value for magnetic field [T] real(GP) :: x_xpoint !! x-coordinate of X-point (Normalised to xnrm) real(GP) :: y_xpoint !! y-coordinate of X-point (Normalised to xnrm) real(GP) :: psix !! Poloidal flux at separatrix [T/m^2] real(GP) :: psi0 !! Poloidal flux at magnetic axis [T/m^2] real(GP) :: xmin_m !! Slope of linear fit for xmin(psi) w/ xmin in [R_0] !! and psi in [Wb/2pi] [2pi R_0/Wb] !! TODO: WHAT ARE THE UNITS HERE, HOWW TO COMPUTE THIS ??? @Robert ???? real(GP) :: xmin_b !! Intercept for linear fit for xmin(psi) w/ xmin in [R_0] !! and psi in [Wb/2pi] [R_0] !! TODO: WHAT ARE THE UNITS HERE, HOWW TO COMPUTE THIS ??? @Robert ???? real(GP) :: invs_aspc !! Aspect ratio a/R_0 real(GP) :: ion_mass !! Ion mass with respect to hydrogen real(GP) :: ion_charge !! ion charge with respect to elementary charge real(GP) :: er_mult !! Scaling factor for radial electric field real(GP) :: ti_mult !! Scaling factor for temperature real(GP) :: ne_nrm !! Normalisation of density [m^-3] real(GP) :: ti_nrm !! Normalisation of ion temperature [eV] real(GP) :: tau_nrm !! Time normalisation R_0/c_s [s] real(GP) :: tau_av_frame !! Time that fields are averaged over in GRILLIX units i.e. [R_0/c_s] logical :: init_continue !! Switch if simulation continued from existing state integer :: kmax !! Number of points in pitch angle space integer :: lmax !! Number of points in energy space integer :: mmax !! Number of points to consider as IOL for R_loss > R_X integer :: nsep_corr_cw !! Number of points that are corrected around the X-point !! in clockwise direction at finding separatrix integer :: nsep_corr_ccw !! Number of points that are corrected around the X-point !! in counter clockwise direction at finding separatrix integer :: io_error character(len=256) :: io_errmsg namelist / iol_params / signBT, signIP, xnrm, bnrm, x_xpoint,y_xpoint, & psix, psi0, xmin_m, xmin_b, & invs_aspc, ion_mass, ion_charge, & er_mult, ti_mult, ne_nrm, ti_nrm,tau_nrm namelist / iol_grid_params / kmax, lmax, mmax, & nsep_corr_cw, nsep_corr_ccw namelist / iol_tau_average / tau_av_frame, init_continue open(unit = 20, file = filename, status = 'old', action = 'read',& iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif ! read iol_params signBt = 1 signIp = 1 xnrm = GP_NAN bnrm = GP_NAN x_xpoint = GP_NAN y_xpoint = GP_NAN psix = GP_NAN psi0 = GP_NAN xmin_m = GP_NAN xmin_b = GP_NAN invs_aspc = GP_NAN ion_mass = 1.0_GP ion_charge = 1.0_GP er_mult = 1.0_GP ti_mult = 1.0_GP ne_nrm = GP_NAN ti_nrm = GP_NAN tau_nrm = GP_NAN read(20, nml = iol_params, iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif self%bnrm = bnrm self%signBt = signBt self%signIp = signIp self%R_magaxis = xnrm self%Rmin_m = xnrm * xmin_m self%Rmin_b = xnrm * xmin_b self%RX = xnrm * x_xpoint self%zX = xnrm * y_xpoint self%BphiX = signBt * equi%btor(x_xpoint, y_xpoint, phi_axsym) * bnrm self%psix = psix self%psirange_conf_red = psi0 - psix self%invs_aspc = invs_aspc self%mass = ion_mass * mass_hydrogen_SI self%charge = ion_charge * elementary_charge_SI self%er_mult = er_mult self%ti_mult = ti_mult self%ne_nrm = ne_nrm self%ti_nrm = ti_nrm self%tau_nrm = tau_nrm ! read iol_grid_params kmax = 0 lmax = 0 mmax = 0 nsep_corr_cw = 0 nsep_corr_ccw = 0 read(20, nml = iol_grid_params, iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif self%kmax = kmax self%lmax = lmax self%mmax = mmax self%nsep_corr_cw = nsep_corr_cw self%nsep_corr_ccw = nsep_corr_ccw ! read iol_tau_average parameters tau_av_frame = GP_LARGEST init_continue = .false. read(20, nml = iol_tau_average, iostat = io_error, iomsg = io_errmsg) if (io_error /= 0) then call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__) endif self%tau_av_frame = tau_av_frame self%init_continue = init_continue ! close file close(20) end subroutine module subroutine grid_init(self, equi, polars) class(iol_source_t), intent(inout) :: self class(equilibrium_t), intent(inout) :: equi type(polars_t), intent(in) :: polars integer :: jp, ip, isp real(GP) :: rho, theta, xnrm, ynrm, xp, yp, psi0 real(GP) :: xx_nrm, yx_nrm real(GP), allocatable, dimension(:) :: thetas ! Set iol grid parameters from polar grid self%imax = polars%grid%get_nrho() self%jmax = polars%grid%get_ntheta() ! Build iol grid allocate(self%Rcoord(self%jmax, self%imax)) allocate(self%zcoord(self%jmax, self%imax)) allocate(self%Babs(self%jmax, self%imax)) allocate(self%Bphi(self%jmax, self%imax)) allocate(self%psi(self%jmax, self%imax)) psi0 = self%psirange_conf_red + self%psix do ip = 1, self%imax do jp = 1, self%jmax rho = polars%grid%get_rho(ip) theta = polars%grid%get_theta(jp) call polar_to_cart(equi, rho, theta, phi_axsym, xnrm, ynrm) xp = xnrm * self%R_magaxis yp = ynrm * self%R_magaxis self%Rcoord(jp, ip) = xp self%zcoord(jp, ip) = yp self%Babs(jp,ip) = equi%absb(xnrm, ynrm, phi_axsym) * self%bnrm self%Bphi(jp,ip) = self%bnrm * self%signBt * & equi%btor(xnrm, ynrm, phi_axsym) ! ??? Does sign convention match with GRILLIX ??? self%psi(jp, ip) = rho**2 * (self%psix-psi0)+psi0 - self%psix enddo enddo ! Compute safety factor allocate(self%safetyf(self%imax)) do ip = 1, self%imax rho = polars%grid%get_rho(ip) self%safetyf(ip) = safety_factor(equi, rho) enddo ! Get separatrix allocate(self%sep_Rcoord(self%jmax)) allocate(self%sep_Zcoord(self%jmax)) xx_nrm = self%RX/self%R_magaxis yx_nrm = self%ZX/self%R_magaxis allocate(thetas(self%jmax)) do jp = 1, self%jmax thetas(jp) = polars%grid%get_theta(jp) enddo call find_separatrix(equi, xx_nrm, yx_nrm, thetas, & self%sep_Rcoord, self%sep_Zcoord, & self%idx_X, self%nsep_corr_cw + 1, self%nsep_corr_ccw + 1) do jp = 1, self%jmax self%sep_Rcoord(jp) = self%sep_Rcoord(jp)*self%R_magaxis self%sep_Zcoord(jp) = self%sep_Zcoord(jp)*self%R_magaxis enddo deallocate(thetas) self%idx_imp = minloc(self%sep_Rcoord, dim = 1) ! Get idx_imp: position of minimum R on separatrix ! Compute poloidal distance to X-point allocate(self%sep_cw(self%jmax)) allocate(self%sep_ccw(self%jmax)) do jp =1, self%jmax self%sep_cw(jp) = integrate_to_xpoint(jp, .false.) self%sep_ccw(jp) = integrate_to_xpoint(jp, .true. ) enddo ! Compute modified poloidal distances to X-point allocate(self%sep_cw_mod(self%jmax)) allocate(self%sep_ccw_mod(self%jmax)) do jp =1, self%jmax self%sep_cw_mod(jp) = integrate_to_xpoint_mod(jp, .false.) self%sep_ccw_mod(jp) = integrate_to_xpoint_mod(jp, .true.) enddo ! Correct entries at X-pint for sep_.. and sep_..._mod self%sep_cw(self%idx_X) = GP_LARGEST self%sep_cw(self%idx_X) = minval(self%sep_cw) self%sep_ccw(self%idx_X) = GP_LARGEST self%sep_ccw(self%idx_X) = minval(self%sep_ccw) self%sep_cw_mod(self%idx_X) = GP_LARGEST self%sep_cw_mod(self%idx_X) = minval(self%sep_cw_mod) self%sep_ccw_mod(self%idx_X) = GP_LARGEST self%sep_ccw_mod(self%idx_X) = minval(self%sep_ccw_mod) contains real(GP) function integrate_to_xpoint(jp, counter) !! Computes poloidal distance from given poloidal position (jp) !! to X-point (idx_X) in co or counter-clockwisedirection integer, intent(in) :: jp !! Poloidal index for which distance is computed logical, intent(in) :: counter !! Integration in counterclockwise direction (false: co-clockwise) logical :: found integer :: jit, k real(GP) :: dst, x0, y0, x1, y1 found = .false. jit = jp dst = 0.0_GP intloop: do k = 1, self%jmax if (jit == self%idx_X) then found = .true. exit intloop endif x0 = self%sep_Rcoord(jit) / self%R_magaxis y0 = self%sep_Zcoord(jit) / self%R_magaxis if (counter) then jit = jit + 1 if (jit > self%jmax) then jit = jit - self%jmax endif else jit = jit - 1 if (jit < 1) then jit = jit + self%jmax endif endif x1 = self%sep_Rcoord(jit) / self%R_magaxis y1 = self%sep_Zcoord(jit) / self%R_magaxis dst = dst + sqrt((x1-x0)**2+(y1-y0)**2) enddo intloop if (found) then integrate_to_xpoint = dst * self%R_magaxis else call handle_error('Integrate_to_xpoint not successful', & GRILLIX_ERR_OTHER, __LINE__, __FILE__) endif end function real(GP) function integrate_to_xpoint_mod(jp, counter) !! Computes modified poloidal distance !! from given poloidal position (jp) !! to X-point (idx_X) in co or counter-clockwisedirection !! Ask @Robert Brzozowski for what is actually computed here integer, intent(in) :: jp !! Poloidal index for which distance is computed logical, intent(in) :: counter !! Integration in counterclockwise direction (false: co-clockwise) logical :: found integer :: jit, k real(GP) :: dst, x0, y0, x1, y1, bpol0, rmin found = .false. jit = jp dst = 0.0_GP intloopmod: do k = 1, self%jmax if (jit == self%idx_X) then found = .true. exit intloopmod endif x0 = self%Rcoord(jit, self%imax) / self%R_magaxis y0 = self%zcoord(jit, self%imax) / self%R_magaxis if (counter) then jit = jit + 1 if (jit > self%jmax) then jit = jit - self%jmax endif else jit = jit - 1 if (jit < 1) then jit = jit + self%jmax endif endif x1 = self%Rcoord(jit, self%imax) / self%R_magaxis y1 = self%zcoord(jit, self%imax) / self%R_magaxis rmin = sqrt( (x0 -equi%x0)**2 + (y0-equi%y0)**2) bpol0 = sqrt( equi%bx(x0,y0,phi_axsym)**2 + & equi%by(x0,y0,phi_axsym)**2 ) dst = dst + sqrt( (x1 -x0)**2 + (y1 - y0)**2 ) & * 1.0_GP / (x0 * bpol0) * (rmin / self%safetyf(self%imax) ) enddo intloopmod if (found) then integrate_to_xpoint_mod = dst * self%R_magaxis else call handle_error('Integrate_to_xpoint not successful', & GRILLIX_ERR_OTHER, __LINE__, __FILE__) endif end function end subroutine module subroutine spline_grid_init(self, equi, mesh, polars) class(iol_source_t), intent(inout) :: self class(equilibrium_t), intent(inout) :: equi type(mesh_cart_t), intent(in) :: mesh type(polars_t), intent(in) :: polars real(GP) :: x, y, rho, theta_idx_imp, x_idx_imp, y_idx_imp integer :: lomp, lc, k, nrho_sp, l ! Get mesh point closest to outboard midplane core call polar_to_cart(equi, polars%grid%get_rhopol_min(), 0.0_GP, & phi_axsym, x, y) call mesh%get_surrounding_indices(x, y, 2, ind_nearest = lomp) lc = lomp do while (lc > 0) lc = mesh%index_neighbor(-1,0,lc) if (lc>0) then lomp = lc endif enddo ! Get indices of outboard midplane mesh points lc = lomp nrho_sp = 1 do while (lc > 0) lc = mesh%index_neighbor(1,0,lc) if (lc > 0) then nrho_sp = nrho_sp + 1 endif enddo self%nrho_sp = nrho_sp allocate(self%indmesh_sp(nrho_sp)) lc = lomp self%indmesh_sp(1) = lomp k = 1 do while (lc > 0) lc = mesh%index_neighbor(1,0,lc) if (lc > 0) then nrho_sp = nrho_sp + 1 k = k+1 self%indmesh_sp(k) = lc endif enddo ! Compute corresponding rho allocate(self%rhopts_sp(self%nrho_sp)) do k = 1, self%nrho_sp l = self%indmesh_sp(k) x = mesh%get_x(l) y = mesh%get_y(l) rho = equi%rho(x,y,phi_axsym) self%rhopts_sp(k) = rho enddo ! Compute Rmin corresponding to flux surface allocate(self%Rmincontour_sp(self%nrho_sp), source = 0.0_GP) do k = 1, self%nrho_sp l = self%indmesh_sp(k) x = mesh%get_x(l) y = mesh%get_y(l) rho = equi%rho(x,y,phi_axsym) theta_idx_imp= polars%grid%get_theta(self%idx_imp) call polar_to_cart(equi, rho, theta_idx_imp, & phi_axsym, x_idx_imp, y_idx_imp) self%Rmincontour_sp(k) = x_idx_imp*self%R_magaxis enddo end subroutine module subroutine xplane_init(self, equi, mesh) class(iol_source_t), intent(inout) :: self class(equilibrium_t), intent(inout) :: equi type(mesh_cart_t), intent(in) :: mesh real(GP) :: rx_nrm, zx_nrm, x, y, psi0 integer :: lx, jx, lc, ixp_min, ixp_max, ixp, ixpn ! Get index of point closest to X-point rx_nrm = self%RX / self%R_magaxis zx_nrm = self%zX / self%R_magaxis call mesh%get_surrounding_indices(rx_nrm, zx_nrm, 2, ind_nearest = lx) ! get vertical index of X-point jx = mesh%get_cart_j(lx) ! get minumum and maximum horizontal index of plane ixp_min = mesh%get_cart_i(lx) lc = lx do while (lc > 0) ixp_min = ixp_min - 1 lc = mesh%index_neighbor(-1,0,lc) enddo ixp_min = ixp_min+1 ixp_max = mesh%get_cart_i(lx) lc = lx do while (lc > 0) ixp_max = ixp_max + 1 lc = mesh%index_neighbor(1,0,lc) enddo ixp_max = ixp_max-1 self%xplane_n = ixp_max - ixp_min +1 allocate(self%xplane_R(self%xplane_n)) allocate(self%xplane_Babs(self%xplane_n)) allocate(self%xplane_bphi(self%xplane_n)) allocate(self%xplane_psi(self%xplane_n)) psi0 = self%psirange_conf_red + self%psix do ixp = ixp_min, ixp_max ixpn = ixp - ixp_min+1 lc = mesh%cart_to_mesh_index(ixp, jx) x = mesh%get_x(lc) y = mesh%get_y(lc) self%xplane_R(ixpn) = x*self%R_magaxis self%xplane_Babs(ixpn) = equi%absb(x, y, phi_axsym) * self%bnrm self%xplane_bphi(ixpn) = self%bnrm * self%signBt * & equi%btor(x, y, phi_axsym) self%xplane_psi(ixpn) = equi%rho(x, y, phi_axsym)**2 * & (self%psix - psi0) + psi0 - self%psix enddo end subroutine module subroutine flip(self) class(iol_source_t), intent(inout) :: self integer :: ip real(GP), dimension(self%jmax) :: wrk ! Flip vertically self%zcoord = -self%zcoord ! Reshuffle for counterclockwise orientation of data do ip = 1, self%imax call flip_array_direction(self%jmax, self%Rcoord(:,ip)) call flip_array_direction(self%jmax, self%zcoord(:,ip)) call flip_array_direction(self%jmax, self%Babs(:,ip)) call flip_array_direction(self%jmax, self%Bphi(:,ip)) call flip_array_direction(self%jmax, self%psi(:,ip)) enddo call flip_array_direction(self%jmax, self%sep_Rcoord) call flip_array_direction(self%jmax, self%sep_cw) call flip_array_direction(self%jmax, self%sep_cw_mod) call flip_array_direction(self%jmax, self%sep_ccw) call flip_array_direction(self%jmax, self%sep_ccw_mod) ! Update indices of IMP and X-point accordingly self%idx_imp = self%jmax - self%idx_imp + 2 self%idx_X = self%jmax - self%idx_X + 2 ! Exchange sep_cw with sep_ccw wrk = self%sep_cw self%sep_cw = self%sep_ccw self%sep_ccw = wrk wrk = self%sep_cw_mod self%sep_cw_mod = self%sep_ccw_mod self%sep_ccw_mod = wrk end subroutine subroutine find_separatrix(equi, x_xpoint, y_xpoint, & thetas, xsep, ysep, jx, ncm, ncp) !! Finds cartesian coordinates of separatrix & !! for list of given geometric poloidal angles class(equilibrium_t), intent(inout) :: equi !! Equilibrium real(GP), intent(in) :: x_xpoint !! x-coordinate of X-point real(GP), intent(in) :: y_xpoint !! y-coordinate of X-point real(GP), dimension(:), intent(inout) :: thetas !! Geometric poloidal angles, will be slightly modified near X-point real(GP), dimension(size(thetas)), intent(out) :: xsep !! X-coordinates of separatrix real(GP), dimension(size(thetas)), intent(out) :: ysep !! y-coordinates of separatrix integer, intent(out) :: jx !! Index of X-point integer, intent(in), optional :: ncm !! Number of points - 1 that are corrected around the X-point !! in clockwise direction integer, intent(in), optional :: ncp !! Number of points - 1 that are corrected around the X-point !! in counter clockwise direction integer :: ntheta, jp real(GP) :: rho_dummy, rho, theta_x, dst, dst_x integer :: jcm, jcp, ncma, ncpa ncma = 0 if (present(ncm)) then ncma = ncm endif ncpa = 0 if (present(ncp)) then ncpa = ncp endif ntheta = size(thetas) ! Find and set X-point call cart_to_polar(equi, x_xpoint, y_xpoint, phi_axsym, rho_dummy, theta_x) dst_x = GP_LARGEST do jp = 1, ntheta dst = abs(thetas(jp)-theta_x) if (dst < dst_x) then jx = jp dst_x = dst endif enddo xsep(jx) = x_xpoint ysep(jx) = y_xpoint thetas(jx) = theta_x rho = 1.0_GP do jp = 1, ntheta if (jp == jx) then cycle elseif ((jp > jx - ncma) .and. jp < (jx + ncpa)) then ! There might be issues to find the points near the X-point, ! so exclude these yet and determine them only lateron (directly after this loop) cycle else call polar_to_cart(equi, rho, thetas(jp), phi_axsym, xsep(jp), ysep(jp)) endif enddo ! Correct for points near X-point, ! put those along straight line between first reasonable point and X-point jcm = jx - ncma do jp = jcm+1, jx-1 xsep(jp) = xsep(jcm) + (xsep(jx)-xsep(jcm))*(jp - jcm)/(1.0_GP*(ncma)) ysep(jp) = ysep(jcm) + (ysep(jx)-ysep(jcm))*(jp - jcm)/(1.0_GP*(ncma)) call cart_to_polar(equi, xsep(jp), ysep(jp), & phi_axsym, rho_dummy, thetas(jp)) enddo jcp = jx + ncpa do jp = jx+1, jcp - 1 xsep(jp) = xsep(jx) + (xsep(jcp)-xsep(jx))*(jp - jx)/(1.0_GP*(ncpa)) ysep(jp) = ysep(jx) + (ysep(jcp)-ysep(jx))*(jp - jx)/(1.0_GP*(ncpa)) call cart_to_polar(equi, xsep(jp), ysep(jp), & phi_axsym, rho_dummy, thetas(jp)) enddo end subroutine end submodule