.

include/sctl/boundary_quadrature.hpp

@@ -2893,20 +2893,86 @@ template <class Real, Integer ORDER=10> class Stellarator {
         return x;
       };
 
-      auto compute_dg_dsigma = [&S, &normal, &compute_dot_prod](const Vector<ElemBasis>& B) { // dg_dsigma = \int 2 B \cdot (\nabla G + n/2)
+      /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+      /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+      const auto pressure = area_elem;
+      auto compute_gvec = [&S,&pressure] (const Vector<ElemBasis>& B) {
+        const Long Nelem = S.GetElemList().NElem();
+        const Long Nnodes = ElemBasis::Size();
+        SCTL_ASSERT(B.Dim() == Nelem * COORD_DIM);
+        SCTL_ASSERT(pressure.Dim() == Nelem);
+
+        Vector<ElemBasis> gvec(Nelem);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            Real B2 = 0;
+            B2 += B[i*COORD_DIM+0][j] * B[i*COORD_DIM+0][j];
+            B2 += B[i*COORD_DIM+1][j] * B[i*COORD_DIM+1][j];
+            B2 += B[i*COORD_DIM+2][j] * B[i*COORD_DIM+2][j];
+            gvec[i][j] = (B2*0.5 - pressure[i][j]) * (B2*0.5 - pressure[i][j]);
+          }
+        }
+        return gvec;
+      };
+      auto compute_dgdB = [&S,&pressure,&compute_inner_prod] (const Vector<ElemBasis>& B) {
+        const Long Nelem = S.GetElemList().NElem();
+        const Long Nnodes = ElemBasis::Size();
+        SCTL_ASSERT(B.Dim() == Nelem * COORD_DIM);
+        SCTL_ASSERT(pressure.Dim() == Nelem);
+
+        Vector<ElemBasis> dgdB(Nelem*COORD_DIM);
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            Real B2 = 0;
+            B2 += B[i*COORD_DIM+0][j] * B[i*COORD_DIM+0][j];
+            B2 += B[i*COORD_DIM+1][j] * B[i*COORD_DIM+1][j];
+            B2 += B[i*COORD_DIM+2][j] * B[i*COORD_DIM+2][j];
+            dgdB[i*COORD_DIM+0][j] = 2 * (B2*0.5 - pressure[i][j]) * B[i*COORD_DIM+0][j];
+            dgdB[i*COORD_DIM+1][j] = 2 * (B2*0.5 - pressure[i][j]) * B[i*COORD_DIM+1][j];
+            dgdB[i*COORD_DIM+2][j] = 2 * (B2*0.5 - pressure[i][j]) * B[i*COORD_DIM+2][j];
+          }
+        }
+        return dgdB;
+      };
+
+      Real flux = 1.0, alpha;
+      Vector<ElemBasis> sigma(S.GetElemList().NElem());
+      compute_invA(sigma, alpha, flux);
+      Vector<ElemBasis> B = compute_half_n_plus_dG(sigma) + compute_B0(alpha);
+      Real g = compute_inner_prod(compute_gvec(B), area_elem*0+1);
+      std::cout<<"g = "<<g<<'\n';
+
+      { // Write VTU
+        VTUData vtu;
+        vtu.AddElems(S.GetElemList(), sigma, ORDER);
+        vtu.WriteVTK("sigma", comm);
+      }
+      { // Write VTU
+        VTUData vtu;
+        vtu.AddElems(S.GetElemList(), B, ORDER);
+        vtu.WriteVTK("B", comm);
+      }
+
+      /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+      /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+      auto compute_dg_dsigma = [&S,&compute_dgdB,&normal,&compute_dot_prod](const Vector<ElemBasis>& B) { // dg_dsigma = \int 2 B \cdot (\nabla G + n/2)
         Vector<ElemBasis> B_dot_gradG;
-        S.quadrature_dUxF.Eval(B_dot_gradG, S.GetElemList(), B, S.Laplace_dUxF);
-        return B_dot_gradG * (-2.0) + compute_dot_prod(B,normal);
+        Vector<ElemBasis> dgdB = compute_dgdB(B);
+        S.quadrature_dUxF.Eval(B_dot_gradG, S.GetElemList(), dgdB, S.Laplace_dUxF);
+        return B_dot_gradG * (-1.0) + compute_dot_prod(dgdB,normal) * 0.5;
       };
-      auto compute_dg_dalpha = [&S,&compute_B0,&compute_inner_prod] (const Vector<ElemBasis>& B) {
+      auto compute_dg_dalpha = [&S,&compute_dgdB,&compute_B0,&compute_inner_prod] (const Vector<ElemBasis>& B) {
         auto dB_dalpha = compute_B0(1);
-        return 2*compute_inner_prod(B,dB_dalpha);
+        Vector<ElemBasis> dgdB = compute_dgdB(B);
+        return compute_inner_prod(dgdB,dB_dalpha);
       };
-
-      auto compute_dg_dnu = [&S,&comm,&normal,&compute_inner_prod,&area_elem,&compute_dB0](const Vector<ElemBasis>& sigma, Real alpha, const Vector<ElemBasis>& B) { // dg_dnu = (B*B) 2H - (2 B) \cdot (n \cdnot nabla) \nabla G[sigma] + (2 B) \alpha dB0_dnu \hat{\theta} + sigma (\nabla D)^T [2 B] + (2H) sigma (\nabla G)^T [2 B]
+      auto compute_dg_dnu = [&S,&compute_gvec,&compute_dgdB,&comm,&normal,&compute_inner_prod,&area_elem,&compute_dB0](const Vector<ElemBasis>& sigma, Real alpha, const Vector<ElemBasis>& B) { // dg_dnu = (B*B) 2H - (2 B) \cdot (n \cdnot nabla) \nabla G[sigma] + (2 B) \alpha dB0_dnu \hat{\theta} + sigma (\nabla D)^T [2 B] + (2H) sigma (\nabla G)^T [2 B]
         const Long Nelem = S.GetElemList().NElem();
         const Long Nnodes = ElemBasis::Size();
-        Vector<ElemBasis> v = B * 2.0;
+        Vector<ElemBasis> gvec = compute_gvec(B);
+        Vector<ElemBasis> v = compute_dgdB(B);
 
         Vector<ElemBasis> dg_dnu0(Nelem), dg_dnu1(Nelem), dg_dnu2(Nelem), dg_dnu3(Nelem), dg_dnu4(Nelem);
         dg_dnu0 = 0;
@@ -2959,9 +3025,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
         for (Long i = 0; i < Nelem; i++) {
           for (Long j = 0; j < Nnodes; j++) {
             dg_dnu0[i][j] = 0;
-            dg_dnu0[i][j] += B[i*COORD_DIM+0][j] * B[i*COORD_DIM+0][j] * (2.0*H[i][j]);
-            dg_dnu0[i][j] += B[i*COORD_DIM+1][j] * B[i*COORD_DIM+1][j] * (2.0*H[i][j]);
-            dg_dnu0[i][j] += B[i*COORD_DIM+2][j] * B[i*COORD_DIM+2][j] * (2.0*H[i][j]);
+            dg_dnu0[i][j] += gvec[i][j] * (2.0*H[i][j]);
           }
         }
 
@@ -3036,29 +3100,8 @@ template <class Real, Integer ORDER=10> class Stellarator {
       /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
       /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-      Real flux = 1.0, alpha;
-      Vector<ElemBasis> sigma(S.GetElemList().NElem());
-      compute_invA(sigma, alpha, flux);
-      Vector<ElemBasis> B = compute_half_n_plus_dG(sigma) + compute_B0(alpha);
-      Real g = compute_inner_prod(B, B);
-      std::cout<<"g = "<<g<<'\n';
-
-      { // Write VTU
-        VTUData vtu;
-        vtu.AddElems(S.GetElemList(), sigma, ORDER);
-        vtu.WriteVTK("sigma", comm);
-      }
-      { // Write VTU
-        VTUData vtu;
-        vtu.AddElems(S.GetElemList(), B, ORDER);
-        vtu.WriteVTK("B", comm);
-      }
-
-      /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-      /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
       if (0) { // test dg_dnu
-        auto compute_g = [&S,&comm,&normal,&area_elem,&sigma,&alpha,&compute_norm_area_elem,&compute_B0,&compute_inner_prod](const Vector<ElemBasis>& nu, Real eps) {
+        auto compute_g = [&compute_gvec,&S,&comm,&normal,&area_elem,&sigma,&alpha,&compute_norm_area_elem,&compute_B0,&compute_inner_prod](const Vector<ElemBasis>& nu, Real eps) {
           const Long Nelem = S.GetElemList().NElem();
           const Long Nnodes = ElemBasis::Size();
 
@@ -3090,7 +3133,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           quadrature_FxdU.template Setup<ElemBasis>(S.GetElemList(), Xt, S.Laplace_FxdU, order_singular, order_direct, -1, comm);
           quadrature_FxdU.Eval(B1, S.GetElemList(), sigma, S.Laplace_FxdU);
 
-          Real g = compute_inner_prod(B0+B1, B0+B1);
+          Real g = compute_inner_prod(compute_gvec(B0+B1),area_elem*0+1);
 
           for (Long i = 0; i < Nelem; i++) {
             for (Long j = 0; j < Nnodes; j++) {
@@ -5295,7 +5338,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           dg_dnu = f;
         }
 
-        auto compute_g = [&sigma,&alpha,&S,&area_elem,&normal,&compute_norm_area_elem,&compute_invA,&compute_half_n_plus_dG,&compute_B0,&compute_inner_prod,&comm] (const Vector<ElemBasis>& nu, Real eps) {
+        auto compute_g = [&compute_gvec,&sigma,&alpha,&S,&area_elem,&normal,&compute_norm_area_elem,&compute_invA,&compute_half_n_plus_dG,&compute_B0,&compute_inner_prod,&comm] (const Vector<ElemBasis>& nu, Real eps) {
           const Long Nelem = S.GetElemList().NElem();
           const Long Nnodes = ElemBasis::Size();
 
@@ -5320,7 +5363,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           Vector<ElemBasis> sigma;
           compute_invA(sigma, alpha, flux);
           Vector<ElemBasis> B = compute_half_n_plus_dG(sigma) + compute_B0(alpha);
-          Real g = compute_inner_prod(B, B);
+          Real g = compute_inner_prod(compute_gvec(B), area_elem*0+1);
 
           for (Long i = 0; i < Nelem; i++) {
             for (Long j = 0; j < Nnodes; j++) {