edits

CMakeLists.txt

@@ -117,6 +117,11 @@ set_target_properties(${PROJECT_NAME}
     CXX_STANDARD 14
     CXX_STANDARD_REQUIRED ON)
 
+target_include_directories(${PROJECT_NAME} PRIVATE "thirdparty/LinkingNumber/toolbox/src/")
+target_link_directories(${PROJECT_NAME} PRIVATE "thirdparty/LinkingNumber/build/")
+target_link_libraries(${PROJECT_NAME} PRIVATE LinkingNumber_linux)
+
+
 target_include_directories(${PROJECT_NAME} PRIVATE "thirdparty/xtensor/include" "thirdparty/xtensor-python/include" "thirdparty/xsimd/include" "thirdparty/xtl/include" "thirdparty/eigen" ${Python_NumPy_INCLUDE_DIRS} "src/simsoptpp/")
 target_link_libraries(${PROJECT_NAME} PRIVATE fmt::fmt)
 

src/simsopt/geo/boozersurface.py

@@ -62,18 +62,18 @@ def run_code(self, boozer_type, iota, G=None, verbose=True):
             return
 
         if boozer_type == 'exact':
-            res = self.solve_residual_equation_exactly_newton(tol=1e-13, maxiter=20, iota=iota, G=G, verbose=verbose)
+            res = self.solve_residual_equation_exactly_newton(tol=1e-13, maxiter=40, iota=iota, G=G, verbose=verbose)
             return res
 
         elif boozer_type == 'ls':
             # first try BFGS
             #res = self.minimize_boozer_penalty_constraints_ls(tol=1e-10, maxiter=100, constraint_weight=self.constraint_weight, iota=iota, G=G, method='manual')
-            res = self.compute_boozerls_BFGS(tol=1e-10, maxiter=1e3, constraint_weight=self.constraint_weight, iota=iota, G=G, verbose=verbose)
+            res = self.compute_boozerls_BFGS(tol=1e-10, maxiter=1500, constraint_weight=self.constraint_weight, iota=iota, G=G, verbose=verbose)
             iota, G = res['iota'], res['G']
 
             ## polish off using Newton's method
             self.need_to_run_code = True
-            res = self.compute_boozerls_newton(tol=1e-11, maxiter=20, constraint_weight=self.constraint_weight, iota=iota, G=G, verbose=verbose)
+            res = self.compute_boozerls_newton(tol=1e-11, maxiter=40, constraint_weight=self.constraint_weight, iota=iota, G=G, verbose=verbose)
             return res
         else:
             raise Exception(f"boozer_type not supported: {boozer_type}")
@@ -609,7 +609,7 @@ def solve_residual_equation_exactly_newton(self, tol=1e-10, maxiter=10, iota=0.,
 
         P, L, U = lu(J)
         res = {
-            "residual": r, "jacobian": J, "iter": i, "success": norm <= tol, "G": G, "s": s, "iota": iota, "PLU": (P, L, U),
+            "residual": r, "jacobian": J, "iter": i, "success": norm <= tol, "G": G, "iota": iota, "PLU": (P, L, U),
             "mask": mask, 'type': 'exact', "vjp": boozer_surface_dexactresidual_dcoils_dcurrents_vjp
         }
 
@@ -794,7 +794,7 @@ def compute_boozerls_newton(self, tol=1e-12, maxiter=10, constraint_weight=1., i
                 x, derivatives=2, constraint_weight=constraint_weight, optimize_G=G is not None)
             norm = np.linalg.norm(dval, ord=np.inf)
 
-            if norm > 1e2:
+            if norm > 1e1:
                 break
 
             i = i+1

src/simsopt/geo/curve.py

@@ -878,4 +878,4 @@ def create_equally_spaced_curves(ncurves, nfp, stellsym, R0=1.0, R1=0.5, order=6
         curve.set("zs(1)", R1)
         curve.x = curve.x  # need to do this to transfer data to C++
         curves.append(curve)
-    return curves
+    return curves

src/simsopt/geo/surfacexyztensorfourier.py

@@ -166,4 +166,8 @@ def npsame(a, b):
         if npsame(phis, np.linspace(0, 1/(2*self.nfp), ntor+1, endpoint=False)) and \
                 npsame(thetas, np.linspace(0, 1, 2*mpol+1, endpoint=False)):
             mask[0, mpol+1:] = False
+        if npsame(phis, np.linspace(0, 1/self.nfp, 2*ntor+1, endpoint=False)[:ntor+1]) and \
+                npsame(thetas, np.linspace(0, 1, 2*mpol+1, endpoint=False)):
+            mask[0, mpol+1:] = False
+
         return mask

src/simsoptpp/python.cpp

@@ -18,6 +18,11 @@ typedef xt::pyarray<double> PyArray;
 #include "boozerradialinterpolant.h"
 #include "boozerresidual.h"
 
+#include "linking_number.hh"
+typedef LK::LinkingNumber<double> LK_class;
+
+
+
 namespace py = pybind11;
 
 using std::vector;
@@ -143,6 +148,32 @@ PYBIND11_MODULE(simsoptpp, m) {
         return 0.5 * res / (nphi*ntheta);
     });
 
+    m.def("ln", [](PyArray& A, PyArray& B) {
+            LK_class lk(2);
+            int nseg1 = A.shape(0);
+            int nseg2 = B.shape(0);
+
+            double c1[10000][3];
+            double c2[10000][3];
+            for(int i; i < nseg1; i++){
+                c1[i][0] = A(i, 0);
+                c1[i][1] = A(i, 1);
+                c1[i][2] = A(i, 2);
+            }
+
+            for(int i; i < nseg2; i++){
+                c2[i][0] = B(i, 0);
+                c2[i][1] = B(i, 1);
+                c2[i][2] = B(i, 2);
+            }
+
+            return lk.eval(c1, nseg1, c2, nseg2) ;
+        });
+
+
+
+
+
 #ifdef VERSION_INFO
     m.attr("__version__") = VERSION_INFO;
 #else