.

include/sctl/boundary_quadrature.hpp

@@ -2716,6 +2716,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
           Ax2 = (S.Nsurf() >= 2 ? compute_inner_prod(area_elem, compute_dot_prod(Bp0, normal), x0) : 0);
         }
 
+        // TODO: precompute A21adj, A22adj
         auto compute_A21adj = [&S,&normal,&area_elem] (bool toroidal_flux) {
           const Long Nelem = S.NElem();
           const Long Nnodes = ElemBasis::Size();
@@ -2776,7 +2777,6 @@ template <class Real, Integer ORDER=10> class Stellarator {
         };
         if (S.Nsurf() >= 1) Ax0 += compute_A21adj( true) * x1;
         if (S.Nsurf() >= 2) Ax0 += compute_A21adj(false) * x2;
-
         if (S.Nsurf() == 1) { // Add flux part of Ax1, Ax2
           Real flux_tor, flux_pol;
           compute_flux(flux_tor, flux_pol, Bt0, normal);
@@ -2816,6 +2816,65 @@ template <class Real, Integer ORDER=10> class Stellarator {
         return x;
       };
 
+      if (0) { // Check u_t A v == v_t Aadj u
+        auto pack = [&S](Vector<Real>& x, const Vector<ElemBasis>& x0, Real x1, Real x2) {
+          const Long Nelem = S.NElem();
+          const Long Nnodes = ElemBasis::Size();
+          x.ReInit(Nelem*Nnodes+S.Nsurf());
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              x[i*Nnodes+j] = x0[i][j];
+            }
+          }
+          if (S.Nsurf() >= 1) x[Nelem*Nnodes+0] = x1;
+          if (S.Nsurf() >= 2) x[Nelem*Nnodes+1] = x2;
+        };
+        auto unpack = [&S](Vector<ElemBasis>& x0, Real& x1, Real& x2, const Vector<Real>& x) {
+          const Long Nelem = S.NElem();
+          const Long Nnodes = ElemBasis::Size();
+          x0.ReInit(Nelem);
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              x0[i][j] = x[i*Nnodes+j];
+            }
+          }
+          if (S.Nsurf() >= 1) x1 = x[Nelem*Nnodes+0];
+          if (S.Nsurf() >= 2) x2 = x[Nelem*Nnodes+1];
+        };
+        const Long Nelem = S.NElem();
+        const Long Nnodes = ElemBasis::Size();
+
+        Vector<ElemBasis> normal, area_elem;
+        compute_norm_area_elem(S, normal, area_elem);
+
+        Vector<Real> u, v;
+        Vector<ElemBasis> u0 = area_elem*0; Real u1 = 3.141, u2 = 0.4142;
+        Vector<ElemBasis> v0 = area_elem*0; Real v1 = 1.645, v2 = 3.6055;
+        { // Set u0, v0
+          auto X = S.GetElemList().ElemVector();
+          for (Long i = 0; i < Nelem; i++) {
+            for (Long j = 0; j < Nnodes; j++) {
+              u0[i][j] = X[i*COORD_DIM+0][j] + X[i*COORD_DIM+1][j] + X[i*COORD_DIM+2][j];
+              v0[i][j] = X[i*COORD_DIM+0][j] + X[i*COORD_DIM+1][j] + X[i*COORD_DIM+2][j];
+            }
+          }
+        }
+        pack(u,u0,u1,u2);
+        pack(v,v0,v1,v2);
+
+        Vector<Real> Av = compute_A(v);
+        Vector<Real> uA = compute_Aadj(u);
+
+        Vector<ElemBasis> Av0; Real Av1, Av2;
+        Vector<ElemBasis> uA0; Real uA1, uA2;
+        unpack(Av0, Av1, Av2, Av);
+        unpack(uA0, uA1, uA2, uA);
+
+        std::cout << compute_inner_prod(area_elem, u0, Av0) + u1*Av1 + (S.Nsurf()>=2?u2*Av2:0) << '\n';
+        std::cout << compute_inner_prod(area_elem, uA0, v0) + uA1*v1 + (S.Nsurf()>=2?uA2*v2:0) << '\n';
+        exit(0);
+      }
+
       /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
       /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
@@ -3618,199 +3677,6 @@ 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);
       };
