Erased the prints. 3d version compiling and running.

3 files changed, 50 insertions, 20 deletions

diff --git a/Source/BoundaryConditions/PML_routines.F90 b/Source/BoundaryConditions/PML_routines.F90
index e6e7f25b7..23cd05e89 100644
--- a/Source/BoundaryConditions/PML_routines.F90
+++ b/Source/BoundaryConditions/PML_routines.F90
@@ -297,8 +297,8 @@ contains
     integer :: i, j, k
     real(amrex_real) :: alpha_xy, alpha_xz, alpha_yx, alpha_yz, alpha_zx, alpha_zy
 
-    PRINT *, "hello world!"
-    PRINT *, "pml_type = ", pml_type(1)
+    ! PRINT *, ">>> PML_routine"
+    ! PRINT *, "pml_type = ", pml_type(1)
 
     if (flag == 0) then
       do       k = xlo(3), xhi(3)
@@ -485,8 +485,17 @@ contains
            do    j = xlo(2), xhi(2)
               do i = xlo(1), xhi(1)
                  ! compute current coefficients alpha
-                 alpha_xy = sigjy(j)/(sigjy(j)+sigjz(k))
-                 alpha_xz = sigjz(k)/(sigjy(j)+sigjz(k))
+                 if ((sigjy(j)==0.) .AND. (sigjz(k)==0.)) then
+                    alpha_xy = 0.5
+                    alpha_xz = 0.5
+                 else
+                    alpha_xy = sigjy(j)/(sigjy(j)+sigjz(k))
+                    alpha_xz = sigjz(k)/(sigjy(j)+sigjz(k))
+                 end if
+                !  PRINT *, "alpha_xy = ", alpha_xy
+                !  PRINT *, "alpha_xz = ", alpha_xz
+                !  PRINT *, "sigsjy = ", sigsjy(j)
+                !  PRINT *, "sigsjz = ", sigsjz(k)
                  Ex(i,j,k,1) = Ex(i,j,k,1) + dtsdy*(Bz(i,j  ,k  ,1)+Bz(i,j  ,k  ,2) &
                       &                            -Bz(i,j-1,k  ,1)-Bz(i,j-1,k  ,2))&
                       &                            -mudt*alpha_xy*jx(i,j,k)
@@ -500,8 +509,18 @@ contains
         do       k = ylo(3), yhi(3)
            do    j = ylo(2), yhi(2)
               do i = ylo(1), yhi(1)
-                 alpha_yx = sigjx(i)/(sigjx(i)+sigjz(k))
-                 alpha_yz = sigjz(k)/(sigjx(i)+sigjz(k))
+                 if ((sigjx(i)==0.) .AND. (sigjz(k)==0.)) then
+                    alpha_yx = 0.5
+                    alpha_yz = 0.5
+                 else
+                    alpha_yx = sigjx(i)/(sigjx(i)+sigjz(k))
+                    alpha_yz = sigjz(k)/(sigjx(i)+sigjz(k))
+                 end if
+
+                !  PRINT *, "alpha_yx = ", alpha_yx
+                !  PRINT *, "alpha_yz = ", alpha_yz
+                !  PRINT *, "sigsjx = ", sigsjx(i)
+                !  PRINT *, "sigsjz = ", sigsjz(k)
                  Ey(i,j,k,1) = Ey(i,j,k,1) + dtsdz*(Bx(i  ,j,k  ,1)+Bx(i  ,j,k  ,2) &
                       &                            -Bx(i  ,j,k-1,1)-Bx(i  ,j,k-1,2))&
                       &                            -mudt*alpha_yx*jy(i,j,k)
@@ -515,8 +534,18 @@ contains
         do       k = zlo(3), zhi(3)
            do    j = zlo(2), zhi(2)
               do i = zlo(1), zhi(1)
-                 alpha_zx = sigjx(i)/(sigjx(i)+sigjy(j))
-                 alpha_zy = sigjy(j)/(sigjx(i)+sigjy(j))
+                 if ((sigjx(i)==0.) .AND. (sigjy(j)==0.)) then
+                    alpha_zx = 0.5
+                    alpha_zy = 0.5
+                 else
+                    alpha_zx = sigjx(i)/(sigjx(i)+sigjy(j))
+                    alpha_zy = sigjy(j)/(sigjx(i)+sigjy(j))
+                 end if
+
+                !  PRINT *, "alpha_zx = ", alpha_zx
+                !  PRINT *, "alpha_zy = ", alpha_zy
+                !  PRINT *, "sigsjx = ", sigsjx(i)
+                !  PRINT *, "sigsjy = ", sigsjy(j)
                  Ez(i,j,k,1) = Ez(i,j,k,1) + dtsdx*(By(i  ,j  ,k,1)+By(i  ,j  ,k,2) &
                       &                            -By(i-1,j  ,k,1)-By(i-1,j  ,k,2))&
                       &                            -mudt*alpha_zx*jz(i,j,k)
diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp
index 470053652..ab5b28ade 100644
--- a/Source/Evolve/WarpXEvolveEM.cpp
+++ b/Source/Evolve/WarpXEvolveEM.cpp
@@ -298,11 +298,11 @@ WarpX::OneStep_nosub (Real cur_time)
     FillBoundaryE();
     FillBoundaryB();
 #else
-    amrex::Print()<< "WarpXEvolveEM.cpp : before CopyJinPMLs "<<std::endl;
+    // amrex::Print()<< "WarpXEvolveEM.cpp : before CopyJinPMLs "<<std::endl;
     if (do_pml && pml_has_particles){
         // do current deposition in PMLs
         // copy current computed on mother grid to PMLs
-        amrex::Print()<< "WarpXEvolveEM.cpp : IN CopyJinPMLs "<<std::endl;
+        // amrex::Print()<< "WarpXEvolveEM.cpp : IN CopyJinPMLs "<<std::endl;
         for (int lev = 0; lev <= finest_level; ++lev)
         {
             if (pml[lev]->ok()){
@@ -317,7 +317,7 @@ WarpX::OneStep_nosub (Real cur_time)
     }
 
     if (do_pml && do_pml_j_damping){
-        amrex::Print()<< "WarpXEvolveEM.cpp : DampJ "<<std::endl;
+        // amrex::Print()<< "WarpXEvolveEM.cpp : DampJ "<<std::endl;
         // damp current in pmls
         // amrex::Print() << "===== DAMPING IN PMLs =====" << std::endl;
 
@@ -336,11 +336,11 @@ WarpX::OneStep_nosub (Real cur_time)
     FillBoundaryE();
     EvolveF(0.5*dt[0], DtType::SecondHalf);
 
-    amrex::Print()<< "WarpXEvolveEM.cpp : before CopyJinReg "<<std::endl;
+    // amrex::Print()<< "WarpXEvolveEM.cpp : before CopyJinReg "<<std::endl;
     if (do_pml && pml_has_particles){
         // do current deposition in PMLs
         // copy current computed on mother grid to PMLs
-        amrex::Print()<< "WarpXEvolveEM.cpp : IN CopyJinReg "<<std::endl;
+        // amrex::Print()<< "WarpXEvolveEM.cpp : IN CopyJinReg "<<std::endl;
         for (int lev = 0; lev <= finest_level; ++lev)
         {
             if (pml[lev]->ok()){
@@ -361,7 +361,7 @@ WarpX::OneStep_nosub (Real cur_time)
 
     EvolveB(0.5*dt[0]); // We now have B^{n+1}
     if (do_pml) {
-        amrex::Print()<< "WarpXEvolveEM.cpp : DampPML "<<std::endl;
+        // amrex::Print()<< "WarpXEvolveEM.cpp : Damp "<<std::endl;
         DampPML();
         FillBoundaryE();
     }
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index abffb7ab6..298b03dc6 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -332,14 +332,15 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt)
             const Box& tey  = mfi.tilebox(Ey_nodal_flag);
             const Box& tez  = mfi.tilebox(Ez_nodal_flag);
 
-            amrex::Print()<< "== sigba.pml_type_array"<<std::endl;
+            // amrex::Print()<< "== sigba.pml_type_array"<<std::endl;
             const auto& pml_type = sigba[mfi].pml_type_array;
-            amrex::Print()<<"&pml_type = "<< &pml_type <<std::endl;
-            amrex::Print()<<"tex.loVect() = "<< tex.loVect() <<std::endl;
-            amrex::Print()<<"sigba[mfi].sigma[0].m_lo = "<< sigba[mfi].sigma[0].m_lo <<std::endl;
+            // amrex::Print()<<"&pml_type = "<< &pml_type <<std::endl;
+            // amrex::Print()<<"pml_type[0] = "<< pml_type[0] <<std::endl;
+            // amrex::Print()<<"tex.loVect() = "<< tex.loVect() <<std::endl;
+            // amrex::Print()<<"sigba[mfi].sigma[0].m_lo = "<< sigba[mfi].sigma[0].m_lo <<std::endl;
             // amrex::Print()<<"sigba[mfi].sigma[0].m_loVec() = "<< sigba[mfi].sigma[0].lo <<std::endl;
 
-            amrex::Print()<<"===== WRPX_PUSH_PML_EVEC ====="<<std::endl;
+            // amrex::Print()<<"===== WRPX_PUSH_PML_EVEC ====="<<std::endl;
             WRPX_PUSH_PML_EVEC(
 			     tex.loVect(), tex.hiVect(),
 			     tey.loVect(), tey.hiVect(),
@@ -367,7 +368,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt)
                              &mu_c2_dt,
 #endif
                              &dtsdx_c2, &dtsdy_c2, &dtsdz_c2);
-            amrex::Print()<<"===== WRPX_PUSH_PML_EVEC FIN ====="<<std::endl;
+            // amrex::Print()<<"===== WRPX_PUSH_PML_EVEC FIN ====="<<std::endl;
             if (pml_F)
             {
                 WRPX_PUSH_PML_EVEC_F(