From 4586368ec6b264144be45db5076cb5b65fe4c313 Mon Sep 17 00:00:00 2001 From: jonatanschatzlmayr Date: Tue, 28 Jan 2025 15:50:04 +0100 Subject: [PATCH] refine mesh around magnetic axis for small values of sfc_s_min --- SRC/points_2d.f90 | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/SRC/points_2d.f90 b/SRC/points_2d.f90 index 885ec91..375821f 100644 --- a/SRC/points_2d.f90 +++ b/SRC/points_2d.f90 @@ -34,6 +34,9 @@ subroutine create_points_2d(inp_label,n_theta,points,points_s_theta_phi,r_scalin double precision :: s,psi,theta,q,dq_ds,sqrtg,bmod,dbmod_dtheta,R,dR_ds,dR_dtheta,Z,dZ_ds,dZ_dtheta ! integer :: n_center_point + + integer :: n + double precision :: s_ratio, r_frac_n ! n_center_point = 1 if (present(repeat_center_point)) then @@ -53,8 +56,25 @@ subroutine create_points_2d(inp_label,n_theta,points,points_s_theta_phi,r_scalin ! ! loads points that are calculated in preload_for_SYNCH into globals from magdata_in_symfluxcoor_mod and spline interpolats them ! ! this is done so we can then call magdata_in_symfluxcoord_ext on this interpolated fine grid. ! call load_magdata_in_symfluxcoord() + + +! start: make first few entries of r_frac non-equidistant in order to avoid very thin triangles close to the magnetic axis + s_ratio = (1.d0-s_min)/((dble(nlabel))*s_min) + n = int(abs(log(s_ratio)/log(10d0))) + + if ((n.gt.1).and.(nlabel.gt.(n+1))) then + r_frac_n = s_min + (1.d0-s_min)/(dble(nlabel-n+1)) + do i = 1,n-1 + r_frac(i) = exp(log(s_min) + dble(i)*(log(r_frac_n)-log(s_min))/dble(n)) + enddo + endif +! end: make first few entries of r_frac non-equidistant in order to avoid very thin triangles close to the magnetic axis +! when using this version, be sure to comment out the line sstarting with "r_frac =" instead of "r_frac(n:nlabel) =" ! - r_frac = s_min + [(dble(i)*(1.d0-s_min), i=1, size(r_frac), 1)] / (dble(nlabel)) + r_frac(n:nlabel) = s_min + [(dble(i)*(1.d0-s_min), i=1, nlabel-n+1, 1)] / (dble(nlabel-n+1)) + + + !r_frac = s_min + [(dble(i)*(1.d0-s_min), i=1, nlabel, 1)] / (dble(nlabel)) ! !r_frac = [(i, i = 1, nlabel,1)] / dfloat(nlabel) if (present(r_scaling_func)) r_frac = r_scaling_func(r_frac)