flux surface example working

Author: jonathanschilling

.gitignore

@@ -31,6 +31,8 @@ test_1d_rodft11_*
 axis_R.dat
 axis_R_symm.dat
 axis_R_halfGrid.dat
+lcfs_R.dat
+lcfs_Z.dat
 
 # auxiliary tests
 test_util

README.md

@@ -1593,7 +1593,12 @@ Inserting this into the Fourier series for, e.g., the *R* coordinate, leads to a
 The *R* coordinate is a real-valued quantity
 which implies that a two-dimensional `c2r` DFT provided by FFTW can be used to
 perform the backward transform in order to evaluate the flux surface geometry.
-
+One further issue consists in the fact that the definition
+of the Fourier geometry employed in VMEC uses an angle argument *m θ - n ζ*
+where the two-dimensional DFT written out above uses *m θ + n ζ*.
+The Fourier coefficients coming from VMEC
+have to be inserted into the positions corresponding to the reverse sign of *n*
+in the input array given to FFTW in order to resolve this pecularity. 
 
 
 

src/app_flux_surface.c

@@ -47,19 +47,86 @@ void app_flux_surface(int n_theta, int n_zeta) {
         return;
     }
 
+    // number of Fourier modes that can be represented with n_theta and n_zeta
+    // without violating Nyquist criterion
+    int nyq_pol = n_theta / 2 + 1;
+    int nyq_tor = n_zeta / 2 + 1;
 
+    if (nyq_pol < mpol) {
+        printf("n_theta too small\n");
+        return;
+    }
 
+    if (nyq_tor < ntor) {
+        printf("n_zeta too small\n");
+        return;
+    }
 
+    fftw_complex *in_R = fftw_alloc_complex(n_zeta * nyq_pol);
+    fftw_complex *in_Z = fftw_alloc_complex(n_zeta * nyq_pol);
+    double *out_R = fftw_alloc_real(n_zeta * n_theta);
+    double *out_Z = fftw_alloc_real(n_zeta * n_theta);
+
+    fftw_plan p_R = fftw_plan_dft_c2r_2d(n_zeta, n_theta, in_R, out_R, FFTW_ESTIMATE);
+    fftw_plan p_Z = fftw_plan_dft_c2r_2d(n_zeta, n_theta, in_Z, out_Z, FFTW_ESTIMATE);
+
+    fill_zero_2d_cplx(n_zeta, nyq_pol, in_R);
+    fill_zero_2d_cplx(n_zeta, nyq_pol, in_Z);
+
+    // copy LCFS Fourier coefficients into input array
+    int idx_vmec = 0, idx_in;
+    int m = 0;
+    for (int n = 0; n <= ntor; ++n) {
+
+        // compute target index for input array
+        if (n <= 0) {
+            idx_in = -n * nyq_pol + m;
+        } else {
+            idx_in = (n_zeta - n) * nyq_pol + m;
+        }
+
+        // VMEC: for m=0, only positive ntor --> summed up within VMEC definition with negative-n --> 2*0.5 = 1
+        // also for sin-parity quantities: m=0 would imply (m theta -n zeta) == (-n zeta), but VMEC uses (n zeta) for m=0
+        // --> does not matter for cos-parity since cos(-n zeta) = cos(n zeta), but sin(-n zeta) = -sin(n zeta)
+        // and hence the additional -1 cancels out by compensating the -1 from the complex multiply in the DFT definition
+        in_R[idx_in] = rmnc[idx_vmec];
+        in_Z[idx_in] = I * zmns[idx_vmec];
+
+        idx_vmec++;
+    }
 
+    for (m = 1; m < mpol; ++m) {
+        for (int n = -ntor; n <= ntor; ++n) {
 
+            // compute target index for input array
+            if (n <= 0) {
+                idx_in = -n * nyq_pol + m;
+            } else {
+                idx_in = (n_zeta - n) * nyq_pol + m;
+            }
 
+            // scale coefficients by 0.5 since FFTW implies
+            // that coeffs for m<0 were present in the logically equivalent DFT input
+            // which is not the case for VMEC output
+            in_R[idx_in] = 0.5 * rmnc[idx_vmec];
+            in_Z[idx_in] = -0.5 * I * zmns[idx_vmec];
 
+            idx_vmec++;
+        }
+    }
 
+    fftw_execute(p_R);
+    fftw_execute(p_Z);
 
+    dump_2d_real("lcfs_R.dat", n_zeta, n_theta, out_R);
+    dump_2d_real("lcfs_Z.dat", n_zeta, n_theta, out_Z);
 
-
-
-
+    fftw_destroy_plan(p_R);
+    fftw_destroy_plan(p_Z);
+    fftw_free(in_R);
+    fftw_free(in_Z);
+    fftw_free(out_R);
+    fftw_free(out_Z);
 }
 
 int main(int argc, char** argv) {

src/plot_lcfs_realspace.py

@@ -0,0 +1,36 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 16 17:53:03 2021
+
+@author: jonathan
+"""
+
+import os
+import numpy as np
+from mayavi import mlab
+
+folder = ".."
+
+# load output from FFTW example 'app_flux_surface.c'
+lcfs_R = np.loadtxt(os.path.join(folder, "lcfs_R.dat"))
+lcfs_Z = np.loadtxt(os.path.join(folder, "lcfs_Z.dat"))
+n_zeta, n_theta = np.shape(lcfs_R)
+
+# assume W7-X
+nfp = 5
+
+# regular grid in phi = zeta/nfp
+phiGrid = 2.0*np.pi*np.arange(n_zeta)/(n_zeta*nfp)
+
+# compute X,Y from R,phi
+lcfs_X = np.zeros([n_zeta, n_theta])
+lcfs_Y = np.zeros([n_zeta, n_theta])
+for i in range(n_zeta):
+    for j in range(n_theta):
+        lcfs_X[i,j] = lcfs_R[i,j]*np.cos(phiGrid[i])
+        lcfs_Y[i,j] = lcfs_R[i,j]*np.sin(phiGrid[i])
+
+# plot real-space geometry of flux surface using Mayavi
+mlab.mesh(lcfs_X, lcfs_Y, lcfs_Z)
+mlab.show()

src/util.h

@@ -9,6 +9,34 @@
 #include <complex.h>
 #include <fftw3.h>
 
+void fill_zero_1d_real(int n, double *arr) {
+    for (int i = 0; i < n; ++i) {
+            arr[i] = 0.0;
+        }
+}
+
+void fill_zero_1d_cplx(int n, fftw_complex *arr) {
+    for (int i = 0; i < n; ++i) {
+        arr[i] = 0.0;
+    }
+}
+
+void fill_zero_2d_real(int rows, int cols, double *arr) {
+    for (int i = 0; i < rows; ++i) {
+        for (int j = 0; j < cols; ++j) {
+            arr[i * cols + j] = 0.0;
+        }
+    }
+}
+
+void fill_zero_2d_cplx(int rows, int cols, fftw_complex *arr) {
+    for (int i = 0; i < rows; ++i) {
+        for (int j = 0; j < cols; ++j) {
+            arr[i * cols + j] = 0.0;
+        }
+    }
+}
+
 void fill_random_1d_real(int n, double *arr) {
     for (int i = 0; i < n; ++i) {
         arr[i] = rand() / ((double) RAND_MAX);