.

include/sctl/boundary_quadrature.hpp

@@ -2555,7 +2555,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
       if (S.Nsurf() >= 2) Ax[Nelem*Nnodes + 1] = flux_pol;
       return Ax;
     }
-    static void compute_invA(Vector<ElemBasis>& sigma, Real& alpha, Real& beta, const Stellarator<Real,ORDER>& S, Real flux_tor, Real flux_pol, const Comm& comm) {
+    static void compute_invA(Vector<ElemBasis>& sigma, Real& alpha, Real& beta, const Stellarator<Real,ORDER>& S, Vector<ElemBasis>& Bdotn, Real flux_tor, Real flux_pol, const Comm& comm) {
       typename sctl::ParallelSolver<Real>::ParallelOp BIOp = [&S](sctl::Vector<Real>* Ax, const sctl::Vector<Real>& x) {
         (*Ax) = compute_A(S, x);
       };
@@ -2563,7 +2563,11 @@ template <class Real, Integer ORDER=10> class Stellarator {
       const Long Nnodes = ElemBasis::Size();
 
       Vector<Real> rhs_(Nelem * Nnodes + S.Nsurf());
-      rhs_ = 0;
+      for (Long i = 0; i < Nelem; i++) {
+        for (Long j = 0; j < Nnodes; j++) {
+          rhs_[i*Nnodes+j] = Bdotn[i][j];
+        }
+      }
       if (S.Nsurf() >= 1) rhs_[Nelem * Nnodes + 0] = flux_tor;
       if (S.Nsurf() >= 2) rhs_[Nelem * Nnodes + 1] = flux_pol;
 
@@ -2580,6 +2584,11 @@ template <class Real, Integer ORDER=10> class Stellarator {
       }
       alpha = (S.Nsurf() >= 1 ? x_[Nelem * Nnodes + 0] : 0);
       beta  = (S.Nsurf() >= 2 ? x_[Nelem * Nnodes + 1] : 0);
+    };
+    static void compute_invA(Vector<ElemBasis>& sigma, Real& alpha, Real& beta, const Stellarator<Real,ORDER>& S, Real flux_tor, Real flux_pol, const Comm& comm) {
+      Vector<ElemBasis> Bdotn(S.NElem());
+      Bdotn = 0;
+      compute_invA(sigma, alpha, beta, S, Bdotn, flux_tor, flux_pol, comm);
     }
     static Vector<Real> compute_Aadj(const Stellarator<Real,ORDER>& S, const Vector<Real>& x) {
       const Long Nelem = S.NElem();
@@ -4388,6 +4397,240 @@ template <class Real, Integer ORDER=10> class Stellarator {
           EvalQuadrature(S.dBp0, S.quadrature_dBS, S, Jt, S.BiotSavartGrad);
         }
       };
+      auto compute_grad = [] (const Stellarator<Real,ORDER>& S, const Vector<ElemBasis>& V) {
+        const Long Nelem = S.GetElemList().NElem();
+        const Long Nnodes = ElemBasis::Size();
+        const Long dof = V.Dim() / Nelem;
+        SCTL_ASSERT(Nelem * dof == V.Dim());
+
+        Vector<ElemBasis> du_dX(Nelem*COORD_DIM*2);
+        { // Set du_dX
+          Vector<ElemBasis> dX;
+          ElemBasis::Grad(dX, S.GetElemList().ElemVector());
+
+          auto inv2x2 = [](Tensor<Real, true, 2, 2> M) {
+            Tensor<Real, true, 2, 2> Mout;
+            Real oodet = 1 / (M(0,0) * M(1,1) - M(0,1) * M(1,0));
+            Mout(0,0) =  M(1,1) * oodet;
+            Mout(0,1) = -M(0,1) * oodet;
+            Mout(1,0) = -M(1,0) * oodet;
+            Mout(1,1) =  M(0,0) * oodet;
+            return Mout;
+          };
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              Tensor<Real, true, 3, 2> dX_du;
+              dX_du(0,0) = dX[(i*COORD_DIM+0)*2+0][j];
+              dX_du(1,0) = dX[(i*COORD_DIM+1)*2+0][j];
+              dX_du(2,0) = dX[(i*COORD_DIM+2)*2+0][j];
+              dX_du(0,1) = dX[(i*COORD_DIM+0)*2+1][j];
+              dX_du(1,1) = dX[(i*COORD_DIM+1)*2+1][j];
+              dX_du(2,1) = dX[(i*COORD_DIM+2)*2+1][j];
+
+              Tensor<Real, true, 2, 2> G; // = dX_du.Transpose() * dX_du;
+              G(0,0) = dX_du(0,0) * dX_du(0,0) + dX_du(1,0) * dX_du(1,0) + dX_du(2,0) * dX_du(2,0);
+              G(0,1) = dX_du(0,0) * dX_du(0,1) + dX_du(1,0) * dX_du(1,1) + dX_du(2,0) * dX_du(2,1);
+              G(1,0) = dX_du(0,1) * dX_du(0,0) + dX_du(1,1) * dX_du(1,0) + dX_du(2,1) * dX_du(2,0);
+              G(1,1) = dX_du(0,1) * dX_du(0,1) + dX_du(1,1) * dX_du(1,1) + dX_du(2,1) * dX_du(2,1);
+
+              Tensor<Real, true, 2, 2> Ginv = inv2x2(G);
+              du_dX[(i*COORD_DIM+0)*2+0][j] = Ginv(0,0) * dX_du(0,0) + Ginv(0,1) * dX_du(0,1);
+              du_dX[(i*COORD_DIM+1)*2+0][j] = Ginv(0,0) * dX_du(1,0) + Ginv(0,1) * dX_du(1,1);
+              du_dX[(i*COORD_DIM+2)*2+0][j] = Ginv(0,0) * dX_du(2,0) + Ginv(0,1) * dX_du(2,1);
+              du_dX[(i*COORD_DIM+0)*2+1][j] = Ginv(1,0) * dX_du(0,0) + Ginv(1,1) * dX_du(0,1);
+              du_dX[(i*COORD_DIM+1)*2+1][j] = Ginv(1,0) * dX_du(1,0) + Ginv(1,1) * dX_du(1,1);
+              du_dX[(i*COORD_DIM+2)*2+1][j] = Ginv(1,0) * dX_du(2,0) + Ginv(1,1) * dX_du(2,1);
+            }
+          }
+        }
+
+        Vector<ElemBasis> dV;
+        ElemBasis::Grad(dV, V);
+
+        Vector<ElemBasis> gradV(Nelem*dof*COORD_DIM);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            for (Long k = 0; k < dof; k++) {
+              gradV[(i*dof+k)*COORD_DIM+0][j] = dV[(i*dof+k)*2+0][j] * du_dX[(i*COORD_DIM+0)*2+0][j] + dV[(i*dof+k)*2+1][j] * du_dX[(i*COORD_DIM+0)*2+1][j];
+              gradV[(i*dof+k)*COORD_DIM+1][j] = dV[(i*dof+k)*2+0][j] * du_dX[(i*COORD_DIM+1)*2+0][j] + dV[(i*dof+k)*2+1][j] * du_dX[(i*COORD_DIM+1)*2+1][j];
+              gradV[(i*dof+k)*COORD_DIM+2][j] = dV[(i*dof+k)*2+0][j] * du_dX[(i*COORD_DIM+2)*2+0][j] + dV[(i*dof+k)*2+1][j] * du_dX[(i*COORD_DIM+2)*2+1][j];
+            }
+          }
+        }
+        return gradV;
+      };
+      auto compute_surfdiv = [&compute_grad] (const Stellarator<Real,ORDER>& S, const Vector<ElemBasis>& V) {
+        const Long Nelem = S.GetElemList().NElem();
+        const Long Nnodes = ElemBasis::Size();
+        SCTL_ASSERT(V.Dim() == Nelem* COORD_DIM);
+
+        Vector<ElemBasis> gradV = compute_grad(S, V);
+        Vector<ElemBasis> divV(Nelem);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            divV[i][j] = gradV[(i*COORD_DIM+0)*COORD_DIM+0][j] + gradV[(i*COORD_DIM+1)*COORD_DIM+1][j] + gradV[(i*COORD_DIM+2)*COORD_DIM+2][j];
+          }
+        }
+        return divV;
+      };
+      auto compute_g = [](const Stellarator<Real,ORDER>& S, const Vector<ElemBasis>& B) {
+        const Long Nelem = S.NElem();
+        const Long Nnodes = ElemBasis::Size();
+        Vector<ElemBasis> normal, area_elem;
+        compute_norm_area_elem(S, normal, area_elem);
+        Vector<ElemBasis> B2(Nelem);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            B2[i][j] = 0;
+            B2[i][j] += B[i*COORD_DIM+0][j] * B[i*COORD_DIM+0][j];
+            B2[i][j] += B[i*COORD_DIM+1][j] * B[i*COORD_DIM+1][j];
+            B2[i][j] += B[i*COORD_DIM+2][j] * B[i*COORD_DIM+2][j];
+          }
+        }
+        return compute_inner_prod(area_elem,B2, B2) * 0.25;
+      };
+      auto compute_H = [] (const ElemList<COORD_DIM,ElemBasis>& elem_lst, const Vector<ElemBasis>& normal) {
+        const Long Nnodes = ElemBasis::Size();
+        const Long Nelem = elem_lst.NElem();
+
+        const Vector<ElemBasis> X = elem_lst.ElemVector();
+        Vector<ElemBasis> dX, d2X, H(Nelem);
+        ElemBasis::Grad(dX, X);
+        ElemBasis::Grad(d2X, dX);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            Tensor<Real,true,2,2> I, invI, II;
+            for (Long k0 = 0; k0 < 2; k0++) {
+              for (Long k1 = 0; k1 < 2; k1++) {
+                I(k0,k1) = 0;
+                I(k0,k1) += dX[(i*COORD_DIM+0)*2+k0][j] * dX[(i*COORD_DIM+0)*2+k1][j];
+                I(k0,k1) += dX[(i*COORD_DIM+1)*2+k0][j] * dX[(i*COORD_DIM+1)*2+k1][j];
+                I(k0,k1) += dX[(i*COORD_DIM+2)*2+k0][j] * dX[(i*COORD_DIM+2)*2+k1][j];
+
+                II(k0,k1) = 0;
+                II(k0,k1) += d2X[(i*COORD_DIM+0)*4+k0*2+k1][j] * normal[i*COORD_DIM+0][j];
+                II(k0,k1) += d2X[(i*COORD_DIM+1)*4+k0*2+k1][j] * normal[i*COORD_DIM+1][j];
+                II(k0,k1) += d2X[(i*COORD_DIM+2)*4+k0*2+k1][j] * normal[i*COORD_DIM+2][j];
+              }
+            }
+            { // Set invI
+              Real detI = I(0,0)*I(1,1)-I(0,1)*I(1,0);
+              invI(0,0) = I(1,1) / detI;
+              invI(0,1) = -I(0,1) / detI;
+              invI(1,0) = -I(1,0) / detI;
+              invI(1,1) = I(0,0) / detI;
+            }
+            { // Set H
+              H[i][j] = 0;
+              H[i][j] += -0.5 * II(0,0)*invI(0,0);
+              H[i][j] += -0.5 * II(0,1)*invI(0,1);
+              H[i][j] += -0.5 * II(1,0)*invI(1,0);
+              H[i][j] += -0.5 * II(1,1)*invI(1,1);
+            }
+          }
+        }
+        return H;
+      };
+      auto compute_A_ = [](const Stellarator<Real,ORDER>& S, const Vector<Real>& x) {
+        const Long Nelem = S.NElem();
+        const Long Nnodes = ElemBasis::Size();
+        SCTL_ASSERT(x.Dim() == Nelem*Nnodes+S.Nsurf());
+
+        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;
+            }
+          }
+        }
+
+        Vector<ElemBasis> sigma(Nelem);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            sigma[i][j] = x[i*Nnodes+j];
+          }
+        }
+        Real alpha = (S.Nsurf() >= 1 ? x[Nelem*Nnodes + 0] : 0);
+        Real beta  = (S.Nsurf() >= 2 ? x[Nelem*Nnodes + 1] : 0);
+
+        Vector<ElemBasis> B = compute_B(S, sigma, alpha, beta);
+        Vector<ElemBasis> BdotN = compute_dot_prod(B, normal);
+
+        Real flux_tor, flux_pol;
+        compute_flux(flux_tor, flux_pol, S, B, normal);
+        { // update flux_tor
+          Vector<ElemBasis> G_BdotN(Nelem), phi_dot_N_over_R(Nelem);
+          EvalQuadrature(G_BdotN, S.quadrature_FxU, S, BdotN, S.Laplace_FxU);
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              Tensor<Real,true,COORD_DIM> x, b0, axis;
+              x(0) = S.Elem(i,0)[j];
+              x(1) = S.Elem(i,1)[j];
+              x(2) = S.Elem(i,2)[j];
+
+              axis(0) = 0;
+              axis(1) = 0;
+              axis(2) = 1;
+              b0(0) = axis(1) * x(2) - axis(2) * x(1);
+              b0(1) = axis(2) * x(0) - axis(0) * x(2);
+              b0(2) = axis(0) * x(1) - axis(1) * x(0);
+              Real scale = 1 / (b0(0)*b0(0) + b0(1)*b0(1) + b0(2)*b0(2));
+              b0(0) *= scale;
+              b0(1) *= scale;
+              b0(2) *= scale;
+
+              phi_dot_N_over_R[i][j] = 0;
+              phi_dot_N_over_R[i][j] += normal[i*COORD_DIM+0][j] * b0(0);
+              phi_dot_N_over_R[i][j] += normal[i*COORD_DIM+1][j] * b0(1);
+              phi_dot_N_over_R[i][j] += normal[i*COORD_DIM+2][j] * b0(2);
+            }
+          }
+          flux_tor += compute_inner_prod(area_elem, phi_dot_N_over_R, G_BdotN)/(2*const_pi<Real>());
+        }
+
+        Vector<Real> Ax(Nelem*Nnodes+S.Nsurf());
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            Ax[i*Nnodes+j] = BdotN[i][j];
+          }
+        }
+        if (S.Nsurf() >= 1) Ax[Nelem*Nnodes + 0] = flux_tor;
+        if (S.Nsurf() >= 2) Ax[Nelem*Nnodes + 1] = flux_pol;
+        return Ax;
+      };
+      auto compute_invA_ = [&compute_A_](Vector<ElemBasis>& sigma, Real& alpha, Real& beta, const Stellarator<Real,ORDER>& S, Vector<ElemBasis>& Bdotn, Real flux_tor, Real flux_pol, const Comm& comm) {
+        typename sctl::ParallelSolver<Real>::ParallelOp BIOp = [&S,&compute_A_](sctl::Vector<Real>* Ax, const sctl::Vector<Real>& x) {
+          (*Ax) = compute_A_(S, x);
+        };
+        const Long Nelem = S.NElem();
+        const Long Nnodes = ElemBasis::Size();
+
+        Vector<Real> rhs_(Nelem * Nnodes + S.Nsurf());
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            rhs_[i*Nnodes+j] = Bdotn[i][j];
+          }
+        }
+        if (S.Nsurf() >= 1) rhs_[Nelem * Nnodes + 0] = flux_tor;
+        if (S.Nsurf() >= 2) rhs_[Nelem * Nnodes + 1] = flux_pol;
+
+        Vector<Real> x_(Nelem * Nnodes + S.Nsurf());
+        x_ = 0;
+        ParallelSolver<Real> linear_solver(comm, true);
+        linear_solver(&x_, BIOp, rhs_, 1e-6, 100);
+
+        sigma.ReInit(Nelem);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            sigma[i][j] = x_[i*Nnodes+j];
+          }
+        }
+        alpha = (S.Nsurf() >= 1 ? x_[Nelem * Nnodes + 0] : 0);
+        beta  = (S.Nsurf() >= 2 ? x_[Nelem * Nnodes + 1] : 0);
+      };
       Comm comm = Comm::World();
       Profile::Enable(true);
 
