.

include/sctl/boundary_quadrature.hpp

@@ -2370,7 +2370,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
         Jp__ = grid2cheb(Jp_, Nt, Np, S.NTor(k), S.NPol(k));
       }
     }
-    static void compute_norm_area_elem(const Stellarator<Real,10>& S, Vector<ElemBasis>& normal, Vector<ElemBasis>& area_elem){ // Set normal, area_elem
+    static void compute_norm_area_elem(const Stellarator<Real,ORDER>& S, Vector<ElemBasis>& normal, Vector<ElemBasis>& area_elem){ // Set normal, area_elem
       const Vector<ElemBasis>& X = S.GetElemList().ElemVector();
       const Long Nelem = X.Dim() / COORD_DIM;
       const Long Nnodes = ElemBasis::Size();
@@ -4370,6 +4370,116 @@ template <class Real, Integer ORDER=10> class Stellarator {
     }
 
     static void test_askham() {
+      auto Setup = [] (Stellarator<Real,ORDER>& S, const Comm& comm) { // Set quadratures, Bt0, Bp0, ...
+        SetupQuadrature(S.quadrature_dBS  , S, S.BiotSavartGrad, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
+        SetupQuadrature(S.quadrature_BS   , S, S.BiotSavart    , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
+        SetupQuadrature(S.quadrature_FxU  , S, S.Laplace_FxU   , order_singular, order_direct, -1.0, comm);
+        SetupQuadrature(S.quadrature_FxdU , S, S.Laplace_FxdU  , order_singular, order_direct, -1.0, comm);
+        SetupQuadrature(S.quadrature_dUxF , S, S.Laplace_dUxF  , order_singular, order_direct, -1.0, comm);
+        SetupQuadrature(S.quadrature_dUxD , S, S.Laplace_dUxD  , order_singular, order_direct, -1.0, comm,  0.01 * pow<-2,Real>(ORDER));
+        SetupQuadrature(S.quadrature_Fxd2U, S, S.Laplace_Fxd2U , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER));
+
+        { // Set Bt0, Bp0, dBt0, dBp0
+          Vector<ElemBasis> Jt, Jp;
+          compute_harmonic_vector_potentials(Jt, Jp, S);
+          EvalQuadrature(S.Bt0 , S.quadrature_BS , S, Jp, S.BiotSavart);
+          EvalQuadrature(S.Bp0 , S.quadrature_BS , S, Jt, S.BiotSavart);
+          EvalQuadrature(S.dBt0, S.quadrature_dBS, S, Jp, S.BiotSavartGrad);
+          EvalQuadrature(S.dBp0, S.quadrature_dBS, S, Jt, S.BiotSavartGrad);
+        }
+      };
+      Comm comm = Comm::World();
+      Profile::Enable(true);
+
+      Long Nsurf = 1;
+      Stellarator<Real,ORDER> S;
+      Vector<Real> flux_tor(Nsurf), flux_pol(Nsurf);
+      { // Init S, flux_tor, flux_pol, pressure
+        Vector<Long> NtNp;
+        NtNp.PushBack(20);
+        NtNp.PushBack(4);
+        S = Stellarator<Real,ORDER>(NtNp);
+        flux_tor = 1;
+        flux_pol = 1;
+      }
+      Setup(S, comm);
+
+      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> nu(Nnodes);
+      { // Set nu
+        nu = area_elem;
+        //nu = 1;
+      }
+
+      Vector<ElemBasis> B, n_dot_dBdn, dBdnu;
+      { // Set B, 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);
+        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_;
+          }
+        }
+      }
+      { // Set dBdnu
+        Real eps = 1e-4;
+        Stellarator<Real,ORDER> S0 = S, S1 = S;
+        for (Long i = 0; i < Nelem; i++) {
+          for (Long j = 0; j < Nnodes; j++) {
+            S0.Elem(i, 0)[j] -= 0.5 * eps * normal[i*COORD_DIM+0][j] * nu[i][j];
+            S0.Elem(i, 1)[j] -= 0.5 * eps * normal[i*COORD_DIM+1][j] * nu[i][j];
+            S0.Elem(i, 2)[j] -= 0.5 * eps * normal[i*COORD_DIM+2][j] * nu[i][j];
+
+            S1.Elem(i, 0)[j] += 0.5 * eps * normal[i*COORD_DIM+0][j] * nu[i][j];
+            S1.Elem(i, 1)[j] += 0.5 * eps * normal[i*COORD_DIM+1][j] * nu[i][j];
+            S1.Elem(i, 2)[j] += 0.5 * eps * normal[i*COORD_DIM+2][j] * nu[i][j];
+          }
+        }
+        Setup(S0, comm);
+        Setup(S1, comm);
+
+        Real alpha0, alpha1, beta0, beta1;
+        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);
+      }
+
+      { // Write VTU
+        VTUData vtu;
+        vtu.AddElems(S.GetElemList(), B, ORDER);
+        vtu.WriteVTK("B", comm);
+      }
+      { // Write VTU
+        VTUData vtu;
+        vtu.AddElems(S.GetElemList(), n_dot_dBdn, ORDER);
+        vtu.WriteVTK("n_dot_dBdn", comm);
+      }
+      { // Write VTU
+        VTUData vtu;
+        vtu.AddElems(S.GetElemList(), dBdnu, ORDER);
+        vtu.WriteVTK("dBdnu", comm);
+      }
+
+#if 0
       Comm comm = Comm::World();
       Profile::Enable(true);
 
@@ -4646,6 +4756,7 @@ template <class Real, Integer ORDER=10> class Stellarator {
         vtu.AddElems(S.GetElemList(), dg_dnu1, ORDER);
         vtu.WriteVTK("dg_dnu1", comm);
       }
+#endif
     }
 
   private:

src/test.cpp

@@ -60,10 +60,10 @@ void TestMatrix() {
 int main(int argc, char** argv) {
   sctl::Comm::MPI_Init(&argc, &argv);
 
-  sctl::MHDEquilib<double>::test();
-  return 0;
+  //sctl::MHDEquilib<double>::test();
+  //return 0;
 
-  sctl::Stellarator<double>::test();
+  sctl::Stellarator<double>::test_askham();
   return 0;
 
   sctl::Spheres<double, 25>::test();