.

include/sctl/boundary_quadrature.hpp

@@ -2050,7 +2050,7 @@ template <class Real> class Quadrature {
 
 template <class Real, Integer ORDER=10> class Stellarator {
   private:
-    static constexpr Integer order_singular = 15;
+    static constexpr Integer order_singular = 25;
     static constexpr Integer order_direct = 35;
     static constexpr Integer COORD_DIM = 3;
     static constexpr Integer ELEM_DIM = COORD_DIM-1;
@@ -3686,7 +3686,131 @@ template <class Real, Integer ORDER=10> class Stellarator {
       }
 
       { // find equilibrium flux surfaces
-        auto filter = [](const Stellarator<Real,ORDER>& S, Vector<ElemBasis>& f) {
+        {
+        //auto filter = [](const Stellarator<Real,ORDER>& S, Vector<ElemBasis>& f) {
+        //  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);
+
+        //    Vector<Real> Xf(dof*Nt*Np); Xf = 0;
+        //    const Long Nnodes = ElemBasis::Size();
+        //    const Matrix<Real>& Mnodes = Basis<Real,1,ORDER>::Nodes();
+        //    for (Long t = 0; t < Nt; t++) {
+        //      for (Long p = 0; p < Np; p++) {
+        //        Real theta = t / (Real)Nt;
+        //        Real phi   = p / (Real)Np;
+        //        Long i = (Long)(theta * Mt);
+        //        Long j = (Long)(phi   * Mp);
+        //        Real x = theta * Mt - i;
+        //        Real y = phi   * Mp - j;
+        //        Long elem_idx = i * Mp + j;
+
+        //        Vector<Real> Interp0(ORDER);
+        //        Vector<Real> Interp1(ORDER);
+        //        { // Set Interp0, Interp1
+        //          auto node = [&Mnodes] (Long i) {
+        //            return Mnodes[0][i];
+        //          };
+        //          for (Long i = 0; i < ORDER; i++) {
+        //            Real wt_x = 1, wt_y = 1;
+        //            for (Long j = 0; j < ORDER; j++) {
+        //              if (j != i) {
+        //                wt_x *= (x - node(j)) / (node(i) - node(j));
+        //                wt_y *= (y - node(j)) / (node(i) - node(j));
+        //              }
+        //              Interp0[i] = wt_x;
+        //              Interp1[i] = wt_y;
+        //            }
+        //          }
+        //        }
+
+        //        for (Long ii = 0; ii < ORDER; ii++) {
+        //          for (Long jj = 0; jj < ORDER; jj++) {
+        //            Long node_idx = jj * ORDER + ii;
+        //            for (Long k = 0; k < dof; k++) {
+        //              Xf[(k*Nt+t)*Np+p] += X[elem_idx*dof+k][node_idx] * Interp0[ii] * Interp1[jj];
+        //            }
+        //          }
+        //        }
+        //      }
+        //    }
+        //    return Xf;
+        //  };
+        //  auto grid2cheb = [] (const Vector<Real>& Xf, Long Nt, Long Np, Long Mt, Long Mp) {
+        //    Long dof = Xf.Dim() / (Nt*Np);
+        //    SCTL_ASSERT(Xf.Dim() == dof*Nt*Np);
+        //    Vector<ElemBasis> X(Mt*Mp*dof);
+        //    constexpr Integer INTERP_ORDER = 12;
+
+        //    for (Long tt = 0; tt < Mt; tt++) {
+        //      for (Long pp = 0; pp < Mp; pp++) {
+        //        for (Long t = 0; t < ORDER; t++) {
+        //          for (Long p = 0; p < ORDER; p++) {
+        //            Matrix<Real> Mnodes = Basis<Real,1,ORDER>::Nodes();
+        //            Real theta = (tt + Mnodes[0][t]) / Mt;
+        //            Real phi   = (pp + Mnodes[0][p]) / Mp;
+        //            Long i = (Long)(theta * Nt);
+        //            Long j = (Long)(phi   * Np);
+        //            Real x = theta * Nt - i;
+        //            Real y = phi   * Np - j;
+
+        //            Vector<Real> Interp0(INTERP_ORDER);
+        //            Vector<Real> Interp1(INTERP_ORDER);
+        //            { // Set Interp0, Interp1
+        //              auto node = [] (Long i) {
+        //                return (Real)i - (INTERP_ORDER-1)/2;
+        //              };
+        //              for (Long i = 0; i < INTERP_ORDER; i++) {
+        //                Real wt_x = 1, wt_y = 1;
+        //                for (Long j = 0; j < INTERP_ORDER; j++) {
+        //                  if (j != i) {
+        //                    wt_x *= (x - node(j)) / (node(i) - node(j));
+        //                    wt_y *= (y - node(j)) / (node(i) - node(j));
+        //                  }
+        //                  Interp0[i] = wt_x;
+        //                  Interp1[i] = wt_y;
+        //                }
+        //              }
+        //            }
+
+        //            for (Long k = 0; k < dof; k++) {
+        //              Real X0 = 0;
+        //              for (Long ii = 0; ii < INTERP_ORDER; ii++) {
+        //                for (Long jj = 0; jj < INTERP_ORDER; jj++) {
+        //                  Long idx_i = (i + ii-(INTERP_ORDER-1)/2 + Nt) % Nt;
+        //                  Long idx_j = (j + jj-(INTERP_ORDER-1)/2 + Np) % Np;
+        //                  X0 += Interp0[ii] * Interp1[jj] * Xf[(k*Nt+idx_i)*Np+idx_j];
+        //                }
+        //              }
+        //              Long elem_idx = tt * Mp + pp;
+        //              Long node_idx = p * ORDER + t;
+        //              X[elem_idx*dof+k][node_idx] = X0;
+        //            }
+        //          }
+        //        }
+        //      }
+        //    }
+        //    return X;
+        //  };
+
+        //  Long dof = f.Dim() / S.NElem();
+        //  SCTL_ASSERT(f.Dim() == S.NElem() * dof);
+        //  for (Long i = 0; i < S.Nsurf(); i++) {
+        //    const Long Mt = S.NTor(i);
+        //    const Long Mp = S.NPol(i);
+        //    const Long Nelem = Mt * Mp;
+        //    const Long offset = S.ElemDsp(i);
+        //    const Long Nt = Mt * ORDER / 5;
+        //    const Long Np = Mp * ORDER / 5;
+
+        //    Vector<ElemBasis> f_(Nelem*dof, f.begin() + offset*dof, false);
+        //    Vector<Real> f_fourier = cheb2grid(f_, Mt, Mp, Nt, Np);
+        //    f_ = grid2cheb(f_fourier, Nt, Np, Mt, Mp);
+        //  }
+        //};
+        }
+
+        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);
@@ -3791,6 +3915,40 @@ template <class Real, Integer ORDER=10> class Stellarator {
             }
             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);
@@ -3799,11 +3957,12 @@ template <class Real, Integer ORDER=10> class Stellarator {
             const Long Mp = S.NPol(i);
             const Long Nelem = Mt * Mp;
             const Long offset = S.ElemDsp(i);
-            const Long Nt = Mt * ORDER / 5;
-            const Long Np = Mp * ORDER / 5;
+            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);
           }
         };
