Tweaks for looking at padding cells

examples/EMSchur3D.cpp

 int main( int argc, char** argv ) {
 
     // BiCGSTAB Diagonal preconditioner
-    //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::IncompleteLUT<Complex> > >::NewSP();
+    auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::IncompleteLUT<Complex> > >::NewSP();
 
     //auto EM3D = EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::NewSP();
 
     //auto EM3D = EMSchur3D< Eigen::SparseLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::NewSP();
 
     // UmfPack
-    auto EM3D = EMSchur3D< Eigen::UmfPackLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::NewSP();
+    //auto EM3D = EMSchur3D< Eigen::UmfPackLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR> > >::NewSP();
 
     // PaStiX
     //auto EM3D = EMSchur3D< Eigen::PastixLU<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, true > >::NewSP();

examples/block.py

 
 def calcH( oxx, nxx, dxx, pad ):
     """Calculates grid spacing along a dimension, if padding cells are requested those are incorporated too"""
-    pp = 1.1 # fixed padding expansion factor 
+    px = np.arange(oxx, oxx+(nxx+1)*dxx, dxx)
+    hx = dxx*np.ones(nxx)
     if not pad:
-        px = np.arange(oxx-dxx/2, oxx-dxx/2+(nxx+1)*dxx, dxx)
-        hx = dxx*np.ones(nxx)
+        #px = np.arange(oxx-dxx/2, oxx-dxx/2+(nxx+1)*dxx, dxx)
         return hx, px 
     else:
-        hx = np.zeros( nxx )
-        for ip in range(0,pad):
-            hx[ip] = dxx*pp**(pad-ip)      
-        hx[pad:nxx-pad] = dxx
-        ip = 1 
-        for ii in range (nxx-pad, nxx):
-            hx[ii] = dxx*pp**ip
-            ip += 1  
-
-        px = np.zeros( nxx + 1 )  # interfaces 
-        px[0] = oxx - (dxx*pp**pad)/2
-        for i, h in enumerate(hx):
-            px[i+1] = px[i]+h
+        pp = 1.1 # fixed padding expansion factor 
+        px_padhi = np.zeros( pad+1 )
+        px_padlo = np.zeros( pad+1 )
+        hx_pad = np.zeros( pad )
+        hx_pad[0] = dxx*pp
+        px_padhi[0] = px[-1] # hi end 
+        px_padlo[0] = px[0]  # lo end 
+        for ip in range(1,pad):
+            px_padhi[ip] = px_padhi[ip-1] + hx_pad[ip-1]    
+            px_padlo[ip] = px_padlo[ip-1] - hx_pad[ip-1]    
+            hx_pad[ip] = pp*hx_pad[ip-1]   
+        px_padhi[-1] = px_padhi[-2] + hx_pad[-1]    
+        px_padlo[-1] = px_padlo[-2] - hx_pad[-1]    
+
+        hx = np.concatenate((hx_pad[::-1],hx, hx_pad))
+        px = np.concatenate((px_padlo[1::][::-1], px, px_padhi[1::]))
 
         return hx, px
 
             fout.write(" ")
 
 def writeGRD(fout, args):
-    
-    fout.write( "%i\t%i\t%i\t// nx ny nz\n" %(args.nx[0],args.ny[0],args.nz[0]) )
-    fout.write( "%f\t%f\t%f\t// ox oy oz\n" %(args.ox[0],args.oy[0],args.oz[0]) )
 
-    fout.write("\n// hx\n")
     hx,px = calcH( args.ox[0], args.nx[0], args.dx[0], args.pad )
-    writeH(fout, hx, px)
+    hy,py = calcH( args.oy[0], args.ny[0], args.dy[0], args.pad )
+    hz,pz = calcH( args.oz[0], args.nz[0], args.dz[0], args.pad )
+    
+    fout.write( "%i %i %i // nx ny nz\n" %(args.nx[0]+2*args.pad,args.ny[0]+2*args.pad,args.nz[0]+2*args.pad) )
+    #fout.write( "%f %f %f // ox oy oz (cell centres)\n" %(args.ox[0]+args.dx[0]/2, args.oy[0]+args.dy[0]/2, args.oz[0]+args.dz[0]/2))
+    fout.write( "%f %f %f // ox oy oz (cell centres)\n" %(px[0]+.5*hx[0], py[0]+.5*hy[0], pz[0]+.5*hz[0]))
     
+    fout.write("\n// hx\n")
+    writeH(fout, hx, px)
     fout.write("\n// hy\n")
-    hy,py = calcH( args.oy[0], args.ny[0], args.dy[0], args.pad )
     writeH(fout, hy, py)
-    
     fout.write("\n// hz\n")
-    hz,pz = calcH( args.oz[0], args.nz[0], args.dz[0], args.pad )
     writeH(fout, hz, pz)
 
-    fout.write("\n\n")
-    #fout.write( "// vim: set tabstop=4 expandtab:\n")
+    fout.write("\n")
     fout.write( "// vim: set filetype=cpp:\n")
 
     return px,py,pz
 
-def writeMOD(fout, args, px, py, pz):
-    sigma = np.zeros( (args.nx[0], args.ny[0], args.nz[0]) )  
+def writeMOD(fout, args, px, py, pz, pad=0):
+    sigma = np.zeros( (args.nx[0]+2*pad, args.ny[0]+2*pad, args.nz[0]+2*pad) )  
     sigma[:,:,pz[1:]>0] = 0.05
     xlo = np.where( px[0:-1] >= args.bxlo[0] )[0][0] 
     xhi = np.where( px[1:] <= args.bxhi[0] )[0][-1] + 1
 parser = argparse.ArgumentParser(description='Set up test EMSchur3D block problem.')
 
 parser.add_argument('nx', metavar='nx', type=int, nargs=1, 
-                   help='number of cells in x')
+                   help='number of cells in x, does NOT include padding cells')
 
 parser.add_argument('ny', metavar='ny', type=int, nargs=1, 
-                   help='number of cells in y')
+                   help='number of cells in y, does NOT include padding cells')
 
 parser.add_argument('nz', metavar='nz', type=int, nargs=1, 
-                   help='number of cells in z')
+                   help='number of cells in z, does NOT include padding cells')
 
 parser.add_argument('ox', metavar='ox', type=float, nargs=1, 
-                   help='offset of grid in x')
+                   help='offset of grid in x (low edge of cell not centre), if pad is true ox should NOT include padding')
 
 parser.add_argument('oy', metavar='oy', type=float, nargs=1, 
-                   help='offset of grid in y')
+                   help='offset of grid in y (low edge of cell not centre), if pad is true oy should NOT include padding')
 
 parser.add_argument('oz', metavar='oz', type=float, nargs=1, 
-                   help='offset of grid in z')
+                   help='offset of grid in z (low edge of cell not centre), if pad is true oz should NOT include padding')
 
 parser.add_argument('dx', metavar='dx', type=int, nargs=1, 
                    help='minimum spacing in x')
 fout.close()
 
 mout = open("example.mod", 'w')
