.

include/sctl/boundary_quadrature.hpp

@@ -2060,8 +2060,8 @@ template <class Real> class Quadrature {
 
 template <class Real, Integer ORDER=10> class Stellarator {
   private:
-    static constexpr Integer order_singular = 25;
-    static constexpr Integer order_direct = 35;
+    static constexpr Integer order_singular = 20;
+    static constexpr Integer order_direct = 25;
     static constexpr Integer COORD_DIM = 3;
     static constexpr Integer ELEM_DIM = COORD_DIM-1;
     using ElemBasis = Basis<Real, ELEM_DIM, ORDER>;
@@ -2231,7 +2231,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
     }
     static void compute_harmonic_vector_potentials(Vector<ElemBasis>& Jt, Vector<ElemBasis>& Jp, const Stellarator<Real,ORDER>& S) {
       Comm comm = Comm::World();
-      Real gmres_tol = 1e-8;
+      Real gmres_tol = 1e-10;
       Long max_iter = 100;
 
       auto cheb2grid = [] (const Vector<ElemBasis>& X, Long Mt, Long Mp, Long Nt, Long Np) {
@@ -2343,7 +2343,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
       if (Jp.Dim() != Nelem * COORD_DIM) Jp.ReInit(Nelem * COORD_DIM);
       if (Jt.Dim() != Nelem * COORD_DIM) Jt.ReInit(Nelem * COORD_DIM);
       for (Long k = 0; k < S.Nsurf(); k++) {
-        Long Nt = S.NTor(k)*ORDER, Np = S.NPol(k)*ORDER;
+        Long Nt = S.NTor(k)*ORDER*2, Np = S.NPol(k)*ORDER*2;
         const auto& X_ = S.GetElemList().ElemVector();
         Vector<ElemBasis> X(S.NTor(k)*S.NPol(k)*COORD_DIM, (Iterator<ElemBasis>)X_.begin()+S.ElemDsp(k)*COORD_DIM, false);
 
@@ -2584,7 +2584,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-6, 100);
+      linear_solver(&x_, BIOp, rhs_, 1e-10, 100);
 
       sigma.ReInit(Nelem);
       for (Long i = 0; i < Nelem; i++) {
@@ -2739,7 +2739,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-6, 100);
+      linear_solver(&x, BIOp, b, 1e-8, 100);
       return x;
     }
 
@@ -3657,7 +3657,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           }
         }
 
-        SetupQuadrature(S.quadrature_BS   , S, S.BiotSavart    , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
+        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);
 
@@ -3666,6 +3666,38 @@ template <class Real, Integer ORDER=10> class Stellarator {
           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);
+
+          Vector<ElemBasis> normal, area_elem;
+          compute_norm_area_elem(S, normal, area_elem);
+          if (S.Nsurf() == 2) {
+            Long Nelem0 = S.NTor(0)*S.NPol(0);
+            for (Long i = 0; i < Nelem0*COORD_DIM; i++) {
+              for (Long j = 0; j < Nnodes; j++) {
+                normal[i][j] *= -1.0;
+              }
+            }
+          }
+
+          const Long Nelem = S.NElem();
+          const Long Nnodes = ElemBasis::Size();
+          for (Long j = 0; j < Nelem; j++) {
+            for (Long k = 0; k < Nnodes; k++) {
+              Real Jxn[COORD_DIM];
+              Jxn[0] = Jp[j*COORD_DIM+1][k] * normal[j*COORD_DIM+2][k] - Jp[j*COORD_DIM+2][k] * normal[j*COORD_DIM+1][k];
+              Jxn[1] = Jp[j*COORD_DIM+2][k] * normal[j*COORD_DIM+0][k] - Jp[j*COORD_DIM+0][k] * normal[j*COORD_DIM+2][k];
+              Jxn[2] = Jp[j*COORD_DIM+0][k] * normal[j*COORD_DIM+1][k] - Jp[j*COORD_DIM+1][k] * normal[j*COORD_DIM+0][k];
+              S.Bt0[j*COORD_DIM+0][k] += 0.5 * Jxn[0];
+              S.Bt0[j*COORD_DIM+1][k] += 0.5 * Jxn[1];
+              S.Bt0[j*COORD_DIM+2][k] += 0.5 * Jxn[2];
+
+              Jxn[0] = Jt[j*COORD_DIM+1][k] * normal[j*COORD_DIM+2][k] - Jt[j*COORD_DIM+2][k] * normal[j*COORD_DIM+1][k];
+              Jxn[1] = Jt[j*COORD_DIM+2][k] * normal[j*COORD_DIM+0][k] - Jt[j*COORD_DIM+0][k] * normal[j*COORD_DIM+2][k];
+              Jxn[2] = Jt[j*COORD_DIM+0][k] * normal[j*COORD_DIM+1][k] - Jt[j*COORD_DIM+1][k] * normal[j*COORD_DIM+0][k];
+              S.Bp0[j*COORD_DIM+0][k] += 0.5 * Jxn[0];
+              S.Bp0[j*COORD_DIM+1][k] += 0.5 * Jxn[1];
+              S.Bp0[j*COORD_DIM+2][k] += 0.5 * Jxn[2];
+            }
+          }
         }
 
         compute_invA(S.sigma, S.alpha, S.beta, S, flux_tor_[i], flux_pol_[i], comm);
@@ -5581,7 +5613,9 @@ template <class Real, Integer ORDER=10> class Stellarator {
 };
 
 template <class Real, Integer ORDER=10> class MHDEquilib {
-    static constexpr Integer fourier_upsample = 4;
+    static constexpr Integer fourier_dim0 = 50*10*4;
+    static constexpr Integer fourier_dim1 = 10*10*4;
+    //static constexpr Integer fourier_upsample = 4;
     static constexpr Integer COORD_DIM = 3;
     static constexpr Integer ELEM_DIM = COORD_DIM-1;
     using ElemBasis = Basis<Real, ELEM_DIM, ORDER>;
@@ -5712,8 +5746,8 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
         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 tt = (t - (t > Nt / 2 ? Nt : 0)) * (2*Np/(Real)Nt);
+            Real pp = p;
             Real f = filter_fn(tt*tt+pp*pp, sigma);
             coeff[(t * Np + p) * 2 + 0] *= f;
             coeff[(t * Np + p) * 2 + 1] *= f;
@@ -5724,39 +5758,101 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
     };
 
     static void fourier_print_spectrum(const Vector<Real>& X, long Nt_, long Np_) {
-      long dof = X.Dim() / (Nt_ * Np_);
-      SCTL_ASSERT(X.Dim() == dof * Nt_ * Np_);
-
-      FFT<Real> fft_r2c, fft_c2r;
-      StaticArray<Long, 2> fft_dim = {Nt_, Np_};
-      fft_r2c.Setup(FFT_Type::R2C, 1, Vector<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);
-
-      Vector<Real> coeff(fft_r2c.Dim(1));
-      Vector<Real> spectrum(20); spectrum = 0;
-      for (long k = 0; k < dof; k++) {
-        const Vector<Real> X_(Nt_*Np_, (Iterator<Real>)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;
-            Long freq = (Long)(sqrt<Real>(tt*tt+pp*pp)*spectrum.Dim());
-            if (freq<20) spectrum[freq] += coeff[(t*Np+p)*2+0]*coeff[(t*Np+p)*2+0];
-            if (freq<20) spectrum[freq] += coeff[(t*Np+p)*2+1]*coeff[(t*Np+p)*2+1];
-          }
+    //  long dof = X.Dim() / (Nt_ * Np_);
+    //  SCTL_ASSERT(X.Dim() == dof * Nt_ * Np_);
+
+    //  FFT<Real> fft_r2c, fft_c2r;
+    //  StaticArray<Long, 2> fft_dim = {Nt_, Np_};
+    //  fft_r2c.Setup(FFT_Type::R2C, 1, Vector<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);
+
+    //  Vector<Real> coeff(fft_r2c.Dim(1));
+    //  Vector<Real> spectrum(200); spectrum = 0;
+    //  for (long k = 0; k < dof; k++) {
+    //    const Vector<Real> X_(Nt_*Np_, (Iterator<Real>)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;
+    //        Long freq = (Long)(sqrt<Real>(tt*tt+pp*pp)*spectrum.Dim());
+    //        if (freq<spectrum.Dim()) spectrum[freq] += coeff[(t*Np+p)*2+0]*coeff[(t*Np+p)*2+0];
+    //        if (freq<spectrum.Dim()) spectrum[freq] += coeff[(t*Np+p)*2+1]*coeff[(t*Np+p)*2+1];
+    //      }
+    //    }
+    //  }
+    //  for (Long i = 0; i < spectrum.Dim(); i++) {
+    //    spectrum[i] = log<Real>(sqrt<Real>(spectrum[i]))/log<Real>(10.0);
+    //  }
+    //  std::cout<<"spectrum = "<<spectrum<<'\n';
+    };
+    static void fourier_print_spectrum(std::string var_name, const sctl::Vector<Real>& B0, sctl::Long Nt0, sctl::Long Np0) {
+      auto Resample = [](const sctl::Vector<Real>& B, long Nt, long Np, long Nt0, long Np0) {
+        sctl::Vector<Real> B0;
+        biest::SurfaceOp<Real>::Upsample(B,Nt,Np, B0,Nt0,Np0);
+        return B0;
+      };
+      auto max_rel_err = [](const sctl::Vector<Real> A, const sctl::Vector<Real> B) {
+        SCTL_ASSERT(A.Dim() == B.Dim());
+        Real max_err = 0;
+        Real max_val = 1e-20;
+        for (sctl::Long i = 0; i < A.Dim(); i++) {
+          max_err = std::max<Real>(max_err, fabs(A[i]-B[i]));
+          max_val = std::max<Real>(max_val, fabs(A[i]));
         }
-      }
-      for (Long i = 0; i < spectrum.Dim(); i++) {
-        spectrum[i] = log<Real>(sqrt<Real>(spectrum[i]))/log<Real>(10.0);
-      }
-      std::cout<<"spectrum = "<<spectrum<<'\n';
+        return max_err/max_val;
+      };
+      auto max_err = [](const sctl::Vector<Real> A, const sctl::Vector<Real> B) {
+        SCTL_ASSERT(A.Dim() == B.Dim());
+        Real max_err = 0;
+        for (sctl::Long i = 0; i < A.Dim(); i++) {
+          max_err = std::max<Real>(max_err, fabs(A[i]-B[i]));
+        }
+        return max_err;
+      };
+
+      std::cout<<var_name<<"=[";
+      for (sctl::Long i = 1; i < 140; i++) {
+        sctl::Long Nt1 = 6*i, Np1 = i;
+        auto B1 = Resample(Resample(B0, Nt0,Np0, Nt1,Np1), Nt1,Np1, Nt0,Np0);
+        std::cout<<log(max_err(B0, B1))/log(10)<<' ';
+      }
+      std::cout<<"];\n";
+
+      //if (0) {
+      //  auto B1 = Resample(B0, Nt0,Np0, 70*30,14*30);
+      //  B1.Write(var_name.c_str());
+      //} else {
+      //  auto B1 = Resample(B0, Nt0,Np0, 70*30,14*30);
+      //  sctl::Vector<Real> B2; B2.Read(var_name.c_str());
+
+      //  Real max_err = 0, max_val = 0;
+      //  for (sctl::Long i = 0; i < B1.Dim(); i++) {
+      //    max_err = std::max<Real>(max_err, fabs(B1[i]-B2[i]));
+      //    max_val = std::max<Real>(max_val, fabs(B2[i]));
+      //  }
+      //  std::cout<<"Error "<<var_name<<" = "<<max_err/max_val<<'\n';
+      //}
     };
+    static void fourier_print_spectrum(std::string var_name, const Vector<Real>& X_, Stellarator<Real,ORDER> S_) {
+      Long dof = X_.Dim()/(S_.Nsurf()*fourier_dim0*fourier_dim1);
+      SCTL_ASSERT(dof * (S_.Nsurf()*fourier_dim0*fourier_dim1) == X_.Dim());
+      for (Long i = 0; i < S_.Nsurf(); i++) {
+        const Long Mt = S_.NTor(i);
+        const Long Mp = S_.NPol(i);
+        const Long offset = S_.ElemDsp(i);
+        const Long Nt = fourier_dim0;
+        const Long Np = fourier_dim1;
+        const Vector<Real> X(dof*Nt*Np, (Iterator<Real>)X_.begin() + dof * i * (fourier_dim0*fourier_dim1), false);
+        fourier_print_spectrum(var_name+std::to_string(i),X,Nt,Np);
+      }
+      std::cout<<"\n";
+    }
 