-      auto compute_Aadj = [&S,&Bt0,&Bp0,&compute_flux] (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);
-
-        Vector<ElemBasis> x0(Nelem);
-        for (Long i = 0; i < Nelem; i++) {
-          for (Long j = 0; j < Nnodes; j++) {
-            x0[i][j] = x[i*Nnodes+j];
-          }
-        }
-        Real x1 = (S.Nsurf() >= 1 ? x[Nelem*Nnodes + 0] : 0);
-        Real x2 = (S.Nsurf() >= 2 ? x[Nelem*Nnodes + 1] : 0);
-
-        Vector<ElemBasis> Ax0;
-        Real Ax1, Ax2;
-        { // Set Ax0, Ax1, Ax2
-          Vector<ElemBasis> x0_n(Nelem*COORD_DIM);
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              x0_n[i*COORD_DIM+0][j] = x0[i][j] * normal[i*COORD_DIM+0][j];
-              x0_n[i*COORD_DIM+1][j] = x0[i][j] * normal[i*COORD_DIM+1][j];
-              x0_n[i*COORD_DIM+2][j] = x0[i][j] * normal[i*COORD_DIM+2][j];
-            }
-          }
-          EvalQuadrature(Ax0, S.quadrature_dUxF, S, x0_n, S.Laplace_dUxF);
-          Ax0 = x0*(-0.5) - Ax0;
-
-          Ax1 = (S.Nsurf() >= 1 ? compute_inner_prod(area_elem, compute_dot_prod(Bt0, normal), x0) : 0);
-          Ax2 = (S.Nsurf() >= 2 ? compute_inner_prod(area_elem, compute_dot_prod(Bp0, normal), x0) : 0);
-        }
-
-        auto compute_A21adj = [&S,&normal,&area_elem] (bool toroidal_flux) {
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-          Vector<ElemBasis> density(Nelem * COORD_DIM);
-          { // Set density
-            Real scal[2];
-            if (S.Nsurf() == 1) {
-              SCTL_ASSERT(toroidal_flux == true);
-              scal[0] = 1.0 / S.NTor(0);
-              scal[1] = 0;
-            } else if (S.Nsurf() == 2) {
-              if (toroidal_flux == true) {
-                scal[0] = -1.0 / S.NTor(0);
-                scal[1] = 1.0 / S.NTor(1);
-              } else {
-                scal[0] = 1.0 / S.NPol(0);
-                scal[1] = -1.0 / S.NPol(1);
-              }
-            } else {
-              SCTL_ASSERT(false);
-            }
-
-            Vector<ElemBasis> dX;
-            ElemBasis::Grad(dX, S.GetElemList().ElemVector());
-            for (Long k = 0; k < S.Nsurf(); k++) {
-              for (Long i_ = 0; i_ < S.NTor(k)*S.NPol(k); i_++) {
-                Long i = S.ElemDsp(k) + i_;
-                for (Long j = 0; j < Nnodes; j++) {
-                  Real s = scal[k] / area_elem[i][j];
-                  density[i*COORD_DIM+0][j] = dX[i*COORD_DIM*2+0+(toroidal_flux?1:0)][j] * s;
-                  density[i*COORD_DIM+1][j] = dX[i*COORD_DIM*2+2+(toroidal_flux?1:0)][j] * s;
-                  density[i*COORD_DIM+2][j] = dX[i*COORD_DIM*2+4+(toroidal_flux?1:0)][j] * s;
-                }
-              }
-            }
-          }
-
-          Vector<ElemBasis> Gdensity, nxGdensity(Nelem * COORD_DIM), A21adj;
-          EvalQuadrature(Gdensity, S.quadrature_FxU, S, density, S.Laplace_FxU);
-          for (Long i = 0; i < Nelem; i++) { // Set nxGdensity
-            for (Long j = 0; j < Nnodes; j++) {
-              Tensor<Real,true,COORD_DIM> Gdensity_, n;
-              Gdensity_(0) = Gdensity[i*COORD_DIM+0][j];
-              Gdensity_(1) = Gdensity[i*COORD_DIM+1][j];
-              Gdensity_(2) = Gdensity[i*COORD_DIM+2][j];
-
-              n(0) = normal[i*COORD_DIM+0][j];
-              n(1) = normal[i*COORD_DIM+1][j];
-              n(2) = normal[i*COORD_DIM+2][j];
-
-              nxGdensity[i*COORD_DIM+0][j] = n(1) * Gdensity_(2) - n(2) * Gdensity_(1);
-              nxGdensity[i*COORD_DIM+1][j] = n(2) * Gdensity_(0) - n(0) * Gdensity_(2);
-              nxGdensity[i*COORD_DIM+2][j] = n(0) * Gdensity_(1) - n(1) * Gdensity_(0);
-            }
-          }
-          EvalQuadrature(A21adj, S.quadrature_dUxF, S, nxGdensity, S.Laplace_dUxF);
-          return A21adj;
-        };
-        if (S.Nsurf() >= 1) Ax0 += compute_A21adj( true) * x1;
-        if (S.Nsurf() >= 2) Ax0 += compute_A21adj(false) * x2;
-
-        if (S.Nsurf() == 1) { // Add flux part of Ax1, Ax2
-          Real flux_tor, flux_pol;
-          compute_flux(flux_tor, flux_pol, Bt0, normal);
-          Ax1 += flux_tor * x1;
-          Ax2 += 0;
-        } else if (S.Nsurf() == 2) {
-          Real flux_tor, flux_pol;
-          compute_flux(flux_tor, flux_pol, Bt0, normal);
-          Ax1 += flux_tor * x1 + flux_pol * x2;
-          compute_flux(flux_tor, flux_pol, Bp0, normal);
-          Ax2 += flux_tor * x1 + flux_pol * x2;
-        } else {
-          SCTL_ASSERT(false);
-        }
-
-        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] = Ax0[i][j];
-          }
-        }
-        if (S.Nsurf() >= 1) Ax[Nelem*Nnodes + 0] = Ax1;
-        if (S.Nsurf() >= 2) Ax[Nelem*Nnodes + 1] = Ax2;
-        return Ax;
-      };
-      auto compute_invAadj = [&S,&comm,&compute_Aadj] (Vector<Real>& b) {
-        typename sctl::ParallelSolver<Real>::ParallelOp BIOp = [&compute_Aadj](sctl::Vector<Real>* Ax, const sctl::Vector<Real>& x) {
-          (*Ax) = compute_Aadj(x);
-        };
-        const Long Nelem = S.NElem();
-        const Long Nnodes = ElemBasis::Size();
-
-        Vector<Real> x(b.Dim());
-        x = 0;
-        ParallelSolver<Real> linear_solver(comm, true);
-        linear_solver(&x, BIOp, b, 1e-8, 100);
-        return x;
-      };
-
-      if (0) { // Check u_t A v == v_t Aadj u
-        auto pack = [&S](Vector<Real>& x, const Vector<ElemBasis>& x0, Real x1, Real x2) {
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-          x.ReInit(Nelem*Nnodes+S.Nsurf());
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              x[i*Nnodes+j] = x0[i][j];
-            }
-          }
-          if (S.Nsurf() >= 1) x[Nelem*Nnodes+0] = x1;
-          if (S.Nsurf() >= 2) x[Nelem*Nnodes+1] = x2;
-        };
-        auto unpack = [&S](Vector<ElemBasis>& x0, Real& x1, Real& x2, const Vector<Real>& x) {
-          const Long Nelem = S.NElem();
-          const Long Nnodes = ElemBasis::Size();
-          x0.ReInit(Nelem);
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              x0[i][j] = x[i*Nnodes+j];
-            }
-          }
-          if (S.Nsurf() >= 1) x1 = x[Nelem*Nnodes+0];
-          if (S.Nsurf() >= 2) x2 = x[Nelem*Nnodes+1];
-        };
-        const Long Nelem = S.NElem();
-        const Long Nnodes = ElemBasis::Size();
-
-        Vector<ElemBasis> normal, area_elem;
-        compute_norm_area_elem(S, normal, area_elem);
-
-        Vector<Real> u, v;
-        Vector<ElemBasis> u0 = area_elem*0; Real u1 = 3.141, u2 = 0.4142;
-        Vector<ElemBasis> v0 = area_elem*0; Real v1 = 1.645, v2 = 3.6055;
-        { // Set u0, v0
-          auto X = S.GetElemList().ElemVector();
-          for (Long i = 0; i < Nelem; i++) {
-            for (Long j = 0; j < Nnodes; j++) {
-              u0[i][j] = X[i*COORD_DIM+0][j] + X[i*COORD_DIM+1][j] + X[i*COORD_DIM+2][j];
-              v0[i][j] = X[i*COORD_DIM+0][j] + X[i*COORD_DIM+1][j] + X[i*COORD_DIM+2][j];
-            }
-          }
-        }
-        pack(u,u0,u1,u2);
-        pack(v,v0,v1,v2);
-
-        Vector<Real> Av = compute_A(v);
-        Vector<Real> uA = compute_Aadj(u);
-
-        Vector<ElemBasis> Av0; Real Av1, Av2;
-        Vector<ElemBasis> uA0; Real uA1, uA2;
-        unpack(Av0, Av1, Av2, Av);
-        unpack(uA0, uA1, uA2, uA);
-
-        std::cout << compute_inner_prod(area_elem, u0, Av0) + u1*Av1 + (S.Nsurf()>=2?u2*Av2:0) << '\n';
-        std::cout << compute_inner_prod(area_elem, uA0, v0) + uA1*v1 + (S.Nsurf()>=2?uA2*v2:0) << '\n';
-        exit(0);
-      }
 
       Vector<ElemBasis> dg_dnu = compute_gradient(S, pressure, flux_tor, flux_pol);
       { // Write VTU