@@ -4396,7 +4639,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
       Vector<Real> flux_tor(Nsurf), flux_pol(Nsurf);
       { // Init S, flux_tor, flux_pol, pressure
         Vector<Long> NtNp;
-        NtNp.PushBack(20);
+        NtNp.PushBack(30);
         NtNp.PushBack(4);
         S = Stellarator<Real,ORDER>(NtNp);
         flux_tor = 1;
@@ -4409,38 +4652,65 @@ template <class Real, Integer ORDER=10> class Stellarator {
       Vector<ElemBasis> normal, area_elem;
       compute_norm_area_elem(S, normal, area_elem);
 
-      Vector<ElemBasis> nu(Nnodes);
+      Vector<ElemBasis> nu(Nelem);
       { // Set nu
         nu = area_elem;
         //nu = 1;
+        //for (Long i = 0; i < Nelem; i++) {
+        //  for (Long j = 0; j < Nnodes; j++) {
+        //    Tensor<Real,true,COORD_DIM> x;
+        //    x(0) = S.Elem(i,0)[j];
+        //    x(1) = S.Elem(i,1)[j];
+        //    x(2) = S.Elem(i,2)[j];
+        //    nu[i][j] = x(2);
+        //  }
+        //}
       }
+      nu = nu * (1.0/sqrt(compute_inner_prod(area_elem, nu, nu)));
 
-      Vector<ElemBasis> B, n_dot_dBdn, dBdnu;
-      { // Set B, n_dot_dBdn
+      Vector<ElemBasis> B, nu_dBdn, nu_n_dot_dBdn;
+      { // Set B, nu_dBdn, nu_n_dot_dBdn
         Real alpha, beta;
         Vector<ElemBasis> sigma;
         compute_invA(sigma, alpha, beta, S, flux_tor[0], flux_pol[0], comm);
         B = compute_B(S, sigma, alpha, beta);
         Vector<ElemBasis> dB = compute_dB(S, sigma, alpha, beta);
-        n_dot_dBdn.ReInit(Nelem);
+        nu_dBdn.ReInit(Nelem * COORD_DIM);
+        nu_n_dot_dBdn.ReInit(Nelem);
         for (Long i = 0; i < Nelem; i++) {
           for (Long j = 0; j < Nnodes; j++) {
-            Real n_dot_dBdn_ = 0;
-            n_dot_dBdn_ += dB[(i*COORD_DIM+0)*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j];
-            n_dot_dBdn_ += dB[(i*COORD_DIM+1)*COORD_DIM+0][j] * normal[i*COORD_DIM+1][j] * normal[i*COORD_DIM+0][j];
-            n_dot_dBdn_ += dB[(i*COORD_DIM+2)*COORD_DIM+0][j] * normal[i*COORD_DIM+2][j] * normal[i*COORD_DIM+0][j];
-            n_dot_dBdn_ += dB[(i*COORD_DIM+0)*COORD_DIM+1][j] * normal[i*COORD_DIM+0][j] * normal[i*COORD_DIM+1][j];
-            n_dot_dBdn_ += dB[(i*COORD_DIM+1)*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j];
-            n_dot_dBdn_ += dB[(i*COORD_DIM+2)*COORD_DIM+1][j] * normal[i*COORD_DIM+2][j] * normal[i*COORD_DIM+1][j];
-            n_dot_dBdn_ += dB[(i*COORD_DIM+0)*COORD_DIM+2][j] * normal[i*COORD_DIM+0][j] * normal[i*COORD_DIM+2][j];
-            n_dot_dBdn_ += dB[(i*COORD_DIM+1)*COORD_DIM+2][j] * normal[i*COORD_DIM+1][j] * normal[i*COORD_DIM+2][j];
-            n_dot_dBdn_ += dB[(i*COORD_DIM+2)*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j];
-            n_dot_dBdn[i][j] = n_dot_dBdn_;
+            Real nu_dBdn_[COORD_DIM] = {0,0,0};
+            nu_dBdn_[0] -= dB[(i*COORD_DIM+0)*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j] * nu[i][j];
+            nu_dBdn_[0] -= dB[(i*COORD_DIM+0)*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j] * nu[i][j];
+            nu_dBdn_[0] -= dB[(i*COORD_DIM+0)*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j] * nu[i][j];
+            nu_dBdn_[1] -= dB[(i*COORD_DIM+1)*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j] * nu[i][j];
+            nu_dBdn_[1] -= dB[(i*COORD_DIM+1)*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j] * nu[i][j];
+            nu_dBdn_[1] -= dB[(i*COORD_DIM+1)*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j] * nu[i][j];
+            nu_dBdn_[2] -= dB[(i*COORD_DIM+2)*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j] * nu[i][j];
+            nu_dBdn_[2] -= dB[(i*COORD_DIM+2)*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j] * nu[i][j];
+            nu_dBdn_[2] -= dB[(i*COORD_DIM+2)*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j] * nu[i][j];
+            nu_dBdn[i*COORD_DIM+0][j] = nu_dBdn_[0];
+            nu_dBdn[i*COORD_DIM+1][j] = nu_dBdn_[1];
+            nu_dBdn[i*COORD_DIM+2][j] = nu_dBdn_[2];
+
+            Real nu_n_dot_dBdn_ = 0;
+            nu_n_dot_dBdn_ += nu_dBdn_[0] * normal[i*COORD_DIM+0][j];
+            nu_n_dot_dBdn_ += nu_dBdn_[1] * normal[i*COORD_DIM+1][j];
+            nu_n_dot_dBdn_ += nu_dBdn_[2] * normal[i*COORD_DIM+2][j];
+            nu_n_dot_dBdn[i][j] = nu_n_dot_dBdn_;
           }
         }
       }
-      { // Set dBdnu
-        Real eps = 1e-4;
+      { // Write VTU
+        VTUData vtu;
+        vtu.AddElems(S.GetElemList(), B, ORDER);
+        vtu.WriteVTK("B", comm);
+      }
+
+      Real dgdnu;
+      Vector<ElemBasis> dBdnu, n_dot_dBdnu;
+      { // Set dBdnu, n_dot_dBdnu, dgdnu (finite-difference approximation)
+        Real eps = 1e-3;
         Stellarator<Real,ORDER> S0 = S, S1 = S;
         for (Long i = 0; i < Nelem; i++) {
           for (Long j = 0; j < Nnodes; j++) {
@@ -4460,25 +4730,185 @@ template <class Real, Integer ORDER=10> class Stellarator {
         Vector<ElemBasis> sigma0, sigma1;
         compute_invA(sigma0, alpha0, beta0, S0, flux_tor[0], flux_pol[0], comm);
         compute_invA(sigma1, alpha1, beta1, S1, flux_tor[0], flux_pol[0], comm);
-        dBdnu = (compute_B(S1, sigma1, alpha1, beta1) - compute_B(S0, sigma0, alpha0, beta0)) * (1/eps);
+        Vector<ElemBasis> B0 = compute_B(S0, sigma0, alpha0, beta0);
+        Vector<ElemBasis> B1 = compute_B(S1, sigma1, alpha1, beta1);
+        dBdnu = (B1 - B0) * (1/eps);
+        dgdnu = (compute_g(S1,B1) - compute_g(S0,B0)) * (1/eps);
+
+        n_dot_dBdnu.ReInit(Nelem);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            Real n_dot_dBdnu_ = 0;
+            n_dot_dBdnu_ += normal[i*COORD_DIM+0][j] * dBdnu[i*COORD_DIM+0][j];
+            n_dot_dBdnu_ += normal[i*COORD_DIM+1][j] * dBdnu[i*COORD_DIM+1][j];
+            n_dot_dBdnu_ += normal[i*COORD_DIM+2][j] * dBdnu[i*COORD_DIM+2][j];
+            n_dot_dBdnu[i][j] = n_dot_dBdnu_;
+          }
+        }
       }
 
+      Vector<ElemBasis> B_dot_gradnu, nu_surfdivB, surfdivBnu;
+      { // Set B_dot_gradnu
+        Vector<ElemBasis> gradnu = compute_grad(S, nu);
+        B_dot_gradnu.ReInit(Nelem);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            Real B_dot_gradnu_ = 0;
+            B_dot_gradnu_ += B[i*COORD_DIM+0][j] * gradnu[i*COORD_DIM+0][j];
+            B_dot_gradnu_ += B[i*COORD_DIM+1][j] * gradnu[i*COORD_DIM+1][j];
+            B_dot_gradnu_ += B[i*COORD_DIM+2][j] * gradnu[i*COORD_DIM+2][j];
+            B_dot_gradnu[i][j] = B_dot_gradnu_;
+          }
+        }
+      }
+      { // Set nu_surfdivB
+        Vector<ElemBasis> surfdivB = compute_surfdiv(S, B);
+        nu_surfdivB.ReInit(Nelem);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            nu_surfdivB[i][j] = nu[i][j] * surfdivB[i][j];
+          }
+        }
+      }
+      { // Set surfdivBnu
+        Vector<ElemBasis> Bnu(Nelem*COORD_DIM);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            Bnu[i*COORD_DIM+0][j] = B[i*COORD_DIM+0][j] * nu[i][j];
+            Bnu[i*COORD_DIM+1][j] = B[i*COORD_DIM+1][j] * nu[i][j];
+            Bnu[i*COORD_DIM+2][j] = B[i*COORD_DIM+2][j] * nu[i][j];
+          }
+        }
+        surfdivBnu = compute_surfdiv(S, Bnu);
+      }
+
+      // nu_surfdivB == -nu_n_dot_dBdn
+      // B_dot_gradnu == n_dot_dBdnu
+      // surfdivBnu == B_dot_gradnu - nu_n_dot_dBdn
+
+      Vector<ElemBasis> dBdnu_;
+      { // Compute dBdnu_
+        Real alpha, beta;
+        Real flux_tor, flux_pol;
+        { // Set flux_tor, flux_pol
+          Vector<ElemBasis> nxB_nu(Nelem * COORD_DIM);
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              nxB_nu[i*COORD_DIM+0][j] = (B[i*COORD_DIM+1][j] * normal[i*COORD_DIM+2][j] - B[i*COORD_DIM+2][j] * normal[i*COORD_DIM+1][j]) * nu[i][j];
+              nxB_nu[i*COORD_DIM+1][j] = (B[i*COORD_DIM+2][j] * normal[i*COORD_DIM+0][j] - B[i*COORD_DIM+0][j] * normal[i*COORD_DIM+2][j]) * nu[i][j];
+              nxB_nu[i*COORD_DIM+2][j] = (B[i*COORD_DIM+0][j] * normal[i*COORD_DIM+1][j] - B[i*COORD_DIM+1][j] * normal[i*COORD_DIM+0][j]) * nu[i][j];
+            }
+          }
+
+          Vector<Real> circ_pol(S.Nsurf()), circ_tor(S.Nsurf());
+          { // compute circ
+            Vector<ElemBasis> dX;
+            ElemBasis::Grad(dX, S.GetElemList().ElemVector());
+            const auto& quad_wts = ElemBasis::QuadWts();
+            for (Long k = 0; k < S.Nsurf(); k++) {
+              circ_pol[k] = 0;
+              circ_tor[k] = 0;
+              Long Ndsp = S.ElemDsp(k);
+              for (Long i = 0; i < S.NTor(k)*S.NPol(k); i++) {
+                for (Long j = 0; j < Nnodes; j++) {
+                  circ_pol[k] += nxB_nu[(Ndsp+i)*COORD_DIM+0][j] * dX[(Ndsp+i)*COORD_DIM*2+1][j] * quad_wts[j] / S.NTor(k);
+                  circ_pol[k] += nxB_nu[(Ndsp+i)*COORD_DIM+1][j] * dX[(Ndsp+i)*COORD_DIM*2+3][j] * quad_wts[j] / S.NTor(k);
+                  circ_pol[k] += nxB_nu[(Ndsp+i)*COORD_DIM+2][j] * dX[(Ndsp+i)*COORD_DIM*2+5][j] * quad_wts[j] / S.NTor(k);
+
+                  circ_tor[k] += nxB_nu[(Ndsp+i)*COORD_DIM+0][j] * dX[(Ndsp+i)*COORD_DIM*2+0][j] * quad_wts[j] / S.NPol(k);
+                  circ_tor[k] += nxB_nu[(Ndsp+i)*COORD_DIM+1][j] * dX[(Ndsp+i)*COORD_DIM*2+2][j] * quad_wts[j] / S.NPol(k);
+                  circ_tor[k] += nxB_nu[(Ndsp+i)*COORD_DIM+2][j] * dX[(Ndsp+i)*COORD_DIM*2+4][j] * quad_wts[j] / S.NPol(k);
+                }
+              }
+            }
+          }
+          if (S.Nsurf() == 1) {
+            flux_tor = circ_pol[0];
+            flux_pol = 0;
+          } else if (S.Nsurf() == 2) {
+            flux_tor = circ_pol[1] - circ_pol[0];
+            flux_pol = circ_tor[0] - circ_tor[1];
+          } else {
+            SCTL_ASSERT(false);
+          }
+        }
+        Vector<ElemBasis> sigma, Bdotn = B_dot_gradnu - nu_n_dot_dBdn;
+        compute_invA_(sigma, alpha, beta, S, Bdotn, flux_tor, flux_pol, comm);
+        dBdnu_ = compute_B(S, sigma, alpha, beta) + nu_dBdn;
+      }
       { // Write VTU
         VTUData vtu;
-        vtu.AddElems(S.GetElemList(), B, ORDER);
-        vtu.WriteVTK("B", comm);
+        vtu.AddElems(S.GetElemList(), dBdnu, ORDER);
+        vtu.WriteVTK("dBdnu", comm);
       }
       { // Write VTU
         VTUData vtu;
-        vtu.AddElems(S.GetElemList(), n_dot_dBdn, ORDER);
-        vtu.WriteVTK("n_dot_dBdn", comm);
+        vtu.AddElems(S.GetElemList(), dBdnu_, ORDER);
+        vtu.WriteVTK("dBdnu_", comm);
       }
       { // Write VTU
         VTUData vtu;
-        vtu.AddElems(S.GetElemList(), dBdnu, ORDER);
-        vtu.WriteVTK("dBdnu", comm);
+        vtu.AddElems(S.GetElemList(), dBdnu_ - dBdnu, ORDER);
+        vtu.WriteVTK("err", comm);
       }
 
+
+      Real dgdnu0, dgdnu1, dgdnu2;
+      { // Set dgdnu0 = \int_{Gamma} (B^2 - p) B . B'
+        Vector<ElemBasis> dB = dBdnu - nu_dBdn;
+        Vector<ElemBasis> B2_p(Nelem), B_dot_dB(Nelem);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            B2_p[i][j] = 0;
+            B2_p[i][j] += B[i*COORD_DIM+0][j]*B[i*COORD_DIM+0][j];
+            B2_p[i][j] += B[i*COORD_DIM+1][j]*B[i*COORD_DIM+1][j];
+            B2_p[i][j] += B[i*COORD_DIM+2][j]*B[i*COORD_DIM+2][j];
+
+            B_dot_dB[i][j] = 0;
+            B_dot_dB[i][j] += B[i*COORD_DIM+0][j] * dB[i*COORD_DIM+0][j];
+            B_dot_dB[i][j] += B[i*COORD_DIM+1][j] * dB[i*COORD_DIM+1][j];
+            B_dot_dB[i][j] += B[i*COORD_DIM+2][j] * dB[i*COORD_DIM+2][j];
+          }
+        }
+        dgdnu0 = compute_inner_prod(area_elem, B2_p, B_dot_dB);
+      }
+      { // Set dgdnu1 = \int_{Gamma} (B^2-p) B . nu_dBdn
+        Vector<ElemBasis> dB = nu_dBdn;
+        Vector<ElemBasis> B2_p(Nelem), B_dot_dB(Nelem);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            B2_p[i][j] = 0;
+            B2_p[i][j] += B[i*COORD_DIM+0][j]*B[i*COORD_DIM+0][j];
+            B2_p[i][j] += B[i*COORD_DIM+1][j]*B[i*COORD_DIM+1][j];
+            B2_p[i][j] += B[i*COORD_DIM+2][j]*B[i*COORD_DIM+2][j];
+
+            B_dot_dB[i][j] = 0;
+            B_dot_dB[i][j] += B[i*COORD_DIM+0][j] * dB[i*COORD_DIM+0][j];
+            B_dot_dB[i][j] += B[i*COORD_DIM+1][j] * dB[i*COORD_DIM+1][j];
+            B_dot_dB[i][j] += B[i*COORD_DIM+2][j] * dB[i*COORD_DIM+2][j];
+          }
+        }
+        dgdnu1 = compute_inner_prod(area_elem, B2_p, B_dot_dB);
+      }
+      { // Set dgdnu2 = \int_{Gamma} 2H(B^2-p)^2 \nu
+        Vector<ElemBasis> H = compute_H(S.GetElemList(), normal);
+        Vector<ElemBasis> H_B2_p_2(Nelem);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            Real B2_p = 0;
+            B2_p += B[i*COORD_DIM+0][j]*B[i*COORD_DIM+0][j];
+            B2_p += B[i*COORD_DIM+1][j]*B[i*COORD_DIM+1][j];
+            B2_p += B[i*COORD_DIM+2][j]*B[i*COORD_DIM+2][j];
+            H_B2_p_2[i][j] = H[i][j] * B2_p*B2_p;
+          }
+        }
+        dgdnu2 = 0.5 * compute_inner_prod(area_elem,H_B2_p_2, nu);
+      }
+
+      std::cout<<dgdnu0<<' '<<dgdnu1<<' '<<dgdnu2<<' '<<dgdnu0+dgdnu1+dgdnu2<<'\n';
+      std::cout<<dgdnu<<'\n';
+
+
 #if 0
       Comm comm = Comm::World();
       Profile::Enable(true);
@@ -5387,10 +5817,11 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       }
 
       Real fx;
-      if (0) {
+      if (1) {
         LBFGSpp::LBFGSParam<Real> param;
         param.max_iterations = 100;
         param.epsilon = 1e-8;
+        param.m = 20;
         LBFGSpp::LBFGSSolver<Real> solver(param);
         Integer niter = solver.minimize(mhd_equilib, x, fx);
       } else {