-    static void filter(const Stellarator<Real,ORDER>& S, const Comm& comm, Vector<ElemBasis>& f, Real sigma) {
+    static void filter_deprecated(const Stellarator<Real,ORDER>& S, const Comm& comm, Vector<ElemBasis>& f, Real sigma) {
       Long dof = f.Dim() / S.NElem();
       SCTL_ASSERT(f.Dim() == S.NElem() * dof);
       for (Long i = 0; i < S.Nsurf()-1; i++) {
@@ -5780,16 +5876,16 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       const Long dof = f.Dim() / Nelem;
       SCTL_ASSERT(Nelem * dof == f.Dim());
 
-      Vector<Real> f_fourier(dof * Nelem * (ORDER*ORDER*fourier_upsample*fourier_upsample));
+      Vector<Real> f_fourier(dof * S.Nsurf() * (fourier_dim0*fourier_dim1));
       for (Long i = 0; i < S.Nsurf(); i++) {
         const Long Mt = S.NTor(i);
         const Long Mp = S.NPol(i);
         const Long offset = S.ElemDsp(i);
-        const Long Nt = Mt * ORDER * fourier_upsample;
-        const Long Np = Mp * ORDER * fourier_upsample;
+        const Long Nt = fourier_dim0;
+        const Long Np = fourier_dim1;
 
         const Vector<ElemBasis> f_(Mt*Mp*dof, (Iterator<ElemBasis>)f.begin() + offset*dof, false);
-        Vector<Real> f_fourier_(dof*Nt*Np, f_fourier.begin() + dof*offset * (ORDER*ORDER*fourier_upsample*fourier_upsample), false);
+        Vector<Real> f_fourier_(dof*Nt*Np, f_fourier.begin() + dof*i * (fourier_dim0*fourier_dim1), false);
         f_fourier_ = cheb2grid(f_, Mt, Mp, Nt, Np);
         SCTL_ASSERT(f_fourier_.Dim() == dof*Nt*Np);
       }
@@ -5798,26 +5894,26 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
     static Vector<ElemBasis> grid2cheb(const Stellarator<Real,ORDER>& S, const Vector<Real>& f_fourier) {
       const Long Nnodes = ElemBasis::Size();
       const Long Nelem = S.NElem();
-      const Long dof = f_fourier.Dim() / (Nelem * (ORDER*ORDER*fourier_upsample*fourier_upsample));
-      SCTL_ASSERT(dof * Nelem * (ORDER*ORDER*fourier_upsample*fourier_upsample) == f_fourier.Dim());
+      const Long dof = f_fourier.Dim() / (S.Nsurf() * (fourier_dim0*fourier_dim1));
+      SCTL_ASSERT(dof * (S.Nsurf() * (fourier_dim0*fourier_dim1)) == f_fourier.Dim());
 
       Vector<ElemBasis> f(Nelem * dof);
       for (Long i = 0; i < S.Nsurf(); i++) {
         const Long Mt = S.NTor(i);
         const Long Mp = S.NPol(i);
         const Long offset = S.ElemDsp(i);
-        const Long Nt = Mt * ORDER * fourier_upsample;
-        const Long Np = Mp * ORDER * fourier_upsample;
+        const Long Nt = fourier_dim0;
+        const Long Np = fourier_dim1;
 
         Vector<ElemBasis> f_(Mt*Mp*dof, f.begin() + offset*dof, false);
-        const Vector<Real> f_fourier_(dof*Nt*Np, (Iterator<Real>)f_fourier.begin() + dof*offset * (ORDER*ORDER*fourier_upsample*fourier_upsample), false);
+        const Vector<Real> f_fourier_(dof*Nt*Np, (Iterator<Real>)f_fourier.begin() + dof*i * (fourier_dim0*fourier_dim1), false);
         f_ = grid2cheb(f_fourier_, Nt, Np, Mt, Mp);
         SCTL_ASSERT(f_.Dim() == Mt*Mp*dof);
       }
       return f;
     }
 
-    template <class Real, class GradOp> static Long GradientDescent3(GradOp& grad_op, Eigen::VectorXd& x, Real& fx, Long max_iter, Real tol) {
+    template <class Real, class GradOp> static Long GradientDescent3(GradOp& grad_fn, Eigen::VectorXd& x, Real& fx, Long max_iter, Real tol) {
       Real dt = 0.1;
       Real step_tol = 1e-1;
       for (Long iter = 0; iter < max_iter; iter++) {
@@ -5839,7 +5935,7 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
 
         Real error;
         Eigen::VectorXd x_(x.size());
-        time_step(&x_, dt, x, grad_op, &error);
+        time_step(&x_, dt, x, grad_fn, &error);
         if (error > 2.0 * step_tol) {
           dt *= 0.5;
         } else {
@@ -5851,7 +5947,7 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       return max_iter;
     }
 
-    template <class Real, class GradOp> static Long GradientDescent2(GradOp& grad_op, Eigen::VectorXd& x, Real& fx, Long max_iter, Real tol) {
+    template <class Real, class GradOp> static Long GradientDescent2(GradOp& grad_fn, Eigen::VectorXd& x, Real& fx, Long max_iter, Real tol) {
       Real dt = 0.1;
       Real step_tol = 1e-1;
       for (Long iter = 0; iter < max_iter; iter++) {
@@ -5878,7 +5974,7 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
 
         Real error;
         Eigen::VectorXd x_(x.size());
-        time_step(&x_, dt, x, grad_op, &error);
+        time_step(&x_, dt, x, grad_fn, &error);
         if (error > 2.0 * step_tol) {
           dt *= 0.5;
         } else {
@@ -5890,7 +5986,7 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       return max_iter;
     }
 
-    template <class Real, class GradOp> static Long GradientDescent2_(GradOp& grad_op, Eigen::VectorXd& x, Real& fx, Long max_iter, Real tol) {
+    template <class Real, class GradOp> static Long GradientDescent2_(GradOp& grad_fn, Eigen::VectorXd& x, Real& fx, Long max_iter, Real tol) {
       Real dt_ = 0, fx_ = 0;
       Eigen::VectorXd x_ = x;
 
@@ -5898,7 +5994,7 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       for (Long iter = 0; iter < max_iter; iter++) {
         Eigen::VectorXd grad0(x.size()), grad1(x.size()), grad2(x.size());
         fx_ = fx;
-        fx = grad_op(x, grad0);
+        fx = grad_fn(x, grad0);
         if (fx < tol) return iter;
         if (iter && fx > fx_) {
           x = x_;
@@ -5911,10 +6007,10 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
         }
         { // Update dt
           Real fx1, fx2;
-          fx1 = grad_op(x - grad0 * dt * 0.50, grad1);
-          fx1 = grad_op(x - grad0 * dt * 0.25 - grad1 * dt * 0.25, grad1);
-          fx2 = grad_op(x - grad1 * dt, grad2);
-          fx2 = grad_op(x - grad0 * dt * (1.0/6.0) - grad1 * dt * (2.0/3.0) - grad2 * dt * (1.0/6.0), grad2);
+          fx1 = grad_fn(x - grad0 * dt * 0.50, grad1);
+          fx1 = grad_fn(x - grad0 * dt * 0.25 - grad1 * dt * 0.25, grad1);
+          fx2 = grad_fn(x - grad1 * dt, grad2);
+          fx2 = grad_fn(x - grad0 * dt * (1.0/6.0) - grad1 * dt * (2.0/3.0) - grad2 * dt * (1.0/6.0), grad2);
 
           Real s;
           { // Calculate optimal step size dt
@@ -5937,17 +6033,17 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       return max_iter;
     }
 
-    template <class Real, class GradOp> static Long GradientDescent(GradOp& grad_op, Eigen::VectorXd& x, Real& fx, Long max_iter, Real tol) {
-      Real dt = 0.1;
+    template <class Real, class GradOp> static Long GradientDescent(GradOp& grad_fn, Eigen::VectorXd& x, Real& fx, Long max_iter, Real tol) {
+      Real dt = 0.192081;
       for (Long iter = 0; iter < max_iter; iter++) {
         Eigen::VectorXd grad(x.size());
-        fx = grad_op(x, grad);
+        fx = grad_fn(x, grad);
         { // Update dt
           Eigen::VectorXd grad_(x.size());
           Eigen::VectorXd x1 = x - grad * dt * 0.5;
           Eigen::VectorXd x2 = x - grad * dt * 1.0;
-          Real fx1 = grad_op(x1, grad_);
-          Real fx2 = grad_op(x2, grad_);
+          Real fx1 = grad_fn(x1, grad_);
+          Real fx2 = grad_fn(x2, grad_);
 
           { // Calculate optimal step size dt
             Real a = 2*fx - 4*fx1 + 2*fx2;
@@ -5970,13 +6066,298 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       return max_iter;
     }
 
+    template <class ValueType> static ValueType QuadraticInterp(Real t, const ValueType& v0, const ValueType& v1, const ValueType& v2, Real t0, Real t1, Real t2) {
+      ValueType v = v0 * (((t-t1)*(t-t2))/((t0-t1)*(t0-t2)));
+      v += v1 * (((t-t0)*(t-t2))/((t1-t0)*(t1-t2)));
+      v += v2 * (((t-t0)*(t-t1))/((t2-t0)*(t2-t1)));
+      return v;
+    }
+
+    template <class Real, class GradOp> static Long GradientDescentNew(GradOp& grad_fn, Eigen::VectorXd& x0, Real& fx, Long max_iter, Real tol) {
+      auto compute_inner_prod = [](const Eigen::VectorXd& A, const Eigen::VectorXd& B) {
+        Real sum = 0;
+        Long N = A.size();
+        SCTL_ASSERT(B.size() == N);
+        for (Long i = 0; i < N; i++) sum += A[i] * B[i];
+        return sum;
+      };
+      auto compute_dt_scale = [&compute_inner_prod](const Eigen::VectorXd& x, const Eigen::VectorXd& y){
+        Real dot_prod = compute_inner_prod(x, y) / sqrt<Real>(compute_inner_prod(x,x)*compute_inner_prod(y,y));
+        return 0.05 / sqrt<Real>(1-dot_prod*dot_prod);
+      };
+
+      Eigen::VectorXd grad, x = x0;
+      Real g = grad_fn(x, grad);
+
+      Real t = 0, dt = -0.1;
+      for (Long j = 0; j < max_iter; j++) {
+        Eigen::VectorXd grad0, grad1, grad2;
+        Real g0 = grad_fn(x + grad*(dt*0.5) + grad *(dt*0.5), grad0);
+        Real g1 = grad_fn(x + grad*(dt*0.5) + grad0*(dt*0.5), grad1);
+        Real g2 = grad_fn(x + grad*(dt*0.5) + grad1*(dt*0.5), grad2);
+
+
+        Eigen::VectorXd grad4, grad3;
+        Real c = compute_inner_prod(grad2-grad1,grad2-grad1) / compute_inner_prod(grad1-grad0,grad2-grad1);
+        grad3 = (grad2 - grad1*c) * (1/(1-c));
+        Real g3 = grad_fn(x + grad*(dt*0.5) + grad3*(dt*0.5), grad4);
+
+
+
+        std::cout<<sqrt<Real>(compute_inner_prod(grad -grad4,grad -grad4)/compute_inner_prod(grad4,grad4))<<' ';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad -grad3,grad -grad3)/compute_inner_prod(grad4,grad4))<<' ';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad -grad2,grad -grad2)/compute_inner_prod(grad4,grad4))<<' ';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad -grad1,grad -grad1)/compute_inner_prod(grad4,grad4))<<' ';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad -grad0,grad -grad0)/compute_inner_prod(grad4,grad4))<<'\n';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad0-grad4,grad0-grad4)/compute_inner_prod(grad4,grad4))<<' ';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad0-grad3,grad0-grad3)/compute_inner_prod(grad4,grad4))<<' ';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad0-grad2,grad0-grad2)/compute_inner_prod(grad4,grad4))<<' ';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad0-grad1,grad0-grad1)/compute_inner_prod(grad4,grad4))<<'\n';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad1-grad4,grad1-grad4)/compute_inner_prod(grad4,grad4))<<' ';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad1-grad3,grad1-grad3)/compute_inner_prod(grad4,grad4))<<' ';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad1-grad2,grad1-grad2)/compute_inner_prod(grad4,grad4))<<'\n';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad2-grad4,grad2-grad4)/compute_inner_prod(grad4,grad4))<<' ';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad2-grad3,grad2-grad3)/compute_inner_prod(grad4,grad4))<<'\n';
+        std::cout<<sqrt<Real>(compute_inner_prod(grad3-grad4,grad3-grad4)/compute_inner_prod(grad4,grad4))<<'\n';
+
+
+
+
+        Real s0 = sqrt<Real>(compute_inner_prod(grad -grad4,grad -grad4)/compute_inner_prod(grad4,grad4));
+        Real s1 = sqrt<Real>(compute_inner_prod(grad0-grad1,grad0-grad1)/compute_inner_prod(grad4,grad4));
+        Real s2 = sqrt<Real>(compute_inner_prod(grad1-grad2,grad1-grad2)/compute_inner_prod(grad4,grad4));
+        Real s3 = sqrt<Real>(compute_inner_prod(grad3-grad4,grad3-grad4)/compute_inner_prod(grad4,grad4));
+        Real dt_scale = (0.1/s0); // 0.5*(s1/s2);
+        if (s3 > s2) dt_scale = std::min(0.5, dt_scale);
+
+        std::cout<<s0<<' '<<s1<<' '<<s2<<' '<<s3<<"        dt_scale = "<<dt_scale<<'\n'; ////////////////////
+
+        if (dt_scale>0.5) {
+          t += dt;
+          g = g3;
+          x = x + grad *(dt*0.5) + grad3*(dt*0.5);
+          grad = grad4;
+          std::cout<<"iter = "<<j<<"   t = "<<t<<"   g = "<<g<<"   dt = "<<dt<<'\n';
+        } else {
+          j--;
+        }
+        dt *= dt_scale;
+      }
+      return 0;
+    }
+
+    template <class Real, class GradOp> static Long GradientDescentNew__(GradOp& grad_fn, Eigen::VectorXd& x, Real& fx, Long max_iter, Real tol) {
+      auto compute_inner_prod = [](const Eigen::VectorXd& A, const Eigen::VectorXd& B) {
+        Real sum = 0;
+        Long N = A.size();
+        SCTL_ASSERT(B.size() == N);
+        for (Long i = 0; i < N; i++) sum += A[i] * B[i];
+        return sum;
+      };
+      auto compute_dt_scale = [&compute_inner_prod](const Eigen::VectorXd& x, const Eigen::VectorXd& y){
+        Real dot_prod = compute_inner_prod(x, y) / sqrt<Real>(compute_inner_prod(x,x)*compute_inner_prod(y,y));
+        return 0.05 / sqrt<Real>(1-dot_prod*dot_prod);
+      };
+
+      Eigen::VectorXd x0, grad0, grad_op0;
+      Eigen::VectorXd x1, grad1, grad_op1;
+      Eigen::VectorXd x2, grad2, grad_op2;
+
+      x0 = x;
+      Real t = 0, dt = -0.1;
+      for (Long j = 0; j < max_iter; j++) {
+        Real g0 = grad_fn(x0, grad0);
+        Real g1 = grad_fn(x0 + grad0*(dt*0.1), grad1);
+        Real g2 = grad_fn(x0 + grad0*(dt*0.1) + grad1*(dt*0.1), grad2);
+
+        // Fit v = v0 + v1 exp(-alpha * dt)
+        auto v1 = (grad1-grad0) * (1/(dt*0.1));
+        Real alpha = (sqrt<Real>(compute_inner_prod(grad1-grad0,grad1-grad0))/sqrt<Real>(compute_inner_prod(grad2-grad1,grad2-grad1))-1) / (0.1*dt);
+        auto v0 = grad0 - v1;
+
+        x0 = x0 + v0*dt - v1 * ((exp<Real>(-alpha*dt)-1.0)/alpha);
+        t += dt;
+
+        std::cout<<"iter = "<<j<<"   t = "<<t<<"   g = "<<g0<<"   dt = "<<dt<<'\n';
+      }
+      return 0;
+    }
+
+    template <class Real, class GradOp> static Long GradientDescentNew_(GradOp& grad_fn, Eigen::VectorXd& x0, Real& fx, Long max_iter, Real tol) {
+      constexpr Long order = 3;
+      Eigen::VectorXd x[order], grad[order], grad_op[order];
+      Real g[order], t[order];
+
+      auto compute_inner_prod = [](const Eigen::VectorXd& A, const Eigen::VectorXd& B) {
+        Real sum = 0;
+        Long N = A.size();
+        SCTL_ASSERT(B.size() == N);
+        for (Long i = 0; i < N; i++) sum += A[i] * B[i];
+        return sum;
+      };
+      auto compute_dt_scale = [&compute_inner_prod](const Eigen::VectorXd& x, const Eigen::VectorXd& y){
+        Real dot_prod = compute_inner_prod(x, y) / sqrt<Real>(compute_inner_prod(x,x)*compute_inner_prod(y,y));
+        return 0.05 / sqrt<Real>(1-dot_prod*dot_prod);
+      };
+
+      t[1] = 0;
+      x[1] = x0;
+      g[1] = grad_fn(x[1], grad[1], &grad_op[1]);
+
+      Real dt = -0.1;
+      for (Long j = 1; j < order-1; j++) {
+        //t[0] = t[1] + dt;
+        //x[0] = x[1] + grad[1] * dt;
+        //g[0] = grad_fn(x[0], grad[0], &grad_op[0]);
+
+        //Real dot_prod = compute_inner_prod(grad[0], grad[1]) / sqrt<Real>(compute_inner_prod(grad[0], grad[0])*compute_inner_prod(grad[1], grad[1]));
+        //Real dt_scale = 0.1 / sqrt<Real>(1-dot_prod*dot_prod);
+        //dt *= dt_scale;
+
+        //if (dt_scale < 0.5 || g[0] >= g[1]) {
+        //  j--;
+        //  continue;
+        //}
+
+        //std::cout<<"iter = "<<j<<"   t = "<<t[0]<<"   g = "<<g[0]<<"   dt = "<<dt<<'\n';
+        //for (Long i = order-1; i >= 1; i--) {
+        //  x[i] = x[i-1];
+        //  grad[i] = grad[i-1];
+        //  grad_op[i] = grad_op[i-1];
+        //  g[i] = g[i-1];
+        //  t[i] = t[i-1];
+        //}
+      }
+
+      for (Long j = 0; j < max_iter; j++) {
+        Eigen::VectorXd x_, grad_, grad_op_;
+        Real g_ = grad_fn(x[1] + grad[1]*(dt*0.5), grad_, &grad_op_);
+        Real dt_scale = compute_dt_scale(grad_, grad[1]);
+        dt *= dt_scale;
+        if (dt_scale < 0.5 || g[0] >= g[1]) {
+          std::cout<<"dt = "<<dt<<'\n';
+          j--;
+          continue;
+        }
+
+        t[0] = t[1] + dt;
+        x[0] = x[1] + grad_*dt;
+        g[0] = grad_fn(x[0], grad[0], &grad_op[0]);
+
+        dt_scale = compute_dt_scale(grad[0], grad_);
+        dt *= dt_scale;
+        if (dt_scale < 0.5 || g[0] >= g[1]) {
+          std::cout<<"dt =  "<<dt<<'\n';
+          j--;
+          continue;
+        }
+
+        std::cout<<"iter = "<<j<<"   t = "<<t[0]<<"   g = "<<g[0]<<"   dt = "<<dt<<'\n';
+        for (Long i = order-1; i >= 1; i--) {
+          x[i] = x[i-1];
+          grad[i] = grad[i-1];
+          grad_op[i] = grad_op[i-1];
+          g[i] = g[i-1];
+          t[i] = t[i-1];
+        }
+      }
+
+
+      for (Long j = 0; j < max_iter; j++) {
+        Eigen::VectorXd x_, grad_, grad_op_;
+        Real g_ = grad_fn(x[1] + grad[1]*(dt*0.5), grad_, &grad_op_);
+        if (compute_dt_scale(grad_, grad[1])<0.5 ||  g_ > g[1]) {
+          dt = dt * compute_dt_scale(grad_, grad[1]);
+          std::cout<<"dt = "<<dt<<'\n';
+          j--;
+          continue;
+        }
+
+        t[0] = t[1] + dt;
+        x[0] = x[1] + grad_*dt;
+        g[0] = grad_fn(x[0], grad[0], &grad_op[0]);
+
+        Real dt_scale = compute_dt_scale(grad[0], grad[1]); /// todo use grad_ instead of grad[1]
+        dt *= dt_scale;
+        if (dt_scale < 0.5 || g[0] >= g[1]) {
+          std::cout<<"dt =  "<<dt<<'\n';
+          j--;
+          continue;
+        }
+
+        std::cout<<"iter = "<<j<<"   t = "<<t[0]<<"   g = "<<g[0]<<"   dt = "<<dt<<'\n';
+        for (Long i = order-1; i >= 1; i--) {
+          x[i] = x[i-1];
+          grad[i] = grad[i-1];
+          grad_op[i] = grad_op[i-1];
+          g[i] = g[i-1];
+          t[i] = t[i-1];
+        }
+      }
+
+      //for (Long iter = 0; iter < max_iter; iter++) {
+      //  t[0] = t[1] + dt;
+      //  x[0] = x[1] + grad[1] * dt;
+      //  g[0] = grad_fn(x[0], grad[0], &grad_op[0]);
+
+      //  Real dot_prod = compute_inner_prod(grad[0], grad[1]) / sqrt<Real>(compute_inner_prod(grad[0], grad[0])*compute_inner_prod(grad[1], grad[1]));
+      //  if (dot_prod < 0.95 || g[0] < g[1]) {
+      //    dt = dt * 1.5;
+      //  }
+      //  if (dot_prod < 0.85 || g[0] >= g[1]) {
+      //    dt = dt * 0.5;
+      //    continue;
+      //  }
+
+      //  for (Long i = order-2; i >= 0; i--) {
+      //    x[i+1] = x[i];
+      //    grad[i+1] = grad[i];
+      //    grad_op[i+1] = grad_op[i];
+      //    g[i+1] = g[i];
+      //    t[i+1] = t[i];
+      //  }
+      //}
+
+
+      //for (Long iter = 0; iter < max_iter; iter++) {
+      //  fx = grad_fn(x, grad, &grad_op);
+      //  { // Update dt
+      //    Eigen::VectorXd grad_(x.size());
+      //    Eigen::VectorXd x1 = x - grad * dt * 0.5;
+      //    Eigen::VectorXd x2 = x - grad * dt * 1.0;
+      //    Real fx1 = grad_fn(x1, grad_);
+      //    Real fx2 = grad_fn(x2, grad_);
+
+      //    { // Calculate optimal step size dt
+      //      Real a = 2*fx - 4*fx1 + 2*fx2;
+      //      Real b =-3*fx + 4*fx1 -   fx2;
+      //      Real c = fx;
+      //      Real s = -b/(2*a);
+      //      dt *= s;
+
+      //      Real fx_ = a*s*s + b*s + c;
+      //      std::cout<<"g = "<<fx_<<' ';
+      //      std::cout<<fx<<' ';
+      //      std::cout<<fx1<<' ';
+      //      std::cout<<fx2<<' ';
+      //      std::cout<<dt<<'\n';
+      //    }
+      //  }
+      //  x = x - grad * dt;
+      //  if (fx < tol) return iter;
+      //}
+      //return max_iter;
+      return 0;
+    }
+
+
   public:
     MHDEquilib(const Stellarator<Real,ORDER>& S, const Vector<Real>& pressure, const Vector<Real>& flux_tor, const Vector<Real>& flux_pol) {
       S_ = S;
       pressure_ = pressure;
       flux_tor_ = flux_tor;
       flux_pol_ = flux_pol;
-      iter = 0;
+      iter = 57;
     }
 
     Real operator()(const Eigen::VectorXd& x, Eigen::VectorXd& grad_direction, Eigen::VectorXd* grad_op_ptr = nullptr) {
@@ -6047,29 +6428,30 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       Stellarator<Real,ORDER>::compute_norm_area_elem(S_, normal, area_elem);
 
       Real g;
-      Vector<ElemBasis> dgdnu = Stellarator<Real,ORDER>::compute_gradient(S_, pressure_, flux_tor_, flux_pol_, &g);
-      //Vector<ElemBasis> dgdnu = Stellarator<Real,ORDER>::compute_pressure_jump(S_, pressure_, flux_tor_, flux_pol_, &g);
+      //Vector<ElemBasis> dgdnu = Stellarator<Real,ORDER>::compute_gradient(S_, pressure_, flux_tor_, flux_pol_, &g);
+      Vector<ElemBasis> dgdnu = Stellarator<Real,ORDER>::compute_pressure_jump(S_, pressure_, flux_tor_, flux_pol_, &g);
 
-      Vector<Real> grad_direction_, grad_op_;
+      Vector<Real> grad_direction_, grad_op_, grad_direction_orig_;
       { // Set grad_direction_ and filter
         Vector<ElemBasis> H = compute_H(S_.GetElemList(), normal);
         Vector<Real> H_fourier = cheb2grid(S_, H);
 
         grad_op_.ReInit(x.size());
         grad_direction_.ReInit(x.size());
+        grad_direction_orig_.ReInit(x.size());
         Vector<Real> dgdnu_fourier = cheb2grid(S_, dgdnu);
         for (Long i = 0; i < S_.Nsurf(); i++) { // Init grad_direction_, make it divergence-free, filter
           const Long Mt = S_.NTor(i);
           const Long Mp = S_.NPol(i);
           const Long offset = S_.ElemDsp(i);
-          const Long Nt = Mt * ORDER * fourier_upsample;
-          const Long Np = Mp * ORDER * fourier_upsample;
+          const Long Nt = fourier_dim0;
+          const Long Np = fourier_dim1;
 
-          const Vector<Real> dgdnu(Nt*Np, (Iterator<Real>)dgdnu_fourier.begin() + offset * (ORDER*ORDER*fourier_upsample*fourier_upsample), false);
-          const Vector<Real> H(Nt*Np, (Iterator<Real>)H_fourier.begin() + offset * (ORDER*ORDER*fourier_upsample*fourier_upsample), false);
-          const Vector<Real> X(COORD_DIM*Nt*Np, (Iterator<Real>)X_fourier.begin() + COORD_DIM*offset * (ORDER*ORDER*fourier_upsample*fourier_upsample), false);
-          Vector<Real> grad_direction(COORD_DIM*Nt*Np, (Iterator<Real>)grad_direction_.begin() + COORD_DIM*offset * (ORDER*ORDER*fourier_upsample*fourier_upsample), false);
-          Vector<Real> grad_op(COORD_DIM*Nt*Np, (Iterator<Real>)grad_op_.begin() + COORD_DIM*offset * (ORDER*ORDER*fourier_upsample*fourier_upsample), false);
+          const Vector<Real>    dgdnu(          Nt*Np, (Iterator<Real>)dgdnu_fourier.begin()   +           i * (fourier_dim0*fourier_dim1), false);
+          const Vector<Real>        H(          Nt*Np, (Iterator<Real>)H_fourier.begin()       +           i * (fourier_dim0*fourier_dim1), false);
+          const Vector<Real>        X(COORD_DIM*Nt*Np, (Iterator<Real>)X_fourier.begin()       + COORD_DIM*i * (fourier_dim0*fourier_dim1), false);
+          Vector<Real> grad_direction(COORD_DIM*Nt*Np, (Iterator<Real>)grad_direction_.begin() + COORD_DIM*i * (fourier_dim0*fourier_dim1), false);
+          Vector<Real>        grad_op(COORD_DIM*Nt*Np, (Iterator<Real>)grad_op_.begin()        + COORD_DIM*i * (fourier_dim0*fourier_dim1), false);
 
           Vector<Real> dX, d2X, normal, area_elem;
           biest::SurfaceOp<Real> Sop(comm, Nt, Np);
@@ -6085,31 +6467,38 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
 
           if (i < S_.Nsurf() - 1) { // Set grad_direction
             Vector<Real> F(Nt*Np), GradInvLapF;
-            for (Long i = 0; i < Nt*Np; i++) { // Set F <-- 2H * dgdnu
-              F[i] = 2*H[i] * dgdnu[i];
+            for (Long j = 0; j < Nt*Np; j++) { // Set F <-- 2H * dgdnu
+              F[j] = 2*H[j] * dgdnu[j];
             }
             Sop.GradInvSurfLap(GradInvLapF, dX, d2X, F, 1e-8, 100, 1.0);
 
-            for (Long i = 0; i < Nt*Np; i++) { // grad_direction <-- normal * dgdnu - GradInvLapF
-              grad_direction[0*Nt*Np+i] = normal[0*Nt*Np+i] * dgdnu[i] - GradInvLapF[0*Nt*Np+i]*0;
-              grad_direction[1*Nt*Np+i] = normal[1*Nt*Np+i] * dgdnu[i] - GradInvLapF[1*Nt*Np+i]*0;
-              grad_direction[2*Nt*Np+i] = normal[2*Nt*Np+i] * dgdnu[i] - GradInvLapF[2*Nt*Np+i]*0;
+            for (Long j = 0; j < Nt*Np; j++) { // grad_direction <-- normal * dgdnu - GradInvLapF
+              grad_direction[0*Nt*Np+j] = normal[0*Nt*Np+j] * dgdnu[j] - GradInvLapF[0*Nt*Np+j]*0;
+              grad_direction[1*Nt*Np+j] = normal[1*Nt*Np+j] * dgdnu[j] - GradInvLapF[1*Nt*Np+j]*0;
+              grad_direction[2*Nt*Np+j] = normal[2*Nt*Np+j] * dgdnu[j] - GradInvLapF[2*Nt*Np+j]*0;
             }
-            fourier_filter(grad_direction, Nt, Np, 0.1/fourier_upsample);
+            { ////////////////////
+              Vector<Real> grad_direction_orig(COORD_DIM*Nt*Np, (Iterator<Real>)grad_direction_orig_.begin() + COORD_DIM*i * (fourier_dim0*fourier_dim1), false);
+              grad_direction_orig = grad_direction;
+            }
+            fourier_filter(grad_direction, Nt, Np, 2);
 
-            for (Long k = 0; k < 20; k++) { // reparameterize
+            for (Long k = 0; k < 50; k++) { // reparameterize
+              Real max_resid = 0;
               Vector<Real> correction(COORD_DIM*Nt*Np);
               for (Long i = 0; i < Nt*Np; i++) { // Set correction
                 Real resid = dgdnu[i];
                 resid -= normal[0*Nt*Np+i] * grad_direction[0*Nt*Np+i];
                 resid -= normal[1*Nt*Np+i] * grad_direction[1*Nt*Np+i];
                 resid -= normal[2*Nt*Np+i] * grad_direction[2*Nt*Np+i];
+                max_resid = std::max(max_resid, fabs(resid));
 
                 correction[0*Nt*Np+i] = normal[0*Nt*Np+i] * resid;
                 correction[1*Nt*Np+i] = normal[1*Nt*Np+i] * resid;
                 correction[2*Nt*Np+i] = normal[2*Nt*Np+i] * resid;
               }
-              fourier_filter(correction, Nt, Np, 0.1/fourier_upsample);
+              std::cout<<max_resid<<' '; //////////////////////////////////
+              fourier_filter(correction, Nt, Np, 2);
 
               Real alpha = 0;
               { // Set alpha <-- (dgdnu - x.n, c.n) / (c.n, c.n)
@@ -6137,20 +6526,34 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
               }
               grad_direction += correction * alpha;
             }
-            //fourier_print_spectrum(dgdnu, Nt, Np);
-            //fourier_print_spectrum(grad_direction, Nt, Np);
+            //fourier_print_spectrum("dgdnu", dgdnu, Nt, Np);
+            //fourier_print_spectrum("grad_dir", grad_direction, Nt, Np);
           } else {
             grad_direction = 0;
           }
         }
       }
-      SCTL_ASSERT(grad_direction.size() == grad_direction_.Dim());
-      for (Long i = 0; i < grad_direction.size(); i++) { // Set grad_direction
-        grad_direction(i) = grad_direction_[i];
+
+      /////////////////////////////////////////////////////////////////////////
+      fourier_print_spectrum("X", X_fourier, S_);
+      fourier_print_spectrum("normal", cheb2grid(S_,normal), S_);
+
+      fourier_print_spectrum("dgdnu", cheb2grid(S_, dgdnu), S_);
+      fourier_print_spectrum("grad_dir", grad_direction_, S_);
+      fourier_print_spectrum("grad_dir_orig", grad_direction_orig_, S_);
+      /////////////////////////////////////////////////////////////////////////
+
+      { // Set grad_direction <-- grad_direction_
+        if (grad_direction.size() != grad_direction_.Dim()) {
+          grad_direction = Eigen::VectorXd(grad_direction_.Dim());
+        }
+        for (Long i = 0; i < grad_direction.size(); i++) {
+          grad_direction(i) = grad_direction_[i];
+        }
       }
       if (grad_op_ptr) { // Set grad_op_ptr
         grad_op_ptr[0] = Eigen::VectorXd(grad_op_.Dim());
-        for (Long i = 0; i < grad_op_ptr->size(); i++) { // Set grad_direction
+        for (Long i = 0; i < grad_op_ptr->size(); i++) {
           grad_op_ptr[0](i) = grad_op_[i];
         }
       }
@@ -6195,14 +6598,14 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
         }
         Vector<Real> X_fourier = cheb2grid(mhd_equilib.S_, X);
         x.resize(X_fourier.Dim());
-        // X_fourier.Read(("X"+std::to_string(0)).c_str()); // Read from file
+        X_fourier.Read(("X"+std::to_string(mhd_equilib.iter)+"_").c_str()); // Read from file
         for (Long i = 0; i < X_fourier.Dim(); i++) {
           x(i) = X_fourier[i];
         }
       }
 
       Real fx;
-      if (1) {
+      if (0) {
         LBFGSpp::LBFGSParam<Real> param;
         param.max_iterations = 200;
         param.epsilon = 1e-16;
@@ -6210,7 +6613,8 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
         LBFGSpp::LBFGSSolver<Real> solver(param);
         Integer niter = solver.minimize(mhd_equilib, x, fx);
       } else {
-        Integer niter = GradientDescent2_(mhd_equilib, x, fx, 200, 1e-12);
+        //Integer niter = GradientDescentNew(mhd_equilib, x, fx, 200, 1e-12);
+        Integer niter = GradientDescent(mhd_equilib, x, fx, 200, 1e-12);
       }
       { // Set x
         // TODO
@@ -6227,10 +6631,10 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       { // Init S, flux_tor, flux_pol, pressure
         Vector<Long> NtNp;
         for (Long i = 0; i < Nsurf; i++) {
-          //NtNp.PushBack(50);
-          //NtNp.PushBack(8);
-          NtNp.PushBack(30);
-          NtNp.PushBack(4);
+          NtNp.PushBack(70);
+          NtNp.PushBack(14);
+          //NtNp.PushBack(30);
+          //NtNp.PushBack(4);
         }
         S = Stellarator<Real,ORDER>(NtNp);
         flux_tor = 1;
@@ -6249,12 +6653,68 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       ComputeEquilibrium(mhd_equilib);
     }
 
+    static void test_() {
+      Comm comm = Comm::World();
+      Profile::Enable(true);
+
+      Long Nsurf = 2;
+      Stellarator<Real,ORDER> S;
+      Vector<Real> flux_tor(Nsurf), flux_pol(Nsurf), pressure(Nsurf);
+      { // Init S, flux_tor, flux_pol, pressure
+        Vector<Long> NtNp;
+        for (Long i = 0; i < Nsurf; i++) {
+          NtNp.PushBack(50);
+          NtNp.PushBack(10);
+          //NtNp.PushBack(30);
+          //NtNp.PushBack(4);
+        }
+        S = Stellarator<Real,ORDER>(NtNp);
+        flux_tor = 1;
+        flux_pol = 1;
+        pressure = 0;
+
+        //flux_tor[0] = 1; //0.791881512;
+        //flux_tor[1] = 1;
+        //flux_pol[0] = 0;
+        //flux_pol[1] = 0;
+        //pressure[0] = 0;
+        //pressure[1] = 0;
+      }
+      MHDEquilib mhd_equilib(S, pressure, flux_tor, flux_pol);
+
+      const Long Nelem = mhd_equilib.S_.NElem();
+      const Long Nnodes = ElemBasis::Size();
+      Eigen::VectorXd x, grad_direction;
+      { // Read x
+        Vector<ElemBasis> X(Nelem * COORD_DIM);
+        for (Long i = 0; i < Nelem; i++) {
+          X[i*COORD_DIM+0] = S.Elem(i,0);
+          X[i*COORD_DIM+1] = S.Elem(i,1);
+          X[i*COORD_DIM+2] = S.Elem(i,2);
+        }
+        Vector<Real> X_fourier = cheb2grid(S, X);
+        //X_fourier.Read("X_tmp");
+        x.resize(X_fourier.Dim());
+        for (Long i = 0; i < X_fourier.Dim(); i++) {
+          x(i) = X_fourier[i];
+        }
+      }
+      Real g = mhd_equilib(x, grad_direction);
+      { // Write grad_direction
+        Vector<Real> dXdt(grad_direction.size());
+        for (Long i = 0; i < dXdt.Dim(); i++) {
+          dXdt[i] = grad_direction(i);
+        }
+        dXdt.Write("dXdt_tmp");
+      }
+    }
+
   private:
     Stellarator<Real,ORDER> S_;
     Vector<Real> pressure_;
     Vector<Real> flux_tor_;
     Vector<Real> flux_pol_;
-    Long iter = 0;
+    Long iter;
 };