@@ -3827,7 +3986,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
                 dXdt[i*COORD_DIM+2][j] = normal[i*COORD_DIM+2][j] * dgdnu[i][j];
               }
             }
-            filter(S, dXdt);
+            filter(S, comm, dXdt, 0.1);
           }
           { // Update dt
             const Long Nelem = S.NElem();
@@ -3886,11 +4045,11 @@ template <class Real, Integer ORDER=10> class Stellarator {
             const Long Nelem = S.NElem();
             Vector<ElemBasis> X(Nelem*COORD_DIM);
             for (Long i = 0; i < S.NElem(); i++) {
-              X[i*COORD_DIM+0] = S.Elem(i, 0) + dXdt[i*COORD_DIM+0] * dt * 0.5;
-              X[i*COORD_DIM+1] = S.Elem(i, 1) + dXdt[i*COORD_DIM+1] * dt * 0.5;
-              X[i*COORD_DIM+2] = S.Elem(i, 2) + dXdt[i*COORD_DIM+2] * dt * 0.5;
+              X[i*COORD_DIM+0] = S.Elem(i, 0) + dXdt[i*COORD_DIM+0] * dt;
+              X[i*COORD_DIM+1] = S.Elem(i, 1) + dXdt[i*COORD_DIM+1] * dt;
+              X[i*COORD_DIM+2] = S.Elem(i, 2) + dXdt[i*COORD_DIM+2] * dt;
             }
-            filter(S, X);
+            filter(S, comm, X, 0.3);
             for (Long i = 0; i < S.NElem(); i++) {
               S.Elem(i, 0) = X[i*COORD_DIM+0];
               S.Elem(i, 1) = X[i*COORD_DIM+1];
@@ -4915,8 +5074,8 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
           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);};
+          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_);
@@ -4939,7 +5098,7 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
 
         Long dof = f.Dim() / S.NElem();
         SCTL_ASSERT(f.Dim() == S.NElem() * dof);
