.

Makefile

@@ -30,6 +30,8 @@ CXXFLAGS += -lfftw3 -DSCTL_HAVE_FFTW
 CXXFLAGS += -lfftw3f -DSCTL_HAVE_FFTWF
 CXXFLAGS += -lfftw3l -DSCTL_HAVE_FFTWL
 
+CXXFLAGS += -I/home/malhotra/local/eigen-3.3.7
+CXXFLAGS += -I/home/malhotra/local/LBFGSpp/include
 
 RM = rm -f
 MKDIRS = mkdir -p

include/sctl/boundary_quadrature.hpp

@@ -2050,6 +2050,8 @@ template <class Real> class Quadrature {
 
 template <class Real, Integer ORDER=10> class Stellarator {
   private:
+    static constexpr Integer order_singular = 15;
+    static constexpr Integer order_direct = 35;
     static constexpr Integer COORD_DIM = 3;
     static constexpr Integer ELEM_DIM = COORD_DIM-1;
     using ElemBasis = Basis<Real, ELEM_DIM, ORDER>;
@@ -2568,7 +2570,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
       Vector<Real> x_(Nelem * Nnodes + S.Nsurf());
       x_ = 0;
       ParallelSolver<Real> linear_solver(comm, true);
-      linear_solver(&x_, BIOp, rhs_, 1e-8, 100);
+      linear_solver(&x_, BIOp, rhs_, 1e-6, 100);
 
       sigma.ReInit(Nelem);
       for (Long i = 0; i < Nelem; i++) {
@@ -2718,7 +2720,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
       Vector<Real> x(b.Dim());
       x = 0;
       ParallelSolver<Real> linear_solver(comm, true);
-      linear_solver(&x, BIOp, b, 1e-8, 100);
+      linear_solver(&x, BIOp, b, 1e-6, 100);
       return x;
     }
 
@@ -2836,6 +2838,25 @@ template <class Real, Integer ORDER=10> class Stellarator {
         elem_offset += Nelem;
       }
     }
+    static Real compute_g(const Vector<Stellarator<Real,ORDER>>& Svec, const Vector<Real>& pressure) {
+      Real g = 0;
+      compute_gvec(Svec, pressure);
+      for (Long i = 0; i < Svec.Dim(); i++) { // Set gvec
+        Vector<ElemBasis> normal, area_elem, wt(Svec[i].NElem());
+        compute_norm_area_elem(Svec[i], normal, area_elem);
+        wt = 0.5;
+        if (i == Svec.Dim()-1) {
+          Long Nsurf = Svec[i].Nsurf();
+          Long Nelem = Svec[i].NTor(Nsurf-1) * Svec[i].NPol(Nsurf-1);
+          Long offset = Svec[i].ElemDsp(Nsurf-1);
+          for (Long j = 0; j < Nelem; j++) {
+            wt[offset + j] = 1.0;
+          }
+        }
+        g += compute_inner_prod(area_elem, Svec[i].gvec, wt);
+      }
+      return g;
+    }
 
     Stellarator(const Vector<Long>& NtNp = Vector<Long>()) {
       NtNp_ = NtNp;
@@ -2887,9 +2908,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
       return elements.NElem();
     }
 
-    static Vector<ElemBasis> compute_gradient(const Stellarator<Real,ORDER>& S_, const Vector<Real>& pressure, const Vector<Real>& flux_tor_, const Vector<Real>& flux_pol_) {
-      constexpr Integer order_singular = 15;
-      constexpr Integer order_direct = 35;
+    static Vector<ElemBasis> compute_gradient(const Stellarator<Real,ORDER>& S_, const Vector<Real>& pressure, const Vector<Real>& flux_tor_, const Vector<Real>& flux_pol_, Real* g_ptr = nullptr) {
       Comm comm = Comm::World();
 
       Vector<Stellarator<Real,ORDER>> Svec(S_.Nsurf());
@@ -2950,6 +2969,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
       }
       compute_gvec(Svec, pressure);
       compute_dgdB(Svec, pressure);
+      if (g_ptr != nullptr) g_ptr[0] = compute_g(Svec, pressure);
 
       auto compute_gradient = [&comm] (const Stellarator<Real,ORDER>& S) {
         const Long Nnodes = ElemBasis::Size();
@@ -3589,9 +3609,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
       return dgdnu;
     }
 
-    static Vector<ElemBasis> compute_pressure_jump(const Stellarator<Real,ORDER>& S_, const Vector<Real>& pressure, const Vector<Real>& flux_tor_, const Vector<Real>& flux_pol_) {
-      constexpr Integer order_singular = 15;
-      constexpr Integer order_direct = 35;
+    static Vector<ElemBasis> compute_pressure_jump(const Stellarator<Real,ORDER>& S_, const Vector<Real>& pressure, const Vector<Real>& flux_tor_, const Vector<Real>& flux_pol_, Real* g_ptr = nullptr) {
       Comm comm = Comm::World();
 
       Vector<Stellarator<Real,ORDER>> Svec(S_.Nsurf());
@@ -3634,83 +3652,14 @@ template <class Real, Integer ORDER=10> class Stellarator {
         compute_invA(S.sigma, S.alpha, S.beta, S, flux_tor_[i], flux_pol_[i], comm);
         S.B = compute_B(S, S.sigma, S.alpha, S.beta);
       }
+      if (g_ptr != nullptr) g_ptr[0] = compute_g(Svec, pressure);
       return compute_pressure_jump(Svec, pressure);
     }
 
-    static Real compute_g(const Stellarator<Real,ORDER>& S_, const Vector<Real>& pressure, const Vector<Real>& flux_tor_, const Vector<Real>& flux_pol_) {
-      constexpr Integer order_singular = 15;
-      constexpr Integer order_direct = 35;
-      Comm comm = Comm::World();
-
-      Vector<Stellarator<Real,ORDER>> Svec(S_.Nsurf());
-      for (Long i = 0; i < S_.Nsurf(); i++) { // Set Svec[i] (quadratures, B)
-        const Long elem_dsp = (i==0 ? 0 : S_.ElemDsp(i-1));
-        const Long Nnodes = ElemBasis::Size();
-        Stellarator<Real,ORDER>& S = Svec[i];
-        if (i == 0) { // Init S
-          Vector<Long> NtNp;
-          NtNp.PushBack(S_.NTor(i));
-          NtNp.PushBack(S_.NPol(i));
-          S = Stellarator<Real,ORDER>(NtNp);
-        } else {
-          Vector<Long> NtNp;
-          NtNp.PushBack(S_.NTor(i-1));
-          NtNp.PushBack(S_.NPol(i-1));
-          NtNp.PushBack(S_.NTor(i));
-          NtNp.PushBack(S_.NPol(i));
-          S = Stellarator<Real,ORDER>(NtNp);
-        }
-        for (Long j = 0; j < S.NElem(); j++) { // Set S coordinates
-          for (Long k = 0; k < Nnodes; k++) {
-            S.Elem(j,0)[k] = S_.Elem(elem_dsp+j,0)[k];
-            S.Elem(j,1)[k] = S_.Elem(elem_dsp+j,1)[k];
-            S.Elem(j,2)[k] = S_.Elem(elem_dsp+j,2)[k];
-          }
-        }
-
-        SetupQuadrature(S.quadrature_BS   , S, S.BiotSavart    , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
-        SetupQuadrature(S.quadrature_FxU  , S, S.Laplace_FxU   , order_singular, order_direct, -1.0, comm);
-        SetupQuadrature(S.quadrature_FxdU , S, S.Laplace_FxdU  , order_singular, order_direct, -1.0, comm);
-
-        { // Set Bt0, Bp0, dBt0, dBp0
-          Vector<ElemBasis> Jt, Jp;
-          compute_harmonic_vector_potentials(Jt, Jp, S);
-          EvalQuadrature(S.Bt0 , S.quadrature_BS , S, Jp, S.BiotSavart);
-          EvalQuadrature(S.Bp0 , S.quadrature_BS , S, Jt, S.BiotSavart);
-        }
-
-        compute_invA(S.sigma, S.alpha, S.beta, S, flux_tor_[i], flux_pol_[i], comm);
-        S.B = compute_B(S, S.sigma, S.alpha, S.beta);
-      }
-
-      auto compute_g = [] (const Vector<Stellarator<Real,ORDER>>& Svec, const Vector<Real>& pressure) {
-        Real g = 0;
-        compute_gvec(Svec, pressure);
-        for (Long i = 0; i < Svec.Dim(); i++) { // Set gvec
-          Vector<ElemBasis> normal, area_elem, wt(Svec[i].NElem());
-          compute_norm_area_elem(Svec[i], normal, area_elem);
-          wt = 0.5;
-          if (i == Svec.Dim()-1) {
-            Long Nsurf = Svec[i].Nsurf();
-            Long Nelem = Svec[i].NTor(Nsurf-1) * Svec[i].NPol(Nsurf-1);
-            Long offset = Svec[i].ElemDsp(Nsurf-1);
-            for (Long j = 0; j < Nelem; j++) {
-              wt[offset + j] = 1.0;
-            }
-          }
-          g += compute_inner_prod(area_elem, Svec[i].gvec, wt);
-        }
-        return g;
-      };
-      return compute_g(Svec, pressure);
-    }
-
 
 
 
     static void test() {
-      constexpr Integer order_singular = 15;
-      constexpr Integer order_direct = 35;
       Comm comm = Comm::World();
       Profile::Enable(true);
 
@@ -3896,9 +3845,10 @@ template <class Real, Integer ORDER=10> class Stellarator {
               S2.Elem(i, 1) += dXdt[i*COORD_DIM+1] * 1.0 * dt;
               S2.Elem(i, 2) += dXdt[i*COORD_DIM+2] * 1.0 * dt;
             }
-            Real g0 = compute_g(S0, pressure, flux_tor, flux_pol);
-            Real g1 = compute_g(S1, pressure, flux_tor, flux_pol);
-            Real g2 = compute_g(S2, pressure, flux_tor, flux_pol);
+            Real g0, g1, g2;
+            compute_pressure_jump(S0, pressure, flux_tor, flux_pol, &g0);
+            compute_pressure_jump(S1, pressure, flux_tor, flux_pol, &g1);
+            compute_pressure_jump(S2, pressure, flux_tor, flux_pol, &g2);
             { // Calculate optimal step size dt
               Real a = 2*g0 - 4*g1 + 2*g2;
               Real b =-3*g0 + 4*g1 -   g2;
@@ -3980,16 +3930,15 @@ template <class Real, Integer ORDER=10> class Stellarator {
             S1.Elem(i, 2)[j] += 0.5 * eps * normal[i*COORD_DIM+2][j] * nu[i][j];
           }
         }
-        Real g0 = compute_g(S0, pressure, flux_tor, flux_pol);
-        Real g1 = compute_g(S1, pressure, flux_tor, flux_pol);
+        Real g0, g1;
+        compute_pressure_jump(S0, pressure, flux_tor, flux_pol, &g0);
+        compute_pressure_jump(S1, pressure, flux_tor, flux_pol, &g1);
         std::cout<<"g0 = "<<g0<<";  g1 = "<<g1<<";  dgdnu_ = "<<(g1-g0)/eps<<'\n';
         std::cout<<"dgdnu  = "<<compute_inner_prod(area_elem, dgdnu, nu)<<'\n';
       }
     }
 
     static void test_() {
-      constexpr Integer order_singular = 15;
-      constexpr Integer order_direct = 35;
       Comm comm = Comm::World();
       Profile::Enable(true);
 
@@ -4262,8 +4211,6 @@ template <class Real, Integer ORDER=10> class Stellarator {
     }
 
     static void test_askham() {
-      constexpr Integer order_singular = 15;
-      constexpr Integer order_direct = 35;
       Comm comm = Comm::World();
       Profile::Enable(true);
 
@@ -4850,7 +4797,7 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       const Long N = Nelem * COORD_DIM * Nnodes;
       SCTL_ASSERT(x.rows() == N);
 
-      auto filter = [](const Stellarator<Real,ORDER>& S, Vector<ElemBasis>& f) {
+      auto filter = [](const Stellarator<Real,ORDER>& S, const Comm& comm, Vector<ElemBasis>& f, Real sigma) {
         auto cheb2grid = [] (const Vector<ElemBasis>& X, Long Mt, Long Mp, Long Nt, Long Np) {
           const Long dof = X.Dim() / (Mt * Mp);
           SCTL_ASSERT(X.Dim() == Mt * Mp *dof);
@@ -4955,6 +4902,40 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
           }
           return X;
         };
+        auto fourier_filter = [](sctl::Vector<Real>& X, long Nt_, long Np_, Real sigma, const Comm& comm) {
+          long dof = X.Dim() / (Nt_ * Np_);
+          SCTL_ASSERT(X.Dim() == dof * Nt_ * Np_);
+
+          sctl::FFT<Real> fft_r2c, fft_c2r;
+          sctl::StaticArray<sctl::Long, 2> fft_dim = {Nt_, Np_};
+          fft_r2c.Setup(sctl::FFT_Type::R2C, 1, sctl::Vector<sctl::Long>(2, fft_dim, false), omp_get_max_threads());
+          fft_c2r.Setup(sctl::FFT_Type::C2R, 1, sctl::Vector<sctl::Long>(2, fft_dim, false), omp_get_max_threads());
+
+          long Nt = Nt_;
+          long Np = fft_r2c.Dim(1) / (Nt * 2);
+          SCTL_ASSERT(fft_r2c.Dim(1) == Nt * Np * 2);
+
+          //auto filter_fn = [](Real x2, Real sigma) {return exp(-x2/(2*sigma*sigma));};
+          auto filter_fn = [](Real x2, Real sigma) {return (x2<sigma*sigma?1.0:0.0);};
+
+          sctl::Vector<Real> normal, gradX;
+          biest::SurfaceOp<Real> op(comm, Nt_, Np_);
+          sctl::Vector<Real> coeff(fft_r2c.Dim(1));
+          for (long k = 0; k < dof; k++) {
+            sctl::Vector<Real> X_(Nt_*Np_, X.begin() + k*Nt_*Np_, false);
+            fft_r2c.Execute(X_, coeff);
+            for (long t = 0; t < Nt; t++) {
+              for (long p = 0; p < Np; p++) {
+                Real tt = (t - (t > Nt / 2 ? Nt : 0)) / (Real)(Nt / 2);
+                Real pp = p / (Real)Np;
+                Real f = filter_fn(tt*tt+pp*pp, sigma);
+                coeff[(t * Np + p) * 2 + 0] *= f;
+                coeff[(t * Np + p) * 2 + 1] *= f;
+              }
+            }
+            fft_c2r.Execute(coeff, X_);
+          }
+        };
 
         Long dof = f.Dim() / S.NElem();
         SCTL_ASSERT(f.Dim() == S.NElem() * dof);
@@ -4963,15 +4944,17 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
           const Long Mp = S.NPol(i);
           const Long Nelem = Mt * Mp;
           const Long offset = S.ElemDsp(i);
-          const Long Nt = Mt * ORDER / 10;
-          const Long Np = Mp * ORDER / 10;
+          const Long Nt = Mt * ORDER * 4;
+          const Long Np = Mp * ORDER * 4;
 
           Vector<ElemBasis> f_(Nelem*dof, f.begin() + offset*dof, false);
           Vector<Real> f_fourier = cheb2grid(f_, Mt, Mp, Nt, Np);
+          fourier_filter(f_fourier, Nt, Np, 0.25 * sigma, comm);
           f_ = grid2cheb(f_fourier, Nt, Np, Mt, Mp);
         }
       };
 
+      Real g;
       for (Long i = 0; i < Nelem; i++) { // Set S_
         for (Long j = 0; j < Nnodes; j++) {
           S_.Elem(i,0)[j] = x[(i*Nnodes+j)*COORD_DIM+0];
@@ -4979,9 +4962,8 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
           S_.Elem(i,2)[j] = x[(i*Nnodes+j)*COORD_DIM+2];
         }
       }
-      Real g = Stellarator<Real,ORDER>::compute_g(S_, pressure_, flux_tor_, flux_pol_);
-      Vector<ElemBasis> dgdnu = Stellarator<Real,ORDER>::compute_gradient(S_, pressure_, flux_tor_, flux_pol_);
-      Vector<ElemBasis> dXdt(N);
+      Vector<ElemBasis> dgdnu = Stellarator<Real,ORDER>::compute_gradient(S_, pressure_, flux_tor_, flux_pol_, &g);
+      Vector<ElemBasis> dXdt(Nelem*COORD_DIM);
       { // Set dXdt
         dXdt = 0;
         const Long Nnodes = ElemBasis::Size();
@@ -4994,7 +4976,7 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
             dXdt[i*COORD_DIM+2][j] = normal[i*COORD_DIM+2][j] * dgdnu[i][j];
           }
         }
-        filter(S_, dXdt);
+        filter(S_, comm, dXdt, 0.3333);
       }
       for (Long i = 0; i < Nelem; i++) { // Set grad
         for (Long j = 0; j < Nnodes; j++) {
@@ -5009,6 +4991,11 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
         vtu.AddElems(S_.GetElemList(), dgdnu, ORDER);
         vtu.WriteVTK("dgdnu"+std::to_string(iter), comm);
       }
+      if (1) { // Write VTU
+        VTUData vtu;
+        vtu.AddElems(S_.GetElemList(), dXdt, ORDER);
+        vtu.WriteVTK("dXdt"+std::to_string(iter), comm);
+      }
       std::cout<<"iter = "<<iter<<"    g = "<<g<<'\n';
 
       iter++;
@@ -5052,8 +5039,6 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
     }
 
     static void test() {
-      constexpr Integer order_singular = 25;
-      constexpr Integer order_direct = 35;
       Comm comm = Comm::World();
       Profile::Enable(true);
 
@@ -5062,14 +5047,12 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       Vector<Real> flux_tor(Nsurf), flux_pol(Nsurf), pressure(Nsurf);
       { // Init S, flux_tor, flux_pol, pressure
         Vector<Long> NtNp;
-        NtNp.PushBack(50);
-        NtNp.PushBack(8);
-        NtNp.PushBack(50);
-        NtNp.PushBack(8);
-        //for (Long i = 0; i < Nsurf; i++) {
-        //  NtNp.PushBack(30);
-        //  NtNp.PushBack(4);
-        //}
+        for (Long i = 0; i < Nsurf; i++) {
+          //NtNp.PushBack(50);
+          //NtNp.PushBack(8);
+          NtNp.PushBack(30);
+          NtNp.PushBack(4);
+        }
         S = Stellarator<Real,ORDER>(NtNp);
         flux_tor = 1;
         flux_pol = 1;