.

include/sctl/boundary_quadrature.hpp

@@ -1133,7 +1133,7 @@ template <class Real> class Quadrature {
       constexpr Integer KDIM1 = Kernel::TrgDim();
       const Long Nelem = elem_lst.NElem();
 
-      BuildNbrList(pair_lst, Xt_, trg_surf, elem_lst, 2.5/order_direct, period_length, comm);
+      BuildNbrList(pair_lst, Xt_, trg_surf, elem_lst, 5.0/order_direct, period_length, comm);
       const Long Ninterac = pair_lst.Dim();
 
       Vector<Real> Xt;
@@ -1753,7 +1753,7 @@ template <class Real> class Quadrature {
       Profile::Toc();
     }
 
-    template <Integer ORDER = 5> static void test(Integer order_singular = 10, Integer order_direct = 5, const Comm& comm = Comm::World()) {
+    template <Integer ORDER = 10> static void test(Integer order_singular = 20, Integer order_direct = 30, const Comm& comm = Comm::World()) {
       constexpr Integer COORD_DIM = 3;
       constexpr Integer ELEM_DIM = COORD_DIM-1;
       using ElemList = ElemList<COORD_DIM, Basis<Real, ELEM_DIM, ORDER>>;
@@ -1763,14 +1763,76 @@ template <class Real> class Quadrature {
       int np = comm.Size();
       int rank = comm.Rank();
       auto build_torus = [rank,np](ElemList& elements, long Nt, long Np, Real Rmajor, Real Rminor){
-        auto nodes = ElemList::CoordBasis::Nodes();
         auto torus = [](Real theta, Real phi, Real Rmajor, Real Rminor) {
-          Real R = Rmajor + Rminor * cos<Real>(phi);
-          Real X = R * cos<Real>(theta);
-          Real Y = R * sin<Real>(theta);
-          Real Z = Rminor * sin<Real>(phi);
+          Real s = Rminor;
+          Integer Nperiod = 5;
+    #if 0
+          Real Aspect_ratio = 10.27932548522949;
+          Real coeffmat[21][21] = { 0.00000478813217,  0.00000000000000,  0.00000351611652,  0.00000135354389,  0.00000061357832,  0.00000220091101,  0.00000423862912, -0.00003000058678,  0.00000064187111, -0.00024228452821,  0.00003116775770,  0.00000176210710,  0.00000289141326, -0.00000150300525,  0.00000772853855,  0.00000098855242,  0.00000316606793,  0.00000002168364,  0.00000212047939,  0.00000299016097,  0.00000443224508,
+                                    0.00000028202930,  0.00000000000000, -0.00000249222421, -0.00000203136278,  0.00000131104809,  0.00000011987446, -0.00000370760154,  0.00004553918916, -0.00007711342914, -0.00004685295062,  0.00011049838213, -0.00000197486270,  0.00000395827146,  0.00000615046474,  0.00000755337123,  0.00000700606006,  0.00000922725030, -0.00000043310337,  0.00000107416383,  0.00000449787694,  0.00000305137178,
+                                    0.00001226376662,  0.00000000000000,  0.00000270820692,  0.00000208059305,  0.00000521478523,  0.00001779037302,  0.00000846544117,  0.00001120913385, -0.00065816845745, -0.00085107452469, -0.00013171190221, -0.00005540943675, -0.00001835885450,  0.00000101879823,  0.00000209222071,  0.00000091532502, -0.00000521515358, -0.00000209227142, -0.00000678545939, -0.00000034963549, -0.00000015111488,
+                                    0.00001560274177,  0.00000000000000,  0.00000350691471, -0.00001160475040, -0.00001763036562,  0.00003487367940, -0.00002787247831, -0.00000910982726,  0.00008818832430, -0.00524408789352,  0.00009378376126,  0.00004184526188,  0.00002849263365, -0.00002757280527,  0.00003388467667,  0.00000706207265,  0.00000625263419, -0.00003315929280, -0.00001181772132,  0.00000311426015,  0.00001875682574,
+                                   -0.00000398287420,  0.00000000000000, -0.00001524541040,  0.00001724056165,  0.00002245173346,  0.00002806861812, -0.00000388776925,  0.00008143573359, -0.00005900909309,  0.00110496615525,  0.00134626252111,  0.00005128383054, -0.00001372421866,  0.00003612563887,  0.00002236580076, -0.00002728391883,  0.00001981237256,  0.00000655450458,  0.00000985319002,  0.00001347597299,  0.00000645987802,
+                                    0.00003304968050,  0.00000000000000, -0.00000530822217,  0.00001324870937, -0.00003610889689, -0.00005478735329, -0.00005818806312, -0.00037112057908, -0.00017812002625, -0.00093204283621,  0.00115969858598, -0.00033559172880, -0.00010441876657, -0.00001617923044, -0.00000555065844,  0.00007343527250, -0.00004408047607,  0.00000403802142,  0.00001843931204,  0.00001694047933,  0.00001213414362,
+                                   -0.00000751115658,  0.00000000000000,  0.00005457974839, -0.00000334614515,  0.00005845565465,  0.00015000770509,  0.00021849104087,  0.00002724147635,  0.00167233624961,  0.00011666602222,  0.00276563479565, -0.00085952825611, -0.00030217235326, -0.00008841593808,  0.00000997664119, -0.00015285826521,  0.00002517224675,  0.00003009161810,  0.00001883217556,  0.00002146127554,  0.00001822445302,
+                                   -0.00004128706860,  0.00000000000000, -0.00003496417776,  0.00001088761655, -0.00000298955979, -0.00005359326315, -0.00019021633489, -0.00017992728681, -0.00347794801928,  0.00064632791327,  0.00449698418379, -0.00017710507382,  0.00006126180233,  0.00018059254216,  0.00002354096432,  0.00008189838991, -0.00010060678323, -0.00017183290038,  0.00019413756672,  0.00021334811754,  0.00011263617489,
+                                    0.00000853522670, -0.00000000000000, -0.00006544789358,  0.00005424076880, -0.00000679056529, -0.00001249735487, -0.00053082982777,  0.00035396864405, -0.00115020677913,  0.05894451215863,  0.06573092192411,  0.01498018857092,  0.00278125284240,  0.00145188067108,  0.00033717858605,  0.00000800427370, -0.00009335305367,  0.00024286781263, -0.00023916347709,  0.00031213948387,  0.00018134393031,
+                                   -0.00002521496390, -0.00000000000000, -0.00054337945767,  0.00012690725271,  0.00053313979879,  0.00064233405283, -0.00047686311882,  0.00176536326762,  0.00074157933705, -0.02684566564858,  1.00000000000000,  0.07176169008017,  0.00837037432939, -0.00000381640211,  0.00088998704450, -0.00049218931235, -0.00024546548957, -0.00036608282244,  0.00049480766756,  0.00031158892671,  0.00006898906577,
+                                    0.00021280418150,  0.00028127161204, -0.00070030166535,  0.00022237010126, -0.00028713891516, -0.00013800295710,  0.00005912094275,  0.00172126013786, -0.00618684850633,  0.03608432412148,    Aspect_ratio  ,  0.49896776676178,  0.00091372377938, -0.00085712829605, -0.00124801427592, -0.00007427225501, -0.00005245858847,  0.00002841771493,  0.00020249813679, -0.00014303345233,  0.00001406490901,
+                                    0.00023699452868,  0.00008661757602,  0.00025744654704, -0.00022715188970, -0.00076146807987,  0.00055185536621, -0.00012325309217, -0.00072356045712, -0.00160693109501,  0.00246682553552, -0.14175094664097, -0.36207047104836, -0.04089594259858,  0.00060774467420,  0.00088646943914,  0.00004865296432, -0.00041878610500, -0.00023025234987, -0.00009676301852, -0.00000000000000,  0.00008409228758,
+                                    0.00011432896281, -0.00000707848403,  0.00004698805787, -0.00043642931269,  0.00081384339137, -0.00065635429928, -0.00011831733718,  0.00017413357273,  0.00224463525228,  0.00478497287259,  0.03294761106372,  0.01078986655921,  0.10731782764196,  0.00075034319889, -0.00009241879889,  0.00055023463210,  0.00006596000458,  0.00005045382932,  0.00014874986664,  0.00000000000000, -0.00015369028552,
+                                    0.00001037383754,  0.00009250180301,  0.00026204055757,  0.00007424291834, -0.00047751804232,  0.00029184055165,  0.00050921301590, -0.00004825839278, -0.00029933769838,  0.00279659987427,  0.00210463814437, -0.00618590926751, -0.02400829829276, -0.02316811867058, -0.00086368201301, -0.00032258985448, -0.00018304496189,  0.00008438774967, -0.00008305341908,  0.00000000000000,  0.00013047417451,
+                                   -0.00001376930322, -0.00001723831701, -0.00011543079017, -0.00022646733851,  0.00013467084500, -0.00004661652201, -0.00008419520600,  0.00035772417323, -0.00011815709877,  0.00028718306567,  0.00092207465786, -0.00317224999890,  0.00061770365573,  0.01017294172198,  0.00294739892706,  0.00014669894881,  0.00015702951350,  0.00003432080121, -0.00008555022214, -0.00000000000000,  0.00000454909878,
+                                   -0.00000196001542, -0.00003198397462, -0.00004425687075, -0.00004129848094, -0.00003789070615, -0.00027583551127,  0.00025874207495, -0.00002334945384, -0.00007259396807, -0.00008295358566,  0.00011360697681, -0.00101968157105,  0.00046784928418, -0.00208410434425, -0.00313158822246, -0.00046005158219, -0.00010552268213, -0.00005850767775,  0.00003971093611,  0.00000000000000, -0.00005275657168,
+                                   -0.00001065901233, -0.00001934838656, -0.00001220186732, -0.00002060524639, -0.00000225423423, -0.00001894621164, -0.00001533334580, -0.00001791087379,  0.00008156246622, -0.00008441298269,  0.00021060956351, -0.00030303673702,  0.00075949780876, -0.00010539998038,  0.00109045265708,  0.00068949378328,  0.00009268362192,  0.00003471063246,  0.00001204656473, -0.00000000000000,  0.00001500743110,
+                                    0.00000105878155, -0.00000910870767, -0.00000172467264, -0.00000722095228,  0.00000699280463, -0.00002061720625, -0.00000889817693, -0.00001993474507,  0.00000370749740, -0.00000090311920,  0.00002677819793,  0.00043428712524,  0.00210293265991,  0.00018200518389, -0.00009621794743, -0.00035250501242, -0.00012996385340, -0.00002185157609, -0.00001116586463, -0.00000000000000, -0.00000451994811,
+                                    0.00000424055270, -0.00000463139304,  0.00000301006116, -0.00000123974939,  0.00000632465435, -0.00002090823000,  0.00001773388794,  0.00000121050368,  0.00001886057362, -0.00001043497195, -0.00002269273500, -0.00021979617304, -0.00001043962493, -0.00116343051195, -0.00004193381756,  0.00007944958634,  0.00007301353617,  0.00002082651736, -0.00000119863023, -0.00000000000000, -0.00001440504820,
+                                   -0.00000391270805, -0.00000490489265, -0.00000504441778, -0.00000904507579, -0.00000111389932,  0.00000597532107,  0.00000047090245, -0.00001553130096, -0.00001524566323, -0.00000522222899, -0.00007707672921, -0.00004165665086,  0.00015764687851,  0.00035649110214,  0.00038701237645,  0.00002386798405, -0.00001946414341, -0.00000913835174, -0.00000489907188,  0.00000000000000,  0.00000172327657,
+                                   -0.00000015388650, -0.00000603232729, -0.00000397650865,  0.00000280493782,  0.00000463132073, -0.00000788678426, -0.00000471605335, -0.00000283715985, -0.00000422824724,  0.00000366817630, -0.00001159603562, -0.00001625759251,  0.00049116823357,  0.00005048640014, -0.00020234247495, -0.00006341376866, -0.00000807822744,  0.00000070463199,  0.00000014041755,  0.00000000000000, -0.00000718306910};
+    #else
+          Real Aspect_ratio = 5;
+          Real coeffmat[21][21] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            s, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, Aspect_ratio, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0.2*s, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0,
+                                   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,            0, 0    , 0, 0, 0, 0, 0, 0, 0, 0, 0};
+    #endif
+          Real Z = 0;
+          Real R = 0;
+          for (long i = -10; i <= 10; i++) {
+            for (long j = -10; j <= 10; j++) {
+              R += coeffmat[i+10][j+10] * cos(-i*phi + Nperiod*j*theta)*Rmajor;
+              Z += coeffmat[i+10][j+10] * sin(-i*phi + Nperiod*j*theta)*Rmajor;
+            }
+          }
+          Real X = R * cos(theta);
+          Real Y = R * sin(theta);
           return std::make_tuple(X,Y,Z);
         };
+        //auto torus = [](Real theta, Real phi, Real Rmajor, Real Rminor) {
+        //  Real R = Rmajor + Rminor * cos<Real>(phi);
+        //  Real X = R * cos<Real>(theta);
+        //  Real Y = R * sin<Real>(theta);
+        //  Real Z = Rminor * sin<Real>(phi);
+        //  return std::make_tuple(X,Y,Z);
+        //};
+        auto nodes = ElemList::CoordBasis::Nodes();
 
         long start = Nt*Np*(rank+0)/np;
         long end   = Nt*Np*(rank+1)/np;
@@ -1790,8 +1852,8 @@ template <class Real> class Quadrature {
         }
       };
       ElemList elements_src, elements_trg;
-      build_torus(elements_src, 28, 16, 2, 1.0);
-      build_torus(elements_trg, 29, 17, 2, 0.99);
+      build_torus(elements_src, 60, 12, 2, 1.0);
+      build_torus(elements_trg, 61, 13, 2, 0.99);
 
       Vector<Real> Xt;
       Vector<PotentialBasis> U_onsurf, U_offsurf;
@@ -5737,8 +5799,8 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       long Np = fft_r2c.Dim(1) / (Nt * 2);
       SCTL_ASSERT(fft_r2c.Dim(1) == Nt * Np * 2);
 
-      auto filter_fn = [](Real x2, Real sigma) {return exp(-x2/(2*sigma*sigma));};
-      //auto filter_fn = [](Real x2, Real sigma) {return (x2<sigma*sigma?1.0:0.0);};
+      //auto filter_fn = [](Real x2, Real sigma) {return exp(-x2/(2*sigma*sigma));};
+      auto filter_fn = [](Real x2, Real sigma) {return (x2<sigma*sigma?1.0:0.0);};
 
       Vector<Real> coeff(fft_r2c.Dim(1));
       for (long k = 0; k < dof; k++) {
@@ -6357,7 +6419,7 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       pressure_ = pressure;
       flux_tor_ = flux_tor;
       flux_pol_ = flux_pol;
-      iter = 57;
+      iter = 0;
     }
 
     Real operator()(const Eigen::VectorXd& x, Eigen::VectorXd& grad_direction, Eigen::VectorXd* grad_op_ptr = nullptr) {
@@ -6441,9 +6503,6 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
         grad_direction_orig_.ReInit(x.size());
         Vector<Real> dgdnu_fourier = cheb2grid(S_, dgdnu);
         for (Long i = 0; i < S_.Nsurf(); i++) { // Init grad_direction_, make it divergence-free, filter
-          const Long Mt = S_.NTor(i);
-          const Long Mp = S_.NPol(i);
-          const Long offset = S_.ElemDsp(i);
           const Long Nt = fourier_dim0;
           const Long Np = fourier_dim1;
 
@@ -6473,32 +6532,37 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
             Sop.GradInvSurfLap(GradInvLapF, dX, d2X, F, 1e-8, 100, 1.0);
 
             for (Long j = 0; j < Nt*Np; j++) { // grad_direction <-- normal * dgdnu - GradInvLapF
-              grad_direction[0*Nt*Np+j] = normal[0*Nt*Np+j] * dgdnu[j] - GradInvLapF[0*Nt*Np+j]*0;
-              grad_direction[1*Nt*Np+j] = normal[1*Nt*Np+j] * dgdnu[j] - GradInvLapF[1*Nt*Np+j]*0;
-              grad_direction[2*Nt*Np+j] = normal[2*Nt*Np+j] * dgdnu[j] - GradInvLapF[2*Nt*Np+j]*0;
+              grad_direction[0*Nt*Np+j] = normal[0*Nt*Np+j] * dgdnu[j] - GradInvLapF[0*Nt*Np+j];
+              grad_direction[1*Nt*Np+j] = normal[1*Nt*Np+j] * dgdnu[j] - GradInvLapF[1*Nt*Np+j];
+              grad_direction[2*Nt*Np+j] = normal[2*Nt*Np+j] * dgdnu[j] - GradInvLapF[2*Nt*Np+j];
             }
             { ////////////////////
               Vector<Real> grad_direction_orig(COORD_DIM*Nt*Np, (Iterator<Real>)grad_direction_orig_.begin() + COORD_DIM*i * (fourier_dim0*fourier_dim1), false);
-              grad_direction_orig = grad_direction;
+              //grad_direction_orig = grad_direction;
+              for (Long j = 0; j < Nt*Np; j++) { // grad_direction <-- normal * dgdnu - GradInvLapF
+                grad_direction_orig[0*Nt*Np+j] = normal[0*Nt*Np+j] * dgdnu[j] - GradInvLapF[0*Nt*Np+j]*0;
+                grad_direction_orig[1*Nt*Np+j] = normal[1*Nt*Np+j] * dgdnu[j] - GradInvLapF[1*Nt*Np+j]*0;
+                grad_direction_orig[2*Nt*Np+j] = normal[2*Nt*Np+j] * dgdnu[j] - GradInvLapF[2*Nt*Np+j]*0;
+              }
             }
-            fourier_filter(grad_direction, Nt, Np, 2);
+            fourier_filter(grad_direction, Nt, Np, S_.NPol(i)*1.2);
 
             for (Long k = 0; k < 50; k++) { // reparameterize
               Real max_resid = 0;
               Vector<Real> correction(COORD_DIM*Nt*Np);
-              for (Long i = 0; i < Nt*Np; i++) { // Set correction
-                Real resid = dgdnu[i];
-                resid -= normal[0*Nt*Np+i] * grad_direction[0*Nt*Np+i];
-                resid -= normal[1*Nt*Np+i] * grad_direction[1*Nt*Np+i];
-                resid -= normal[2*Nt*Np+i] * grad_direction[2*Nt*Np+i];
+              for (Long j = 0; j < Nt*Np; j++) { // Set correction
+                Real resid = dgdnu[j];
+                resid -= normal[0*Nt*Np+j] * grad_direction[0*Nt*Np+j];
+                resid -= normal[1*Nt*Np+j] * grad_direction[1*Nt*Np+j];
+                resid -= normal[2*Nt*Np+j] * grad_direction[2*Nt*Np+j];
                 max_resid = std::max(max_resid, fabs(resid));
 
-                correction[0*Nt*Np+i] = normal[0*Nt*Np+i] * resid;
-                correction[1*Nt*Np+i] = normal[1*Nt*Np+i] * resid;
-                correction[2*Nt*Np+i] = normal[2*Nt*Np+i] * resid;
+                correction[0*Nt*Np+j] = normal[0*Nt*Np+j] * resid;
+                correction[1*Nt*Np+j] = normal[1*Nt*Np+j] * resid;
+                correction[2*Nt*Np+j] = normal[2*Nt*Np+j] * resid;
               }
               std::cout<<max_resid<<' '; //////////////////////////////////
-              fourier_filter(correction, Nt, Np, 2);
+              fourier_filter(correction, Nt, Np, S_.NPol(i)*1.2);
 
               Real alpha = 0;
               { // Set alpha <-- (dgdnu - x.n, c.n) / (c.n, c.n)
@@ -6597,8 +6661,8 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
           X[i*COORD_DIM+2] = mhd_equilib.S_.Elem(i,2);
         }
         Vector<Real> X_fourier = cheb2grid(mhd_equilib.S_, X);
+        //X_fourier.Read(("X"+std::to_string(mhd_equilib.iter)+"_").c_str()); // Read from file
         x.resize(X_fourier.Dim());
-        X_fourier.Read(("X"+std::to_string(mhd_equilib.iter)+"_").c_str()); // Read from file
         for (Long i = 0; i < X_fourier.Dim(); i++) {
           x(i) = X_fourier[i];
         }
@@ -6631,8 +6695,8 @@ template <class Real, Integer ORDER=10> class MHDEquilib {
       { // Init S, flux_tor, flux_pol, pressure
         Vector<Long> NtNp;
         for (Long i = 0; i < Nsurf; i++) {
-          NtNp.PushBack(70);
-          NtNp.PushBack(14);
+          NtNp.PushBack(35);
+          NtNp.PushBack(7);
           //NtNp.PushBack(30);
           //NtNp.PushBack(4);
         }