-        for (Long i = 0; i < S.Nsurf(); i++) {
+        for (Long i = 0; i < S.Nsurf()-1; i++) {
           const Long Mt = S.NTor(i);
           const Long Mp = S.NPol(i);
           const Long Nelem = Mt * Mp;
@@ -4962,6 +5121,23 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
           S_.Elem(i,2)[j] = x[(i*Nnodes+j)*COORD_DIM+2];
         }
       }
+      Stellarator<Real,ORDER> SS; //////////////////////////
+      { // Update S <-- filter(S)
+        const Long Nelem = S_.NElem();
+        Vector<ElemBasis> X(Nelem*COORD_DIM);
+        for (Long i = 0; i < S_.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);
+        }
+        SS = S_;
+        filter(S_, comm, X, 0.1);
+        for (Long i = 0; i < S_.NElem(); i++) {
+          S_.Elem(i, 0) = X[i*COORD_DIM+0];
+          S_.Elem(i, 1) = X[i*COORD_DIM+1];
+          S_.Elem(i, 2) = X[i*COORD_DIM+2];
+        }
+      }
       Vector<ElemBasis> dgdnu = Stellarator<Real,ORDER>::compute_gradient(S_, pressure_, flux_tor_, flux_pol_, &g);
       Vector<ElemBasis> dXdt(Nelem*COORD_DIM);
       { // Set dXdt
@@ -4976,7 +5152,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_, comm, dXdt, 0.3333);
+        //filter(S_, comm, dXdt, 0.1);
       }
       for (Long i = 0; i < Nelem; i++) { // Set grad
         for (Long j = 0; j < Nnodes; j++) {
@@ -4996,6 +5172,11 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
         vtu.AddElems(S_.GetElemList(), dXdt, ORDER);
         vtu.WriteVTK("dXdt"+std::to_string(iter), comm);
       }
+      if (1) { // Write VTU
+        VTUData vtu;
+        vtu.AddElems(SS.GetElemList(), dgdnu, ORDER);
+        vtu.WriteVTK("S"+std::to_string(iter), comm);
+      }
       std::cout<<"iter = "<<iter<<"    g = "<<g<<'\n';
 
       iter++;
@@ -5008,7 +5189,7 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       const Long N = Nelem * COORD_DIM * Nnodes;
 
       LBFGSpp::LBFGSParam<Real> param;
-      param.epsilon = 1e-6;
+      param.epsilon = 1e-8;
       param.max_iterations = 100;
 
       // Create solver and function object