.

include/sctl/boundary_quadrature.hpp

@@ -2841,11 +2841,11 @@ template <class Real, Integer ORDER=10> class Stellarator {
       SCTL_ASSERT(Nsurf*2 == NtNp_.Dim());
 
       Long Nelem = 0;
-      elem_dsp.ReInit(Nsurf);
-      if (elem_dsp.Dim()) elem_dsp[0] = 0;
+      elem_dsp.ReInit(Nsurf+1);
+      elem_dsp[0] = 0;
       for (Long i = 0; i < Nsurf; i++) {
         Nelem += NtNp_[i*2+0]*NtNp_[i*2+1];
-        if (i+1 < Nsurf) elem_dsp[i+1] = Nelem;
+        elem_dsp[i+1] = Nelem;
       }
       elements.ReInit(Nelem);
       for (Long i = 0; i < Nsurf; i++) {
@@ -2854,7 +2854,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
     }
 
     Long ElemIdx(Long s, Long t, Long p) {
-      SCTL_ASSERT(0 <= s && s < elem_dsp.Dim());
+      SCTL_ASSERT(0 <= s && s < Nsurf());
       SCTL_ASSERT(0 <= t && t < NtNp_[s*2+0]);
       SCTL_ASSERT(0 <= p && p < NtNp_[s*2+1]);
       return elem_dsp[s] + t*NtNp_[s*2+1] + p;
@@ -2870,7 +2870,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
     }
 
     Long Nsurf() const {
-      return elem_dsp.Dim();
+      return elem_dsp.Dim()-1;
     }
     Long ElemDsp(Long s) const {
       return elem_dsp[s];
@@ -2935,12 +2935,12 @@ 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);
-        { // Write VTU
+        if (0) { // Write VTU
           VTUData vtu;
           vtu.AddElems(S.GetElemList(), S.sigma, ORDER);
           vtu.WriteVTK("sigma"+std::to_string(i), comm);
         }
-        { // Write VTU
+        if (0) { // Write VTU
           VTUData vtu;
           vtu.AddElems(S.GetElemList(), S.B, ORDER);
           vtu.WriteVTK("B"+std::to_string(i), comm);
@@ -3567,7 +3567,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           const Long elem_dsp = (i==0 ? 0 : S_.ElemDsp(i-1));
           const Long Nnodes = ElemBasis::Size();
           auto dgdnu_ = compute_gradient(Svec[i]);
-          { // Write VTU
+          if (0) { // Write VTU
             VTUData vtu;
             vtu.AddElems(Svec[i].GetElemList(), dgdnu_, ORDER);
             vtu.WriteVTK("dgdnu-"+std::to_string(i), comm);
@@ -3731,10 +3731,91 @@ template <class Real, Integer ORDER=10> class Stellarator {
         //pressure[1] = 0;
       }
 
-      //{ // find equilibrium flux surfaces
+      { // find equilibrium flux surfaces
+        const Long Nnodes = ElemBasis::Size();
+        const Long Nelem = S.NElem();
+        const Long Nsurf = S.Nsurf();
 
-      //  return;
-      //}
+        Long iter = 0;
+        Real dt = 0.1;
+        while (1) { // time-step
+          Vector<ElemBasis> normal, area_elem;
+          compute_norm_area_elem(S, normal, area_elem);
+          Vector<ElemBasis> dgdnu = compute_gradient(S, pressure, flux_tor, flux_pol)*(-1);
+          //Vector<ElemBasis> dgdnu = compute_pressure_jump(S, pressure, flux_tor, flux_pol)*(-1);
+
+          Vector<ElemBasis> dXdt(dgdnu.Dim()*COORD_DIM);
+          { // Set dXdt
+            dXdt = 0;
+            for (Long i = 0; i < S.ElemDsp(Nsurf-1); i++) {
+              for (Long j = 0; j < Nnodes; j++) {
+                dXdt[i*COORD_DIM+0][j] = normal[i*COORD_DIM+0][j] * dgdnu[i][j];
+                dXdt[i*COORD_DIM+1][j] = normal[i*COORD_DIM+1][j] * dgdnu[i][j];
+                dXdt[i*COORD_DIM+2][j] = normal[i*COORD_DIM+2][j] * dgdnu[i][j];
+              }
+            }
+          }
+
+          { // Update S, dt
+            Stellarator<Real,ORDER> S0 = S, S1 = S, S2 = S;
+            for (Long i = 0; i < S.NElem(); i++) {
+              for (Long j = 0; j < Nnodes; j++) {
+                S0.Elem(i, 0)[j] += 0.0 * dt * dXdt[i*COORD_DIM+0][j];
+                S0.Elem(i, 1)[j] += 0.0 * dt * dXdt[i*COORD_DIM+1][j];
+                S0.Elem(i, 2)[j] += 0.0 * dt * dXdt[i*COORD_DIM+2][j];
+
+                S1.Elem(i, 0)[j] += 0.5 * dt * dXdt[i*COORD_DIM+0][j];
+                S1.Elem(i, 1)[j] += 0.5 * dt * dXdt[i*COORD_DIM+1][j];
+                S1.Elem(i, 2)[j] += 0.5 * dt * dXdt[i*COORD_DIM+2][j];
+
+                S2.Elem(i, 0)[j] += 1.0 * dt * dXdt[i*COORD_DIM+0][j];
+                S2.Elem(i, 1)[j] += 1.0 * dt * dXdt[i*COORD_DIM+1][j];
+                S2.Elem(i, 2)[j] += 1.0 * dt * dXdt[i*COORD_DIM+2][j];
+              }
+            }
+            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);
+            { // Calculate optimal step size dt
+              Real a = 2*g0 - 4*g1 + 2*g2;
+              Real b =-3*g0 + 4*g1 -   g2;
+              Real c = g0;
+              Real s = -b/(2*a);
+              dt *= s;
+
+              Real g_ = a*s*s + b*s + c;
+              std::cout<<"g = "<<g_<<' ';
+              std::cout<<g0<<' ';
+              std::cout<<g1<<' ';
+              std::cout<<g2<<' ';
+              std::cout<<dt<<'\n';
+            }
+
+            for (Long i = 0; i < S.NElem(); i++) {
+              for (Long j = 0; j < Nnodes; j++) {
+                S.Elem(i, 0)[j] += dt * dXdt[i*COORD_DIM+0][j];
+                S.Elem(i, 1)[j] += dt * dXdt[i*COORD_DIM+1][j];
+                S.Elem(i, 2)[j] += dt * dXdt[i*COORD_DIM+2][j];
+              }
+            }
+          }
+          iter++;
+
+          Vector<ElemBasis> pressure_jump = compute_pressure_jump(S, pressure, flux_tor, flux_pol);
+          { // Write VTU
+            VTUData vtu;
+            vtu.AddElems(S.GetElemList(), dgdnu, ORDER);
+            vtu.WriteVTK("dgdnu"+std::to_string(iter), comm);
+          }
+          { // Write VTU
+            VTUData vtu;
+            vtu.AddElems(S.GetElemList(), pressure_jump, ORDER);
+            vtu.WriteVTK("pressure_jump"+std::to_string(iter), comm);
+          }
+        }
+
+        return;
+      }
 
       { // Verify using finite difference approximation
         Vector<ElemBasis> dgdnu = compute_gradient(S, pressure, flux_tor, flux_pol);