lazersos
commented
6 years ago In the STELLOPT virtual_casing_mod.f90 module used in DIAGNO, FIELDLINES, and BEAMS3D, we construct splines over the surface of the required quantities. This alleviates some of the overhead of evaluating the trigonometric functions but isn't necessarily optimal for this problem.
To calculate the field from the virtual casing we need three vector quantities. The field on the plasma vacuum interface (or rather the related surface current), the position where we want to evaluate the field, and the position of the interface. Now in SPEC the second quantity does not change during the simulation so this can be computed once on a finite grid. The position of the interface and the surface current change quite a bit.
Now the DCUHRE algorithm is nice in that it will evaluate the integral to a defined tolerance by refining the mesh of points over which it's sampling. This is good if for example the surface current is approaching the point at which we wish to evaluate the field. However, this can be problematic as well, as the whole processes can become computationally expensive.
A hack is to apply stellarator symmetry to reduce the number of points at which the field needs to be evaluated. As it turns out for stellarators we only need to evaluate the field on a grid from phi=[0,pi/NFP] and theta=[0,pi]. This is because of the mirror symmetry. Thus compared to a non-stellarator symmetric simulation the number of points at which the field needs to be evaluated is reduced by a factor of 4. BTW, VMEC does this when sampling the vacuum field.
However, this is just a 'trick.' Another method one could employ is to pre-calculate the Fourier terms into a matrix. If we always want to evaluate quantities on a known (u,v) mesh then evaluating the Fourier double integral is just a matrix multiply and only requires an initial evaluate of the trig coefficients, see this following code (from virtual_casing_mod.f90)
SUBROUTINE uvtomn_local(k1,k,mnmax,nu,nv,xu,xv,fmn,xm,xn,f,signs,calc_trig)
IMPLICIT NONE
! INPUT VARIABLES
INTEGER, INTENT(in) :: k1
INTEGER, INTENT(in) :: k
INTEGER, INTENT(in) :: mnmax
INTEGER, INTENT(in) :: nu
INTEGER, INTENT(in) :: nv
DOUBLE PRECISION, INTENT(in) :: xu(1:nu)
DOUBLE PRECISION, INTENT(in) :: xv(1:nv)
DOUBLE PRECISION, INTENT(inout) :: fmn(1:mnmax,k1:k)
INTEGER, INTENT(in) :: xm(1:mnmax)
INTEGER, INTENT(in) :: xn(1:mnmax)
DOUBLE PRECISION, INTENT(in) :: f(1:nu,1:nv,k1:k)
INTEGER, INTENT(in) :: signs
INTEGER, INTENT(in) :: calc_trig
! LOCAL VARIABLES
INTEGER :: mn, i, j, ier, ik
DOUBLE PRECISION :: xn_temp(1:mnmax,1)
DOUBLE PRECISION :: xm_temp(1:mnmax,1)
DOUBLE PRECISION :: pi2_l
DOUBLE PRECISION :: mt(1:mnmax,1:nu)
DOUBLE PRECISION :: nz(1:mnmax,1:nv)
DOUBLE PRECISION :: xu_temp(1,1:nu)
DOUBLE PRECISION :: xv_temp(1,1:nv)
DOUBLE PRECISION, ALLOCATABLE, SAVE :: fnuv(:)
DOUBLE PRECISION, ALLOCATABLE, SAVE :: cosmt(:,:)
DOUBLE PRECISION, ALLOCATABLE, SAVE :: sinmt(:,:)
DOUBLE PRECISION, ALLOCATABLE, SAVE :: cosnz(:,:)
DOUBLE PRECISION, ALLOCATABLE, SAVE :: sinnz(:,:)
! BEGIN SUBROUTINE
pi2_l = 8.0D+0 * ATAN(1.)
IF (calc_trig == 1) THEN
IF (ALLOCATED(cosmt)) DEALLOCATE(cosmt)
IF (ALLOCATED(sinmt)) DEALLOCATE(sinmt)
IF (ALLOCATED(cosnz)) DEALLOCATE(cosnz)
IF (ALLOCATED(sinnz)) DEALLOCATE(sinnz)
IF (ALLOCATED(fnuv)) DEALLOCATE(fnuv)
ALLOCATE(fnuv(1:mnmax),STAT=ier)
ALLOCATE(cosmt(1:mnmax,1:nu),sinmt(1:mnmax,1:nu),&
cosnz(1:mnmax,1:nv),sinnz(1:mnmax,1:nv),STAT=ier)
FORALL(i=1:mnmax) xm_temp(i,1)=DBLE(xm(i))
FORALL(i=1:mnmax) xn_temp(i,1)=DBLE(xn(i))
FORALL(i=1:mnmax) fnuv(i) = 2.0D+0/DBLE(nu*nv)
WHERE(xm == 0) fnuv = 5.0D-1*fnuv
FORALL(i=1:nu) xu_temp(1,i)=xu(i)
FORALL(i=1:nv) xv_temp(1,i)=xv(i)
mt = MATMUL(xm_temp,xu_temp)
nz = MATMUL(xn_temp,xv_temp)
FORALL(mn=1:mnmax,i=1:nu) cosmt(mn,i) = dcos(pi2_l*mt(mn,i))
FORALL(mn=1:mnmax,i=1:nu) sinmt(mn,i) = dsin(pi2_l*mt(mn,i))
FORALL(mn=1:mnmax,i=1:nv) cosnz(mn,i) = dcos(pi2_l*nz(mn,i))
FORALL(mn=1:mnmax,i=1:nv) sinnz(mn,i) = dsin(pi2_l*nz(mn,i))
END IF
fmn = zero
IF (signs == 0) THEN
DO mn = 1, mnmax
DO i = 1, nu
DO j = 1, nv
fmn(mn,k1:k) = fmn(mn,k1:k) + f(i,j,k1:k)*(cosmt(mn,i)*cosnz(mn,j) &
- sinmt(mn,i)*sinnz(mn,j))*fnuv(mn)
END DO
END DO
END DO
ELSE IF (signs == 1) THEN
DO mn = 1, mnmax
DO i = 1, nu
DO j = 1, nv
fmn(mn,k1:k) = fmn(mn,k1:k) + f(i,j,k1:k)*(sinmt(mn,i)*cosnz(mn,j) &
+ cosmt(mn,i)*sinnz(mn,j))*fnuv(mn)
END DO
END DO
END DO
END IF
! END SUBROUTINE
END SUBROUTINE uvtomn_localIf we get clever about the grid spacing then it may be possible to just to integrate over fixed grids removing the overhead of DCUHRE.
The virtual casing method in SPEC is far from optimal for the problem at hand. We should address this.
Background links
The integration of the equations of propagation of electric waves Use of the virtual-casing principle in calculating the containing magnetic field in toroidal plasma systems The virtual-casing principle for 3D toroidal systems A magnetic diagnostic code for 3D fusion equilibria The virtual-casing principle and Helmholtz's theorem
Background info from @SRHudson
It seems that there was some bug in replacing the previous NAG routine with DCUHRE, but this is not the main problem. If anyone wants to fix this bug (I started, and have committed/pushed my changes, but I need to go home now.)
In SPEC we need to calculate the magnetic field on a regular grid, Nt in theta and Nz in zeta, on the computational boundary. For each one of these Ntz = Nt Nz virtual casing integrals, we need to perform a surface integral over the plasma boundary. I think that the best algorithm in SPEC would be to compute the geometry and tangential field on the plasma boundary just once, and then to use this information for all of the virtual casing integrals required on the Nt Nz grid on the computational boundary. The point is that there is a lot of information that is common, so if we compute the common information only once then SPEC will be faster.
Presently I am using a slow Fourier transform to compute the magnetic field at whatever point that DCUHRE requires, and because DCUHRE is adaptive, the points on the surface where DCUHRE requires the geometry/tangential field of the plasma boundary will vary as the position on the computational boundary varies. One approach, which I first suggested, is to use/develop a quadrature algorithm that only requires the integrand on a fixed regular grid. This will probably be fasted, but we will potentially lose the accuracy that adaptive routines usually guarantee. Another approach is to first compute the geometry/tangential field on the plasma boundary on a regular grid and thereafter compute the geometry/tangential field using splines (or something fast). I am not sure that this will really solve the problem, because if we need such adaptivity it probably means we need to run SPEC at higher Fourier resolution in any case.