-writeMOD(mout, args, px, py, pz)
+writeMOD(mout, args, px, py, pz, args.pad)
 
 #print(args.nx)
 #print(args.accumulate(args.nx))

include/EMSchur3D.h

         success = implicitbicgstab(D, idx, ms, ADiv, Phi, CSolver[iw], max_it, tol, errorn, iter_done, logio);
         //Phi.array() *= MAC.array().cast<Complex>(); // remove phi from air regions
 
+
         /* Restart if necessary */
+/*
         int nrestart(1);
         // TODO send MAC to implicitbicgstab?
         while (success == 2 && nrestart < 18 && iter_done > 1) {
             //Phi.array() *= MAC.array().cast<Complex>(); // remove phi from air regions
             nrestart += 1;
         }
+*/
 
         logio << "Implicit BiCGStab solution in " << iter_done << " iterations."
                 << " with error " << std::setprecision(8) << std::scientific << errorn << std::endl;
     void EMSchur3D< Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::IncompleteLUT<Complex> > > ::BuildCDirectSolver() {
         CSolver = new Eigen::BiCGSTAB<Eigen::SparseMatrix<Complex, SPARSEMAJOR>, Eigen::IncompleteLUT<Complex> > [Omegas.size()];
         for (int iw=0; iw<Omegas.size(); ++iw) {
-            Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
-            CSolver[iw].preconditioner().setDroptol(1e-6);       //1e-5);      // 1e-12
-            CSolver[iw].preconditioner().setFillfactor(5e1);     // 1e2
+            //Csym = Cvec[iw].selfadjointView<Eigen::Lower>();
+            CSolver[iw].preconditioner().setDroptol(1e-5);       //1e-5);      // 1e-12
+            CSolver[iw].preconditioner().setFillfactor(1e1);     // 1e2
             jsw_timer timer;
             timer.begin();
             /*  Complex system */
             std::cout << "BiCGSTAB(ILU) pattern analyzing C_" << iw << ",";
             std::cout.flush();
-            CSolver[iw].analyzePattern( Csym );
-            //CSolver[iw].analyzePattern( Cvec[iw]);
+            //CSolver[iw].analyzePattern( Csym );
+            CSolver[iw].analyzePattern( Cvec[iw]);
             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
             /* factorize */
             timer.begin();
             std::cout << "BiCGSTAB(ILU) factorising C_" << iw << ", ";
             std::cout.flush();
-            CSolver[iw].factorize( Csym );
-            //CSolver[iw].factorize( Cvec[iw] );
+            //CSolver[iw].factorize( Csym );
+            CSolver[iw].factorize( Cvec[iw] );
             std::cout << " done in " << timer.end() / 60. << " [m]" << std::endl;
         }
     }

include/EMSchur3DBase.h

 
 namespace Lemma {
 
-// Pardiso likes LOWER for Sym problem
-#define UPPER 1  // LOWER WAS 0
+// Pardiso likes LOWER for Sym problem,
+//         for non-symmetric, set both to 1
+#define UPPER 1  // 1=true, 0=false
 #define LOWER 1  // 1=true, 0=false
 
 enum SOLVER{ SPARSELU, SimplicialLLT, SimplicialLDLT, BiCGStab, SparseQR };
     std::vector<int> idx;
 
     /** Dirichlet on low X */
-    //bool DirichletXLO = true;
-    bool DirichletXLO = false;
+    bool DirichletXLO = true;
+    //bool DirichletXLO = false;
 
-    //bool DirichletXHI = true;
-    bool DirichletXHI = false;
+    bool DirichletXHI = true;
+    //bool DirichletXHI = false;
 
-    //bool DirichletYLO = true;
-    bool DirichletYLO = false;
+    bool DirichletYLO = true;
+    //bool DirichletYLO = false;
 
-    //bool DirichletYHI = true;
-    bool DirichletYHI = false;
+    bool DirichletYHI = true;
+    //bool DirichletYHI = false;
 
-//    bool DirichletZLO = true;
-    bool DirichletZLO = false;
+    bool DirichletZLO = true;
+    //bool DirichletZLO = false;
 
     bool DirichletZHI = true;
-//    bool DirichletZHI = false;
+    //bool DirichletZHI = false;
 
     private: