From 92c4e93cdb0fb34423218d20b47fd265a7a614b8 Mon Sep 17 00:00:00 2001 From: ablelly Date: Mon, 10 Jun 2019 20:23:10 +0200 Subject: First modifications for current in PMLs (modification of Maxwell's equations in PMLs) --- Source/FieldSolver/WarpXPushFieldsEM.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index c53e13f8f..c84a23e51 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -137,6 +137,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); + const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp(); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -338,6 +339,12 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), +#if (AMREX_SPACEDIM==2) + BL_TO_FORTRAN_3D((*pml_j[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_j[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_j[2])[mfi]), + &mu_c2_dt, +#endif &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); if (pml_F) @@ -435,4 +442,3 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type) } } } - -- cgit v1.2.3 From 24ae8389d7b526a2c2042d90f539fb145e05ea7e Mon Sep 17 00:00:00 2001 From: ablelly Date: Mon, 10 Jun 2019 20:38:44 +0200 Subject: Add current in PML. Code is compiling. --- Source/BoundaryConditions/PML_routines.F90 | 16 ++++++++-------- Source/FieldSolver/WarpXPushFieldsEM.cpp | 3 ++- 2 files changed, 10 insertions(+), 9 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML_routines.F90 b/Source/BoundaryConditions/PML_routines.F90 index 9d28c491a..91911febc 100644 --- a/Source/BoundaryConditions/PML_routines.F90 +++ b/Source/BoundaryConditions/PML_routines.F90 @@ -440,16 +440,16 @@ contains Exlo(2), Exhi(2), Eylo(2), Eyhi(2), Ezlo(2), Ezhi(2), & Bxlo(2), Bxhi(2), Bylo(2), Byhi(2), Bzlo(2), Bzhi(2), & jxlo(2), jxhi(2), jylo(2), jyhi(2), jzlo(2), jzhi(2) - real(amrex_real), intent(in) :: dtsdx, dtsdy, dtsdz + real(amrex_real), intent(in) :: mudt, dtsdx, dtsdy, dtsdz real(amrex_real), intent(inout) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),2) real(amrex_real), intent(inout) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),2) real(amrex_real), intent(inout) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),2) real(amrex_real), intent(in ) :: Bx (Bxlo(1):Bxhi(1),Bxlo(2):Bxhi(2),2) real(amrex_real), intent(in ) :: By (Bylo(1):Byhi(1),Bylo(2):Byhi(2),2) real(amrex_real), intent(in ) :: Bz (Bzlo(1):Bzhi(1),Bzlo(2):Bzhi(2),2) - real(amrex_real), intent(in ) :: jx (jxlo(1):jxhi(1),jxlo(2):jxhi(2),1) - real(amrex_real), intent(in ) :: jy (jylo(1):jyhi(1),jylo(2):jyhi(2),1) - real(amrex_real), intent(in ) :: jz (jzlo(1):jzhi(1),jzlo(2):jzhi(2),1) + real(amrex_real), intent(in ) :: jx (jxlo(1):jxhi(1),jxlo(2):jxhi(2)) !jx (jxlo(1):jxhi(1),jxlo(2):jxhi(2),1) + real(amrex_real), intent(in ) :: jy (jylo(1):jyhi(1),jylo(2):jyhi(2)) !jy (jylo(1):jyhi(1),jylo(2):jyhi(2),1) + real(amrex_real), intent(in ) :: jz (jzlo(1):jzhi(1),jzlo(2):jzhi(2)) !jz (jzlo(1):jzhi(1),jzlo(2):jzhi(2),1) integer :: i, k @@ -457,7 +457,7 @@ contains do i = xlo(1), xhi(1) Ex(i,k,2) = Ex(i,k,2) - dtsdz*(By(i,k ,1)+By(i,k ,2) & & -By(i,k-1,1)-By(i,k-1,2))& - & - mudt * jx(j,k) + & - mudt * jx(i,k) end do end do @@ -465,10 +465,10 @@ contains do i = ylo(1), yhi(1) Ey(i,k,1) = Ey(i,k,1) + dtsdz*(Bx(i ,k ,1)+Bx(i ,k ,2) & & -Bx(i ,k-1,1)-Bx(i ,k-1,2))& - & - mudt * 0.5 * jy(j,k) + & - mudt * 0.5 * jy(i,k) Ey(i,k,2) = Ey(i,k,2) - dtsdx*(Bz(i ,k ,1)+Bz(i ,k ,2) & & -Bz(i-1,k ,1)-Bz(i-1,k ,2))& - & - mudt * 0.5 * jy(j,k) + & - mudt * 0.5 * jy(i,k) end do end do @@ -476,7 +476,7 @@ contains do i = zlo(1), zhi(1) Ez(i,k,1) = Ez(i,k,1) + dtsdx*(By(i ,k,1)+By(i ,k,2) & & -By(i-1,k,1)-By(i-1,k,2))& - & - mudt * jz(j,k) + & - mudt * jz(i,k) end do end do diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index c84a23e51..b468049d4 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -137,7 +137,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp(); + // const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp(); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -319,6 +319,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); + const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp(); const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) -- cgit v1.2.3 From 2fbcf7f95a241922f459871f5b779ba8db82590c Mon Sep 17 00:00:00 2001 From: ablelly Date: Tue, 11 Jun 2019 19:12:55 +0200 Subject: Add changes to be able to use the exectutable without particles in PMLs. Compiling, but produces a Segfault only when `pml_has_particles = 0`. --- Source/BoundaryConditions/PML_routines.F90 | 45 +++++++++++++++++++++--------- Source/FieldSolver/WarpXPushFieldsEM.cpp | 1 + Source/FortranInterface/WarpX_f.H | 1 + 3 files changed, 34 insertions(+), 13 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML_routines.F90 b/Source/BoundaryConditions/PML_routines.F90 index 91911febc..b81b9b0e5 100644 --- a/Source/BoundaryConditions/PML_routines.F90 +++ b/Source/BoundaryConditions/PML_routines.F90 @@ -434,12 +434,14 @@ contains & jx, jxlo, jxhi, & & jy, jylo, jyhi, & & jz, jzlo, jzhi, & + & pml_has_particles, & & mudt, dtsdx, dtsdy, dtsdz) & bind(c,name='warpx_push_pml_evec_2d') integer, intent(in) :: xlo(2), xhi(2), ylo(2), yhi(2), zlo(2), zhi(2), & Exlo(2), Exhi(2), Eylo(2), Eyhi(2), Ezlo(2), Ezhi(2), & Bxlo(2), Bxhi(2), Bylo(2), Byhi(2), Bzlo(2), Bzhi(2), & - jxlo(2), jxhi(2), jylo(2), jyhi(2), jzlo(2), jzhi(2) + jxlo(2), jxhi(2), jylo(2), jyhi(2), jzlo(2), jzhi(2), & + pml_has_particles real(amrex_real), intent(in) :: mudt, dtsdx, dtsdy, dtsdz real(amrex_real), intent(inout) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),2) real(amrex_real), intent(inout) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),2) @@ -455,28 +457,45 @@ contains do k = xlo(2), xhi(2) do i = xlo(1), xhi(1) - Ex(i,k,2) = Ex(i,k,2) - dtsdz*(By(i,k ,1)+By(i,k ,2) & - & -By(i,k-1,1)-By(i,k-1,2))& - & - mudt * jx(i,k) + if (pml_has_particles==1) then + Ex(i,k,2) = Ex(i,k,2) - dtsdz*(By(i,k ,1)+By(i,k ,2) & + & -By(i,k-1,1)-By(i,k-1,2))& + & - mudt * jx(i,k) + else + Ex(i,k,2) = Ex(i,k,2) - dtsdz*(By(i,k ,1)+By(i,k ,2) & + & -By(i,k-1,1)-By(i,k-1,2)) + end if end do end do do k = ylo(2), yhi(2) do i = ylo(1), yhi(1) - Ey(i,k,1) = Ey(i,k,1) + dtsdz*(Bx(i ,k ,1)+Bx(i ,k ,2) & - & -Bx(i ,k-1,1)-Bx(i ,k-1,2))& - & - mudt * 0.5 * jy(i,k) - Ey(i,k,2) = Ey(i,k,2) - dtsdx*(Bz(i ,k ,1)+Bz(i ,k ,2) & - & -Bz(i-1,k ,1)-Bz(i-1,k ,2))& - & - mudt * 0.5 * jy(i,k) + if (pml_has_particles==1) then + Ey(i,k,1) = Ey(i,k,1) + dtsdz*(Bx(i ,k ,1)+Bx(i ,k ,2) & + & -Bx(i ,k-1,1)-Bx(i ,k-1,2))& + & - mudt * 0.5 * jy(i,k) + Ey(i,k,2) = Ey(i,k,2) - dtsdx*(Bz(i ,k ,1)+Bz(i ,k ,2) & + & -Bz(i-1,k ,1)-Bz(i-1,k ,2))& + & - mudt * 0.5 * jy(i,k) + else + Ey(i,k,1) = Ey(i,k,1) + dtsdz*(Bx(i ,k ,1)+Bx(i ,k ,2) & + & -Bx(i ,k-1,1)-Bx(i ,k-1,2)) + Ey(i,k,2) = Ey(i,k,2) - dtsdx*(Bz(i ,k ,1)+Bz(i ,k ,2) & + & -Bz(i-1,k ,1)-Bz(i-1,k ,2)) + end if end do end do do k = zlo(2), zhi(2) do i = zlo(1), zhi(1) - Ez(i,k,1) = Ez(i,k,1) + dtsdx*(By(i ,k,1)+By(i ,k,2) & - & -By(i-1,k,1)-By(i-1,k,2))& - & - mudt * jz(i,k) + if (pml_has_particles==1) then + Ez(i,k,1) = Ez(i,k,1) + dtsdx*(By(i ,k,1)+By(i ,k,2) & + & -By(i-1,k,1)-By(i-1,k,2))& + & - mudt * jz(i,k) + else + Ez(i,k,1) = Ez(i,k,1) + dtsdx*(By(i ,k,1)+By(i ,k,2) & + & -By(i-1,k,1)-By(i-1,k,2)) + end if end do end do diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index b468049d4..4b3be0a43 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -344,6 +344,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) BL_TO_FORTRAN_3D((*pml_j[0])[mfi]), BL_TO_FORTRAN_3D((*pml_j[1])[mfi]), BL_TO_FORTRAN_3D((*pml_j[2])[mfi]), + &pml_has_particles, &mu_c2_dt, #endif &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index ed3edf1a7..00879487b 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -381,6 +381,7 @@ extern "C" BL_FORT_FAB_ARG_3D(jx), BL_FORT_FAB_ARG_3D(jy), BL_FORT_FAB_ARG_3D(jz), + const int* flag, const amrex::Real* mudt, #endif const amrex::Real* dtsdx, -- cgit v1.2.3 From 967ef322d47c11bc2b34149212342bf2e4ced1de Mon Sep 17 00:00:00 2001 From: ablelly Date: Wed, 12 Jun 2019 00:37:16 +0200 Subject: First intent to take ghost cells + communications into account. Compiling, but does not give the expected result. --- Source/BoundaryConditions/PML.cpp | 76 +++++++++++++++++++----------- Source/BoundaryConditions/PML_routines.F90 | 11 ++--- Source/Evolve/WarpXEvolveEM.cpp | 57 +++++++++++----------- Source/FieldSolver/WarpXPushFieldsEM.cpp | 2 +- 4 files changed, 82 insertions(+), 64 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 25ff66e8e..33ba3c44d 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -384,17 +384,23 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, pml_B_fp[1]->setVal(0.0); pml_B_fp[2]->setVal(0.0); - - if (pml_has_particles){ - pml_j_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::jx_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jx_nodal_flag) - pml_j_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::jy_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jy_nodal_flag) - pml_j_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::jz_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jz_nodal_flag) - pml_j_fp[0]->setVal(0.0); - pml_j_fp[1]->setVal(0.0); - pml_j_fp[2]->setVal(0.0); - amrex::Print() << "PML HAS PARTICLES - fine"<< std::endl; - - } + pml_j_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::jx_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jx_nodal_flag) + pml_j_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::jy_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jy_nodal_flag) + pml_j_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::jz_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jz_nodal_flag) + pml_j_fp[0]->setVal(0.0); + pml_j_fp[1]->setVal(0.0); + pml_j_fp[2]->setVal(0.0); + + // if (pml_has_particles){ + // pml_j_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::jx_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jx_nodal_flag) + // pml_j_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::jy_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jy_nodal_flag) + // pml_j_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::jz_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jz_nodal_flag) + // pml_j_fp[0]->setVal(0.0); + // pml_j_fp[1]->setVal(0.0); + // pml_j_fp[2]->setVal(0.0); + // amrex::Print() << "PML HAS PARTICLES - fine"<< std::endl; + // + // } if (do_dive_cleaning) { @@ -438,17 +444,22 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, pml_F_cp->setVal(0.0); } - - if (pml_has_particles) - { - pml_j_cp[0].reset(new MultiFab(amrex::convert(cba,WarpX::jx_nodal_flag), cdm, 1, ngb)); - pml_j_cp[1].reset(new MultiFab(amrex::convert(cba,WarpX::jy_nodal_flag), cdm, 1, ngb)); - pml_j_cp[2].reset(new MultiFab(amrex::convert(cba,WarpX::jz_nodal_flag), cdm, 1, ngb)); - pml_j_cp[0]->setVal(0.0); - pml_j_cp[1]->setVal(0.0); - pml_j_cp[2]->setVal(0.0); - amrex::Print() << "PML HAS PARTICLES - coarse"<< std::endl; - } + pml_j_cp[0].reset(new MultiFab(amrex::convert(cba,WarpX::jx_nodal_flag), cdm, 1, ngb)); + pml_j_cp[1].reset(new MultiFab(amrex::convert(cba,WarpX::jy_nodal_flag), cdm, 1, ngb)); + pml_j_cp[2].reset(new MultiFab(amrex::convert(cba,WarpX::jz_nodal_flag), cdm, 1, ngb)); + pml_j_cp[0]->setVal(0.0); + pml_j_cp[1]->setVal(0.0); + pml_j_cp[2]->setVal(0.0); + // if (pml_has_particles) + // { + // pml_j_cp[0].reset(new MultiFab(amrex::convert(cba,WarpX::jx_nodal_flag), cdm, 1, ngb)); + // pml_j_cp[1].reset(new MultiFab(amrex::convert(cba,WarpX::jy_nodal_flag), cdm, 1, ngb)); + // pml_j_cp[2].reset(new MultiFab(amrex::convert(cba,WarpX::jz_nodal_flag), cdm, 1, ngb)); + // pml_j_cp[0]->setVal(0.0); + // pml_j_cp[1]->setVal(0.0); + // pml_j_cp[2]->setVal(0.0); + // amrex::Print() << "PML HAS PARTICLES - coarse"<< std::endl; + // } // sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba_reduced, cgeom->CellSize(), ncell, delta)); @@ -714,6 +725,7 @@ PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom) const IntVect& ngp = pml.nGrowVect(); const int ncp = pml.nComp(); const auto& period = geom.periodicity(); + int ncells = 10; MultiFab tmpregmf(reg.boxArray(), reg.DistributionMap(), ncp, ngr); @@ -735,11 +747,19 @@ PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom) { const FArrayBox& src = tmpregmf[mfi]; FArrayBox& dst = reg[mfi]; - const BoxList& bl = amrex::boxDiff(dst.box(), mfi.validbox()); - for (const Box& bx : bl) - { - dst.copy(src, bx, 0, bx, 0, 1); - } + // BoxList + // Box box1 = dst.box(); + // amrex::grow(box1, -ncells); + // Box box2 = mfi.validbox(); + // amrex::Print() << "dst.box() = ["<< box1.smallEnd()[0] << ", "<ok()){ - pml[lev]->CopyFieldInReg({ Efield_fp[lev][0].get(), - Efield_fp[lev][1].get(), - Efield_fp[lev][2].get() }, - { Efield_cp[lev][0].get(), - Efield_cp[lev][1].get(), - Efield_cp[lev][2].get() }); - } - } - } + // if (do_pml) { + // // copy E from pmls to mothergrid + // for (int lev = 0; lev <= finest_level; ++lev) + // { + // if (pml[lev]->ok()){ + // pml[lev]->CopyFieldInReg({ Efield_fp[lev][0].get(), + // Efield_fp[lev][1].get(), + // Efield_fp[lev][2].get() }, + // { Efield_cp[lev][0].get(), + // Efield_cp[lev][1].get(), + // Efield_cp[lev][2].get() }); + // } + // } + // } EvolveF(0.5*dt[0], DtType::SecondHalf); EvolveB(0.5*dt[0]); // We now have B^{n+1} if (do_pml) { @@ -342,20 +341,20 @@ WarpX::OneStep_nosub (Real cur_time) FillBoundaryE(); } FillBoundaryB(); - if (do_pml) { - // copy B from pmls to mothergrid - for (int lev = 0; lev <= finest_level; ++lev) - { - if (pml[lev]->ok()){ - pml[lev]->CopyFieldInReg({ Bfield_fp[lev][0].get(), - Bfield_fp[lev][1].get(), - Bfield_fp[lev][2].get() }, - { Bfield_cp[lev][0].get(), - Bfield_cp[lev][1].get(), - Bfield_cp[lev][2].get() }); - } - } - } + // if (do_pml) { + // // copy B from pmls to mothergrid + // for (int lev = 0; lev <= finest_level; ++lev) + // { + // if (pml[lev]->ok()){ + // pml[lev]->CopyFieldInReg({ Bfield_fp[lev][0].get(), + // Bfield_fp[lev][1].get(), + // Bfield_fp[lev][2].get() }, + // { Bfield_cp[lev][0].get(), + // Bfield_cp[lev][1].get(), + // Bfield_cp[lev][2].get() }); + // } + // } + // } #endif } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 4b3be0a43..b3d7b8469 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -348,7 +348,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) &mu_c2_dt, #endif &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); - + if (pml_F) { WRPX_PUSH_PML_EVEC_F( -- cgit v1.2.3 From beec5d75193fc12145ee092ffea049640987692c Mon Sep 17 00:00:00 2001 From: ablelly Date: Thu, 4 Jul 2019 01:39:42 +0200 Subject: Modifications for particles in PMLs in 3d : added variables to push pmls and to damp current. Compiling in 3d. --- Source/BoundaryConditions/PML.H | 9 +++- Source/BoundaryConditions/PML.cpp | 26 +--------- Source/BoundaryConditions/PML_routines.F90 | 78 +++++++++++++++++++++++++++--- Source/FieldSolver/WarpXPushFieldsEM.cpp | 12 ++++- Source/FortranInterface/WarpX_f.H | 24 ++++++++- 5 files changed, 114 insertions(+), 35 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 2ff736bc0..97b0b6b6f 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -16,6 +16,14 @@ (x).sigma_star_fac[1].data(), (x).sigma_star_fac[1].m_lo, (x).sigma_star_fac[1].m_hi, \ (x).sigma_star_fac[2].data(), (x).sigma_star_fac[2].m_lo, (x).sigma_star_fac[2].m_hi +#define WRPX_PML_SIGMAJ_TO_FORTRAN(x) \ + (x).sigma_cum_fac[0].data(), (x).sigma_cum_fac[0].m_lo, (x).sigma_cum_fac[0].m_hi, \ + (x).sigma_cum_fac[1].data(), (x).sigma_cum_fac[1].m_lo, (x).sigma_cum_fac[1].m_hi, \ + (x).sigma_cum_fac[2].data(), (x).sigma_cum_fac[2].m_lo, (x).sigma_cum_fac[2].m_hi, \ + (x).sigma_star_cum_fac[0].data(), (x).sigma_star_cum_fac[0].m_lo, (x).sigma_star_cum_fac[0].m_hi, \ + (x).sigma_star_cum_fac[1].data(), (x).sigma_star_cum_fac[1].m_lo, (x).sigma_star_cum_fac[1].m_hi, \ + (x).sigma_star_cum_fac[2].data(), (x).sigma_star_cum_fac[2].m_lo, (x).sigma_star_cum_fac[2].m_hi + #else #define WRPX_PML_TO_FORTRAN(x) \ @@ -58,7 +66,6 @@ struct SigmaBox SigmaVect sigma_star_fac; SigmaVect sigma_star_cum_fac; - std::array alpha; //{alpha_xy, alpha_xz, alpha_xf, alpha_yx, alpha_yz, alpha_fy, alpha_zx, alpha_zy, alpha_fz} std::string pml_type; // if pml normal to x : if pml Hi : pml_type={1,0,0}; if pml Lo: pml_type={-1,0,0} }; diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index a4e95d74b..f9ef0b632 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -24,7 +24,6 @@ namespace int slo = sigma.m_lo; int shi = sigma.m_hi; int sslo = sigma_star.m_lo; - Real cumsum=0.; Real x = 10.0; Real theta = 10.0; Real coeff_damp = std::sin(theta*MathConst::pi/180.); @@ -53,11 +52,9 @@ namespace int slo = sigma.m_lo; int shi = sigma.m_hi; int sslo = sigma_star.m_lo; - Real cumsum = 0.; Real x = 10.0; Real theta = 10.0; Real coeff_damp = std::sin(theta*MathConst::pi/180.); - // amrex::Print()<<"===== FillHi =====" <(i-ghi-1); @@ -123,7 +120,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n sigma_star_cum_fac[idim].m_lo = lo[idim]; sigma_star_cum_fac[idim].m_hi = hi[idim]; } - alpha = {0.,0.,0.,0.,0.,0.,0.,0.,0.}; + pml_type = ""; Array fac; @@ -346,27 +343,6 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n } amrex::Print()<<"pml_type = "<GetE_fp() : pml[lev]->GetE_cp(); const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp(); const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); + const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp() + : pml[lev]->GetMultiSigmaBox_cp(); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -346,9 +348,17 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) BL_TO_FORTRAN_3D((*pml_j[2])[mfi]), &pml_has_particles, &mu_c2_dt, +#endif +#if (AMREX_SPACEDIM==3) + BL_TO_FORTRAN_3D((*pml_j[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_j[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_j[2])[mfi]), + &pml_has_particles, + WRPX_PML_SIGMAJ_TO_FORTRAN(sigba[mfi]), + &mu_c2_dt, #endif &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); - + if (pml_F) { WRPX_PUSH_PML_EVEC_F( diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index c9ef7cc11..4ca44b92f 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -38,6 +38,7 @@ #define WRPX_INTERPOLATE_CIC_TWO_LEVELS warpx_interpolate_cic_two_levels_3d #define WRPX_PUSH_LEAPFROG warpx_push_leapfrog_3d #define WRPX_PUSH_LEAPFROG_POSITIONS warpx_push_leapfrog_positions_3d +#define WRPX_DAMPJ_PML warpx_dampJ_pml_3d #elif (AMREX_SPACEDIM == 2) @@ -384,6 +385,19 @@ extern "C" BL_FORT_FAB_ARG_3D(jz), const int* flag, const amrex::Real* mudt, +#endif +#if (AMREX_SPACEDIM == 3) + BL_FORT_FAB_ARG_3D(jx), + BL_FORT_FAB_ARG_3D(jy), + BL_FORT_FAB_ARG_3D(jz), + const int* flag, + const amrex::Real* sigjx, int sigjx_lo, int sigjx_hi, + const amrex::Real* sigjy, int sigjy_lo, int sigjy_hi, + const amrex::Real* sigjz, int sigjz_lo, int sigjz_hi, + const amrex::Real* sigsjx, int sigsjx_lo, int sigsjx_hi, + const amrex::Real* sigsjy, int sigsjy_lo, int sigsjy_hi, + const amrex::Real* sigjsz, int sigsjz_lo, int sigsjz_hi, + const amrex::Real* mudt, #endif const amrex::Real* dtsdx, const amrex::Real* dtsdy, @@ -446,7 +460,7 @@ extern "C" #endif const amrex::Real* sigbz, int sigbz_lo, int sigbz_hi); -#if (AMREX_SPACEDIM == 2) + void WRPX_DAMPJ_PML (const int* tjxlo, const int* tjxhi, const int* tjylo, const int* tjyhi, const int* tjzlo, const int* tjzhi, @@ -454,10 +468,16 @@ extern "C" amrex::Real* jy, const int* jylo, const int* jyhi, amrex::Real* jz, const int* jzlo, const int* jzhi, const amrex::Real* sigjx, int sigjx_lo, int sigjx_hi, +#if (AMREX_SPACEDIM == 3) + const amrex::Real* sigjy, int sigjy_lo, int sigjy_hi, +#endif const amrex::Real* sigjz, int sigjz_lo, int sigjz_hi, const amrex::Real* sigsjx, int sigsjx_lo, int sigsjx_hi, - const amrex::Real* sigjsz, int sigsjz_lo, int sigsjz_hi); +#if (AMREX_SPACEDIM == 3) + const amrex::Real* sigsjy, int sigsjy_lo, int sigsjy_hi, #endif + const amrex::Real* sigjsz, int sigsjz_lo, int sigsjz_hi); + void WRPX_SYNC_CURRENT (const int* lo, const int* hi, BL_FORT_FAB_ARG_ANYD(crse), -- cgit v1.2.3 From 40e9b4387e6d617358290bda29124ee2eae1e1ec Mon Sep 17 00:00:00 2001 From: ablelly Date: Fri, 5 Jul 2019 19:04:24 +0200 Subject: Started to implement push_pml_evec_3d for all different PMLs. --- Source/BoundaryConditions/PML.H | 8 + Source/BoundaryConditions/PML_routines.F90 | 334 ++++++++++++++++++++++++--- Source/BoundaryConditions/WarpXEvolvePML.cpp | 2 +- Source/FieldSolver/WarpXPushFieldsEM.cpp | 1 + Source/FortranInterface/WarpX_f.H | 1 + 5 files changed, 312 insertions(+), 34 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 97b0b6b6f..ef1682f92 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -24,6 +24,14 @@ (x).sigma_star_cum_fac[1].data(), (x).sigma_star_cum_fac[1].m_lo, (x).sigma_star_cum_fac[1].m_hi, \ (x).sigma_star_cum_fac[2].data(), (x).sigma_star_cum_fac[2].m_lo, (x).sigma_star_cum_fac[2].m_hi +#define WRPX_PML_SIGMACOEFF_TO_FORTRAN(x) \ + (x).sigma[0].data(), (x).sigma[0].m_lo, (x).sigma[0].m_hi, \ + (x).sigma[1].data(), (x).sigma[1].m_lo, (x).sigma[1].m_hi, \ + (x).sigma[2].data(), (x).sigma[2].m_lo, (x).sigma[2].m_hi, \ + (x).sigma_star[0].data(), (x).sigma_star[0].m_lo, (x).sigma_star[0].m_hi, \ + (x).sigma_star[1].data(), (x).sigma_star[1].m_lo, (x).sigma_star[1].m_hi, \ + (x).sigma_star[2].data(), (x).sigma_star[2].m_lo, (x).sigma_star[2].m_hi + #else #define WRPX_PML_TO_FORTRAN(x) \ diff --git a/Source/BoundaryConditions/PML_routines.F90 b/Source/BoundaryConditions/PML_routines.F90 index 39aa707d3..c17889062 100644 --- a/Source/BoundaryConditions/PML_routines.F90 +++ b/Source/BoundaryConditions/PML_routines.F90 @@ -259,6 +259,7 @@ contains & jy, jylo, jyhi, & & jz, jzlo, jzhi, & & flag, & + & pml_type, & & sigjx, sjxlo, sjxhi, & & sigjy, sjylo, sjyhi, & & sigjz, sjzlo, sjzhi, & @@ -275,7 +276,7 @@ contains sjxlo, sjxhi, sjylo, sjyhi, sjzlo, sjzhi, & ssjxlo, ssjxhi, ssjylo, ssjyhi, ssjzlo, ssjzhi, & flag - real(amrex_real), intent(in) :: mudt, dtsdx, dtsdy, dtsdz + real(amrex_real), intent(in ) :: mudt, dtsdx, dtsdy, dtsdz real(amrex_real), intent(inout) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),Exlo(3):Exhi(3),2) real(amrex_real), intent(inout) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),Eylo(3):Eyhi(3),2) real(amrex_real), intent(inout) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),Ezlo(3):Ezhi(3),2) @@ -291,8 +292,11 @@ contains real(amrex_real), intent(in ) :: sigsjx (ssjxlo:ssjxhi) real(amrex_real), intent(in ) :: sigsjy (ssjylo:ssjyhi) real(amrex_real), intent(in ) :: sigsjz (ssjzlo:ssjzhi) + Character, intent(in ) :: pml_type (9) integer :: i, j, k + real(amrex_real) :: sum_sigma + real(amrex_real) :: alpha_xy, alpha_xz, alpha_yx, alpha_yz, alpha_zx, alpha_zy if (flag == 0) then do k = xlo(3), xhi(3) @@ -328,40 +332,304 @@ contains end do end do - else - do k = xlo(3), xhi(3) - do j = xlo(2), xhi(2) - do i = xlo(1), xhi(1) - 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)) - Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & - & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) - end do - end do - end do + else !flag = 1 + !!!!! CORNER + if (pml_type(1)=='c') then + do k = xlo(3), xhi(3) + do j = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + ! compute current coefficients alpha - do k = ylo(3), yhi(3) - do j = ylo(2), yhi(2) - do i = ylo(1), yhi(1) - 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)) - Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & - & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) - end do - end do - end do + 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)) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) + end do + end do + end do - do k = zlo(3), zhi(3) - do j = zlo(2), zhi(2) - do i = zlo(1), zhi(1) - 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)) - Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & - & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) - end do - end do - end do - end if + do k = ylo(3), yhi(3) + do j = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + 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)) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) + end do + end do + end do + + do k = zlo(3), zhi(3) + do j = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + 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)) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) + end do + end do + end do + + + !!!!! DIRECT FACE + else if (pml_type(1)=='d') then + if (pml_type(3)=='0') then + do k = xlo(3), xhi(3) + do j = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + ! compute current coefficients alpha + + 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)) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) + end do + end do + end do + + do k = ylo(3), yhi(3) + do j = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + 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)) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) + end do + end do + end do + + do k = zlo(3), zhi(3) + do j = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + 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)) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) + end do + end do + end do + + else if (pml_type(3)=='1') then + do k = xlo(3), xhi(3) + do j = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + ! compute current coefficients alpha + + 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)) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) + end do + end do + end do + + do k = ylo(3), yhi(3) + do j = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + 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)) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) + end do + end do + end do + + do k = zlo(3), zhi(3) + do j = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + 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)) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) + end do + end do + end do + + else !(pml_type(3)=='2') + do k = xlo(3), xhi(3) + do j = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + ! compute current coefficients alpha + + 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)) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) + end do + end do + end do + + do k = ylo(3), yhi(3) + do j = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + 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)) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) + end do + end do + end do + + do k = zlo(3), zhi(3) + do j = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + 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)) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) + end do + end do + end do + end if !DIRECT FACE + + !!!!! SIDE EDGE + else !(pml_type(1)=='s') + if ((pml_type(3)=='0') .AND. (pml_type(4)=='1')) then + do k = xlo(3), xhi(3) + do j = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + ! compute current coefficients alpha + + 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)) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) + end do + end do + end do + + do k = ylo(3), yhi(3) + do j = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + 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)) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) + end do + end do + end do + + do k = zlo(3), zhi(3) + do j = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + 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)) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) + end do + end do + end do + else if ((pml_type(3)=='0') .AND. (pml_type(4)=='2')) then + do k = xlo(3), xhi(3) + do j = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + ! compute current coefficients alpha + + 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)) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) + end do + end do + end do + + do k = ylo(3), yhi(3) + do j = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + 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)) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) + end do + end do + end do + + do k = zlo(3), zhi(3) + do j = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + 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)) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) + end do + end do + end do + + else !((pml_type(3)=='1') .AND. (pml_type(4)=='2')) + do k = xlo(3), xhi(3) + do j = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + ! compute current coefficients alpha + + 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)) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) + end do + end do + end do + + do k = ylo(3), yhi(3) + do j = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + 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)) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) + end do + end do + end do + + do k = zlo(3), zhi(3) + do j = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + 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)) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) + end do + end do + end do + end if ! SIDE EDGE + + end if !PML_TYPE + + end if ! FLAG + ! do k = xlo(3), xhi(3) + ! do j = xlo(2), xhi(2) + ! do i = xlo(1), xhi(1) + ! ! compute current coefficients alpha + ! + ! 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)) + ! Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + ! & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) + ! end do + ! end do + ! end do + ! + ! do k = ylo(3), yhi(3) + ! do j = ylo(2), yhi(2) + ! do i = ylo(1), yhi(1) + ! 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)) + ! Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + ! & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) + ! end do + ! end do + ! end do + ! + ! do k = zlo(3), zhi(3) + ! do j = zlo(2), zhi(2) + ! do i = zlo(1), zhi(1) + ! 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)) + ! Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + ! & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) + ! end do + ! end do + ! end do + ! end if end subroutine warpx_push_pml_evec_3d diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp index 8050f4b4e..5f3057514 100644 --- a/Source/BoundaryConditions/WarpXEvolvePML.cpp +++ b/Source/BoundaryConditions/WarpXEvolvePML.cpp @@ -125,7 +125,7 @@ WarpX::DampJPML (int lev, PatchType patch_type) BL_TO_FORTRAN_3D((*pml_j[0])[mfi]), BL_TO_FORTRAN_3D((*pml_j[1])[mfi]), BL_TO_FORTRAN_3D((*pml_j[2])[mfi]), - WRPX_PML_SIGMAJ_TO_FORTRAN(sigba[mfi])); + WRPX_PML_SIGMACOEFF_TO_FORTRAN(sigba[mfi])); //WRPX_PML_SIGMAJ_TO_FORTRAN(sigba[mfi]) } } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 322b38611..38dfeb3f6 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -354,6 +354,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) BL_TO_FORTRAN_3D((*pml_j[1])[mfi]), BL_TO_FORTRAN_3D((*pml_j[2])[mfi]), &pml_has_particles, + sigba[mfi].pml_type, WRPX_PML_SIGMAJ_TO_FORTRAN(sigba[mfi]), &mu_c2_dt, #endif diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 4ca44b92f..c505b598a 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -391,6 +391,7 @@ extern "C" BL_FORT_FAB_ARG_3D(jy), BL_FORT_FAB_ARG_3D(jz), const int* flag, + std::string pml_type, const amrex::Real* sigjx, int sigjx_lo, int sigjx_hi, const amrex::Real* sigjy, int sigjy_lo, int sigjy_hi, const amrex::Real* sigjz, int sigjz_lo, int sigjz_hi, -- cgit v1.2.3 From cce601387d0e0c346d02a63ceb62b387f69e89ae Mon Sep 17 00:00:00 2001 From: ablelly Date: Fri, 5 Jul 2019 23:17:48 +0200 Subject: Added damp J in 3d, corrected dampJ in 2d. --- Source/BoundaryConditions/PML.H | 6 ++++ Source/BoundaryConditions/PML_routines.F90 | 48 +++++++++++++++------------- Source/BoundaryConditions/WarpXEvolvePML.cpp | 2 +- Source/FieldSolver/WarpXPushFieldsEM.cpp | 2 +- 4 files changed, 34 insertions(+), 24 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index ef1682f92..20f7ee4d7 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -46,6 +46,12 @@ (x).sigma_star_cum_fac[0].data(), (x).sigma_star_cum_fac[0].m_lo, (x).sigma_star_cum_fac[0].m_hi, \ (x).sigma_star_cum_fac[1].data(), (x).sigma_star_cum_fac[1].m_lo, (x).sigma_star_cum_fac[1].m_hi +#define WRPX_PML_SIGMACOEFF_TO_FORTRAN(x) \ + (x).sigma[0].data(), (x).sigma[0].m_lo, (x).sigma[0].m_hi, \ + (x).sigma[1].data(), (x).sigma[1].m_lo, (x).sigma[1].m_hi, \ + (x).sigma_star[0].data(), (x).sigma_star[0].m_lo, (x).sigma_star[0].m_hi, \ + (x).sigma_star[1].data(), (x).sigma_star[1].m_lo, (x).sigma_star[1].m_hi + #endif struct Sigma : amrex::Vector diff --git a/Source/BoundaryConditions/PML_routines.F90 b/Source/BoundaryConditions/PML_routines.F90 index 5bb08e693..162fd1ed2 100644 --- a/Source/BoundaryConditions/PML_routines.F90 +++ b/Source/BoundaryConditions/PML_routines.F90 @@ -295,7 +295,6 @@ contains Character, intent(in ) :: pml_type (9) integer :: i, j, k - real(amrex_real) :: sum_sigma real(amrex_real) :: alpha_xy, alpha_xz, alpha_yx, alpha_yz, alpha_zx, alpha_zy if (flag == 0) then @@ -709,7 +708,7 @@ contains end do end do - else + else !no particles in PML do k = xlo(2), xhi(2) do i = xlo(1), xhi(1) Ex(i,k,2) = Ex(i,k,2) - dtsdz*(By(i,k ,1)+By(i,k ,2) & @@ -732,7 +731,8 @@ contains & -By(i-1,k,1)-By(i-1,k,2)) end do end do - end if + + end if !flag end subroutine warpx_push_pml_evec_2d @@ -1154,7 +1154,7 @@ contains do k = tjxlo(2), tjxhi(2) do i = tjxlo(1), tjxhi(1) ! jx(i,k) = jx(i,k) * sigjx(i) !minval((/sigjx(i),sigjz(k)/)) - jx(i,k) = jx(i,k) * sigsjx(i) + jx(i,k) = jx(i,k) * sigsjx(i) * sigjz(k) end do end do @@ -1166,7 +1166,7 @@ contains do k = tjzlo(2), tjzhi(2) do i = tjzlo(1), tjzhi(1) - jz(i,k) = jz(i,k) * sigjx(i) + jz(i,k) = jz(i,k) * sigjx(i) * sigsjz(k) end do end do @@ -1178,12 +1178,12 @@ contains & sigjx, sjxlo, sjxhi, sigjy, sjylo, sjyhi, sigjz, sjzlo, sjzhi, & & sigsjx, ssjxlo, ssjxhi, sigsjy, ssjylo, ssjyhi, sigsjz, ssjzlo, ssjzhi) & bind(c,name='warpx_dampJ_pml_3d') - integer, dimension(2), intent(in) :: tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & + integer, dimension(3), intent(in) :: tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & jxlo, jxhi, jylo, jyhi, jzlo, jzhi integer, intent(in), value :: sjxlo, sjxhi, sjylo, sjyhi, sjzlo, sjzhi, ssjxlo, ssjxhi, ssjylo, ssjyhi, ssjzlo, ssjzhi - real(amrex_real), intent(inout) :: jx(jxlo(1):jxhi(1),jxlo(2):jxhi(2)) - real(amrex_real), intent(inout) :: jy(jylo(1):jyhi(1),jylo(2):jyhi(2)) - real(amrex_real), intent(inout) :: jz(jzlo(1):jzhi(1),jzlo(2):jzhi(2)) + real(amrex_real), intent(inout) :: jx(jxlo(1):jxhi(1),jxlo(2):jxhi(2),jxlo(3):jxhi(3)) + real(amrex_real), intent(inout) :: jy(jylo(1):jyhi(1),jylo(2):jyhi(2),jylo(3):jyhi(3)) + real(amrex_real), intent(inout) :: jz(jzlo(1):jzhi(1),jzlo(2):jzhi(2),jzlo(3):jzhi(3)) real(amrex_real), intent(in) :: sigjx(sjxlo:sjxhi) real(amrex_real), intent(in) :: sigjy(sjylo:sjyhi) real(amrex_real), intent(in) :: sigjz(sjzlo:sjzhi) @@ -1191,28 +1191,32 @@ contains real(amrex_real), intent(in) :: sigsjy(ssjylo:ssjyhi) real(amrex_real), intent(in) :: sigsjz(ssjzlo:ssjzhi) - integer :: i,k - !!!! FOR A PML ALONG X-AXIS !!!!! - do k = tjxlo(2), tjxhi(2) - do i = tjxlo(1), tjxhi(1) - ! jx(i,k) = jx(i,k) * sigjx(i) !minval((/sigjx(i),sigjz(k)/)) - jx(i,k) = jx(i,k) * sigsjx(i) + integer :: i,j,k + + do k = tjxlo(3), tjxhi(3) + do j = tjylo(2), tjyhi(2) + do i = tjxlo(1), tjxhi(1) + jx(i,j,k) = jx(i,j,k) * sigsjx(i) * sigjy(j) * sigjz(k) + end do end do end do - do k = tjylo(2), tjyhi(2) - do i = tjylo(1), tjyhi(1) - jy(i,k) = jy(i,k) !* minval((/sigjx(i),sigjz(k)/)) !sigjz(k) !no current jy... + do k = tjxlo(3), tjxhi(3) + do j = tjylo(2), tjyhi(2) + do i = tjxlo(1), tjxhi(1) + jy(i,j,k) = jy(i,j,k) * sigjx(i) * sigsjy(j) * sigjz(k) + end do end do end do - do k = tjzlo(2), tjzhi(2) - do i = tjzlo(1), tjzhi(1) - jz(i,k) = jz(i,k) * sigjx(i) + do k = tjxlo(3), tjxhi(3) + do j = tjylo(2), tjyhi(2) + do i = tjxlo(1), tjxhi(1) + jz(i,j,k) = jz(i,j,k) * sigjx(i) * sigjy(j) * sigsjz(k) + end do end do end do - end subroutine warpx_dampJ_pml_3d diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp index 5f3057514..08173e424 100644 --- a/Source/BoundaryConditions/WarpXEvolvePML.cpp +++ b/Source/BoundaryConditions/WarpXEvolvePML.cpp @@ -125,7 +125,7 @@ WarpX::DampJPML (int lev, PatchType patch_type) BL_TO_FORTRAN_3D((*pml_j[0])[mfi]), BL_TO_FORTRAN_3D((*pml_j[1])[mfi]), BL_TO_FORTRAN_3D((*pml_j[2])[mfi]), - WRPX_PML_SIGMACOEFF_TO_FORTRAN(sigba[mfi])); //WRPX_PML_SIGMAJ_TO_FORTRAN(sigba[mfi]) + WRPX_PML_SIGMAJ_TO_FORTRAN(sigba[mfi])); //WRPX_PML_SIGMACOEFF_TO_FORTRAN(sigba[mfi]) } } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 38dfeb3f6..afb979996 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -355,7 +355,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) BL_TO_FORTRAN_3D((*pml_j[2])[mfi]), &pml_has_particles, sigba[mfi].pml_type, - WRPX_PML_SIGMAJ_TO_FORTRAN(sigba[mfi]), + WRPX_PML_SIGMACOEFF_TO_FORTRAN(sigba[mfi]), //WRPX_PML_SIGMAJ_TO_FORTRAN( &mu_c2_dt, #endif &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); -- cgit v1.2.3 From 916f8744b353930add6572f711ae9d7a08260d97 Mon Sep 17 00:00:00 2001 From: ablelly Date: Tue, 9 Jul 2019 17:27:16 +0200 Subject: Damp PML modified (corrected the entries) --- Source/BoundaryConditions/PML.H | 13 +++---- Source/BoundaryConditions/PML.cpp | 54 +++++++++++++++++++++++++----- Source/BoundaryConditions/PML_routines.F90 | 15 +++++---- Source/Evolve/WarpXEvolveEM.cpp | 7 +++- Source/FieldSolver/WarpXPushFieldsEM.cpp | 14 ++++++-- Source/FortranInterface/WarpX_f.H | 14 ++++---- 6 files changed, 86 insertions(+), 31 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 20f7ee4d7..39f99dfef 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -25,12 +25,12 @@ (x).sigma_star_cum_fac[2].data(), (x).sigma_star_cum_fac[2].m_lo, (x).sigma_star_cum_fac[2].m_hi #define WRPX_PML_SIGMACOEFF_TO_FORTRAN(x) \ - (x).sigma[0].data(), (x).sigma[0].m_lo, (x).sigma[0].m_hi, \ - (x).sigma[1].data(), (x).sigma[1].m_lo, (x).sigma[1].m_hi, \ - (x).sigma[2].data(), (x).sigma[2].m_lo, (x).sigma[2].m_hi, \ - (x).sigma_star[0].data(), (x).sigma_star[0].m_lo, (x).sigma_star[0].m_hi, \ - (x).sigma_star[1].data(), (x).sigma_star[1].m_lo, (x).sigma_star[1].m_hi, \ - (x).sigma_star[2].data(), (x).sigma_star[2].m_lo, (x).sigma_star[2].m_hi + (x).sigma[0].dataPtr(), &((x).sigma[0].m_lo), &((x).sigma[0].m_hi), \ + (x).sigma[1].dataPtr(), &((x).sigma[1].m_lo), &((x).sigma[1].m_hi), \ + (x).sigma[2].dataPtr(), &((x).sigma[2].m_lo), &((x).sigma[2].m_hi), \ + (x).sigma_star[0].dataPtr(), &((x).sigma_star[0].m_lo), &((x).sigma_star[0].m_hi), \ + (x).sigma_star[1].dataPtr(), &((x).sigma_star[1].m_lo), &((x).sigma_star[1].m_hi), \ + (x).sigma_star[2].dataPtr(), &((x).sigma_star[2].m_lo), &((x).sigma_star[2].m_hi) #else @@ -81,6 +81,7 @@ struct SigmaBox SigmaVect sigma_star_cum_fac; std::string pml_type; // if pml normal to x : if pml Hi : pml_type={1,0,0}; if pml Lo: pml_type={-1,0,0} + std::array pml_type_array; }; namespace amrex { diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index f9ef0b632..ec77ca7cf 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -26,7 +26,7 @@ namespace int sslo = sigma_star.m_lo; Real x = 10.0; Real theta = 10.0; - Real coeff_damp = std::sin(theta*MathConst::pi/180.); + Real coeff_damp = 1.; //std::sin(theta*MathConst::pi/180.); for (int i = olo; i <= ohi+1; ++i) { @@ -54,7 +54,7 @@ namespace int sslo = sigma_star.m_lo; Real x = 10.0; Real theta = 10.0; - Real coeff_damp = std::sin(theta*MathConst::pi/180.); + Real coeff_damp = 1.;//std::sin(theta*MathConst::pi/180.); for (int i = olo; i <= ohi+1; ++i) { Real offset = static_cast(i-ghi-1); @@ -85,7 +85,8 @@ namespace SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int ncell, int delta) { amrex::Print()<< "===== BUILD SigmaBox ====="< fac; for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { @@ -198,10 +200,15 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n FillLo(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], looverlap, grid_box, fac[idim]); if (idim == 0){ pml_type = "crn-++"; // for 2d only use the two first components + pml_type_array[0] = 1; + pml_type_array[1] = 0; + pml_type_array[2] = 1; + pml_type_array[3] = 1; } else { if (pml_type[0] == 'c'){ pml_type[3+idim] = '-'; + pml_type_array[1+idim] = 0; } } } @@ -216,10 +223,15 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n FillHi(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], hioverlap, grid_box, fac[idim]); if (idim == 0){ pml_type = "crn+++"; // for 2d only use the two first components + pml_type_array[0] = 1; + pml_type_array[1] = 1; + pml_type_array[2] = 1; + pml_type_array[3] = 1; } else { if (pml_type[0] == 'c'){ pml_type[3+idim] = '+'; + pml_type_array[1+idim] = 1; } } } @@ -238,6 +250,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n FillZero(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], overlap); if (idim==0){ pml_type="sed33--"; + pml_type_array[0] = 2; } } else { @@ -255,16 +268,23 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n FillLo(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], looverlap, grid_box, fac[idim]); if (idim == 0){ pml_type = "sed03-+"; // for 2d only use the two first components + pml_type_array[0] = 2; + pml_type_array[1] = idim; + pml_type_array[3] = 0; } else { if (pml_type[0] == 's'){ - if (pml_type[3]==3){ + if (pml_type[3]=='3'){ pml_type[3]=idim+'0'; pml_type[5]='-'; + pml_type_array[1] = idim; + pml_type_array[3] = 0; } else{ pml_type[4] = idim+'0'; pml_type[6] = '-'; + pml_type_array[2] = idim; + pml_type_array[4] = 0; } } } @@ -276,16 +296,23 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n FillHi(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], hioverlap, grid_box, fac[idim]); if (idim == 0){ pml_type = "sed03++"; // for 2d only use the two first components + pml_type_array[0] = 2; + pml_type_array[1] = idim; + pml_type_array[3] = 1; } else { if (pml_type[0] == 's'){ - if (pml_type[3]==3){ + if (pml_type[3]=='3'){ pml_type[3]=idim+'0'; pml_type[5]='+'; + pml_type_array[1] = idim; + pml_type_array[3] = 1; } else{ pml_type[4] = idim+'0'; pml_type[6] = '+'; + pml_type_array[2] = idim; + pml_type_array[4] = 1; } } } @@ -321,6 +348,9 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n if (looverlap.ok()) { FillLo(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], looverlap, grid_box, fac[idim]); pml_type = "df0-"; + pml_type_array[0] = 0; + pml_type_array[1] = idim; + pml_type_array[2] = 0; if (idim > 0){pml_type[2]=idim+'0';} //std::to_string(idim) } @@ -329,6 +359,9 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n if (hioverlap.ok()) { FillHi(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], hioverlap, grid_box, fac[idim]); pml_type = "df0+"; + pml_type_array[0] = 0; + pml_type_array[1] = idim; + pml_type_array[2] = 1; if (idim>0){pml_type[2] = idim+'0';} } @@ -342,7 +375,11 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n } } amrex::Print()<<"pml_type = "<Domain(), -ncell); + amrex::Print() << "[" << domain0.smallEnd()[0]<<", "<< domain0.smallEnd()[1]<<", "<< domain0.smallEnd()[2]<< ", "<ok()){ @@ -316,6 +317,7 @@ WarpX::OneStep_nosub (Real cur_time) } if (do_pml && do_pml_j_damping){ + amrex::Print()<< "WarpXEvolveEM.cpp : DampJ "<ok()){ @@ -357,6 +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 "<Getj_fp() : pml[lev]->Getj_cp(); const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp(); const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp() - : pml[lev]->GetMultiSigmaBox_cp(); + : pml[lev]->GetMultiSigmaBox_cp(); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -332,6 +332,14 @@ 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"<* pml_type,//const std::string pml_type, + const amrex::Real* sigjx, const int* sigjx_lo, const int* sigjx_hi, + const amrex::Real* sigjy, const int* sigjy_lo, const int* sigjy_hi, + const amrex::Real* sigjz, const int* sigjz_lo, const int* sigjz_hi, + const amrex::Real* sigsjx, const int* sigsjx_lo, const int* sigsjx_hi, + const amrex::Real* sigsjy, const int* sigsjy_lo, const int* sigsjy_hi, + const amrex::Real* sigjsz, const int* sigsjz_lo, const int* sigsjz_hi, const amrex::Real* mudt, #endif const amrex::Real* dtsdx, -- cgit v1.2.3 From e5c40e77957c0b3b8b8d85f68c75b847fce56a91 Mon Sep 17 00:00:00 2001 From: ablelly Date: Tue, 9 Jul 2019 19:30:35 +0200 Subject: Erased the prints. 3d version compiling and running. --- Source/BoundaryConditions/PML_routines.F90 | 45 ++++++++++++++++++++++++------ Source/Evolve/WarpXEvolveEM.cpp | 12 ++++---- Source/FieldSolver/WarpXPushFieldsEM.cpp | 13 +++++---- 3 files changed, 50 insertions(+), 20 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') 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 "<ok()){ @@ -317,7 +317,7 @@ WarpX::OneStep_nosub (Real cur_time) } if (do_pml && do_pml_j_damping){ - amrex::Print()<< "WarpXEvolveEM.cpp : DampJ "<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 "< Date: Tue, 16 Jul 2019 18:29:37 +0200 Subject: Started to add flag do_pml_in_domain for the merge request. --- Source/BoundaryConditions/PML.H | 18 +++---- Source/BoundaryConditions/PML.cpp | 82 +++++++++++++++++++------------- Source/FieldSolver/WarpXPushFieldsEM.cpp | 2 +- Source/Initialization/WarpXInitData.cpp | 6 +-- Source/Parallelization/WarpXComm.cpp | 18 +++---- 5 files changed, 71 insertions(+), 55 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 62a378986..c64eb8193 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -81,7 +81,7 @@ struct SigmaBox SigmaVect sigma_star_cum_fac; std::string pml_type; // if pml normal to x : if pml Hi : pml_type={1,0,0}; if pml Lo: pml_type={-1,0,0} - std::array pml_type_array; + std::array pml_type_array; }; namespace amrex { @@ -128,7 +128,7 @@ class PML public: PML (const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, const amrex::Geometry* geom, const amrex::Geometry* cgeom, - int ncell, int delta, int ref_ratio, int do_dive_cleaning, int do_moving_window, int pml_has_particles); + int ncell, int delta, int ref_ratio, int do_dive_cleaning, int do_moving_window, int pml_has_particles, int do_pml_in_domain); void ComputePMLFactors (amrex::Real dt); @@ -149,25 +149,25 @@ public: { return *sigba_cp; } void ExchangeB (const std::array& B_fp, - const std::array& B_cp); + const std::array& B_cp, int do_pml_in_domain); void ExchangeE (const std::array& E_fp, - const std::array& E_cp); + const std::array& E_cp, int do_pml_in_domain); void CopyJinPMLs (const std::array& j_fp, const std::array& j_cp); void CopyJinReg (const std::array& j_fp, const std::array& j_cp); void ExchangeB (PatchType patch_type, - const std::array& Bp); + const std::array& Bp, int do_pml_in_domain); void ExchangeE (PatchType patch_type, - const std::array& Ep); + const std::array& Ep, int do_pml_in_domain); void CopyJinPMLs (PatchType patch_type, const std::array& jp); void CopyJinReg (PatchType patch_type, const std::array& jp); - void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp); - void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp); + void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp, int do_pml_in_domain); + void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp, int do_pml_in_domain); void FillBoundary (); void FillBoundaryE (); @@ -205,7 +205,7 @@ private: static amrex::BoxArray MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell); - static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); + static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom, int do_pml_in_domain); static void CopyRegInPMLs (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); static void CopyPMLsInReg (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); }; diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 1f1cc80c3..c8d8585d8 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -474,10 +474,11 @@ MultiSigmaBox::ComputePMLFactorsE (const Real* dx, Real dt) PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, const Geometry* geom, const Geometry* cgeom, - int ncell, int delta, int ref_ratio, int do_dive_cleaning, int do_moving_window, int pml_has_particles) + int ncell, int delta, int ref_ratio, int do_dive_cleaning, int do_moving_window, int pml_has_particles, int do_pml_in_domain) : m_geom(geom), m_cgeom(cgeom) { + Box domain0 = geom->Domain(); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { if ( ! geom->isPeriodic(idim) ) { @@ -486,7 +487,11 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, } const BoxArray grid_ba_reduced = BoxArray(grid_ba.boxList().intersect(domain0)); - const BoxArray& ba = MakeBoxArray(*geom, grid_ba_reduced, ncell); //MakeBoxArray(*geom, grid_ba, ncell); + // const BoxArray& ba = MakeBoxArray(*geom, grid_ba_reduced, ncell); + // + // const BoxArray& ba = MakeBoxArray(*geom, grid_ba, ncell); + const BoxArray& ba = (do_pml_in_domain)? MakeBoxArray(*geom, grid_ba_reduced, ncell) : MakeBoxArray(*geom, grid_ba, ncell); + if (ba.size() == 0) { m_ok = false; return; @@ -529,8 +534,13 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, pml_F_fp->setVal(0.0); } + if (do_pml_in_domain){ + sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba_reduced, geom->CellSize(), ncell, delta)); + } + else { + sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba, geom->CellSize(), ncell, delta)); + } - sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba_reduced, geom->CellSize(), ncell, delta)); if (cgeom) { @@ -541,7 +551,8 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, BoxArray grid_cba = grid_ba; grid_cba.coarsen(ref_ratio); const BoxArray grid_cba_reduced = BoxArray(grid_cba.boxList().intersect(domain0)); - const BoxArray& cba = MakeBoxArray(*cgeom, grid_cba_reduced, ncell); + // const BoxArray& cba = MakeBoxArray(*cgeom, grid_cba_reduced, ncell); + const BoxArray& cba = (do_pml_in_domain) ? MakeBoxArray(*cgeom, grid_cba_reduced, ncell) : MakeBoxArray(*cgeom, grid_cba, ncell); DistributionMapping cdm{cba}; @@ -572,8 +583,13 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, pml_j_cp[1]->setVal(0.0); pml_j_cp[2]->setVal(0.0); - // sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); - sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba_reduced, cgeom->CellSize(), ncell, delta)); + if (do_pml_in_domain){ + sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba_reduced, cgeom->CellSize(), ncell, delta)); + } + else { + sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); + } + } } @@ -701,53 +717,53 @@ PML::GetF_cp () void PML::ExchangeB (const std::array& B_fp, - const std::array& B_cp) + const std::array& B_cp, int do_pml_in_domain) { - ExchangeB(PatchType::fine, B_fp); - ExchangeB(PatchType::coarse, B_cp); + ExchangeB(PatchType::fine, B_fp, do_pml_in_domain); + ExchangeB(PatchType::coarse, B_cp, do_pml_in_domain); } void PML::ExchangeB (PatchType patch_type, - const std::array& Bp) + const std::array& Bp, int do_pml_in_domain) { if (patch_type == PatchType::fine && pml_B_fp[0] && Bp[0]) { - Exchange(*pml_B_fp[0], *Bp[0], *m_geom); - Exchange(*pml_B_fp[1], *Bp[1], *m_geom); - Exchange(*pml_B_fp[2], *Bp[2], *m_geom); + Exchange(*pml_B_fp[0], *Bp[0], *m_geom, do_pml_in_domain); + Exchange(*pml_B_fp[1], *Bp[1], *m_geom, do_pml_in_domain); + Exchange(*pml_B_fp[2], *Bp[2], *m_geom, do_pml_in_domain); } else if (patch_type == PatchType::coarse && pml_B_cp[0] && Bp[0]) { - Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom); - Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom); - Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom); + Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom, do_pml_in_domain); + Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom, do_pml_in_domain); + Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom, do_pml_in_domain); } } void PML::ExchangeE (const std::array& E_fp, - const std::array& E_cp) + const std::array& E_cp, int do_pml_in_domain) { - ExchangeE(PatchType::fine, E_fp); - ExchangeE(PatchType::coarse, E_cp); + ExchangeE(PatchType::fine, E_fp, do_pml_in_domain); + ExchangeE(PatchType::coarse, E_cp, do_pml_in_domain); } void PML::ExchangeE (PatchType patch_type, - const std::array& Ep) + const std::array& Ep, int do_pml_in_domain) { if (patch_type == PatchType::fine && pml_E_fp[0] && Ep[0]) { - Exchange(*pml_E_fp[0], *Ep[0], *m_geom); - Exchange(*pml_E_fp[1], *Ep[1], *m_geom); - Exchange(*pml_E_fp[2], *Ep[2], *m_geom); + Exchange(*pml_E_fp[0], *Ep[0], *m_geom, do_pml_in_domain); + Exchange(*pml_E_fp[1], *Ep[1], *m_geom, do_pml_in_domain); + Exchange(*pml_E_fp[2], *Ep[2], *m_geom, do_pml_in_domain); } else if (patch_type == PatchType::coarse && pml_E_cp[0] && Ep[0]) { - Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom); - Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom); - Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom); + Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom, do_pml_in_domain); + Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom, do_pml_in_domain); + Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom, do_pml_in_domain); } } @@ -804,19 +820,19 @@ PML::CopyJinReg (const std::array& j_fp, } void -PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp) +PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp, int do_pml_in_domain) { - ExchangeF(PatchType::fine, F_fp); - ExchangeF(PatchType::coarse, F_cp); + ExchangeF(PatchType::fine, F_fp, do_pml_in_domain); + ExchangeF(PatchType::coarse, F_cp, do_pml_in_domain); } void -PML::ExchangeF (PatchType patch_type, MultiFab* Fp) +PML::ExchangeF (PatchType patch_type, MultiFab* Fp, int do_pml_in_domain) { if (patch_type == PatchType::fine && pml_F_fp && Fp) { - Exchange(*pml_F_fp, *Fp, *m_geom); + Exchange(*pml_F_fp, *Fp, *m_geom, do_pml_in_domain); } else if (patch_type == PatchType::coarse && pml_F_cp && Fp) { - Exchange(*pml_F_cp, *Fp, *m_cgeom); + Exchange(*pml_F_cp, *Fp, *m_cgeom, do_pml_in_domain); } } @@ -864,7 +880,7 @@ PML::ExchangeF (PatchType patch_type, MultiFab* Fp) // } void -PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom) +PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom, int do_pml_in_domain) { const IntVect& ngr = reg.nGrowVect(); diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 298b03dc6..8b749fc54 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -315,7 +315,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) if (do_pml && pml[lev]->ok()) { - if (F) pml[lev]->ExchangeF(patch_type, F); + if (F) pml[lev]->ExchangeF(patch_type, F, do_pml_in_domain); const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 76704bcfe..5598b29fd 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -134,13 +134,13 @@ WarpX::InitPML () if (do_pml) { pml[0].reset(new PML(boxArray(0), DistributionMap(0), &Geom(0), nullptr, - pml_ncell, pml_delta, 0, do_dive_cleaning, do_moving_window, pml_has_particles)); //pml_has_particles)); + pml_ncell, pml_delta, 0, do_dive_cleaning, do_moving_window, pml_has_particles, do_pml_in_domain)); //pml_has_particles)); for (int lev = 1; lev <= finest_level; ++lev) { pml[lev].reset(new PML(boxArray(lev), DistributionMap(lev), &Geom(lev), &Geom(lev-1), pml_ncell, pml_delta, refRatio(lev-1)[0], do_dive_cleaning, - do_moving_window, pml_has_particles)); //pml_has_particles)); + do_moving_window, pml_has_particles, do_pml_in_domain)); //pml_has_particles)); } } } @@ -322,7 +322,7 @@ WarpX::InitLevelData (int lev, Real time) void WarpX::InitLevelDataFFT (int lev, Real time) { - + Efield_fp_fft[lev][0]->setVal(0.0); Efield_fp_fft[lev][1]->setVal(0.0); Efield_fp_fft[lev][2]->setVal(0.0); diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 9d85783b0..7c00a5297 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -18,7 +18,7 @@ WarpX::ExchangeWithPmlB (int lev) Bfield_fp[lev][2].get() }, { Bfield_cp[lev][0].get(), Bfield_cp[lev][1].get(), - Bfield_cp[lev][2].get() }); + Bfield_cp[lev][2].get() }, do_pml_in_domain); } } @@ -31,7 +31,7 @@ WarpX::ExchangeWithPmlE (int lev) Efield_fp[lev][2].get() }, { Efield_cp[lev][0].get(), Efield_cp[lev][1].get(), - Efield_cp[lev][2].get() }); + Efield_cp[lev][2].get() }, do_pml_in_domain); } } @@ -40,7 +40,7 @@ WarpX::ExchangeWithPmlF (int lev) { if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeF(F_fp[lev].get(), - F_cp[lev].get()); + F_cp[lev].get(), do_pml_in_domain); } } @@ -252,7 +252,7 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) pml[lev]->ExchangeE(patch_type, { Efield_fp[lev][0].get(), Efield_fp[lev][1].get(), - Efield_fp[lev][2].get() }); + Efield_fp[lev][2].get() }, do_pml_in_domain); pml[lev]->FillBoundaryE(patch_type); } @@ -267,7 +267,7 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) pml[lev]->ExchangeE(patch_type, { Efield_cp[lev][0].get(), Efield_cp[lev][1].get(), - Efield_cp[lev][2].get() }); + Efield_cp[lev][2].get() }, do_pml_in_domain); pml[lev]->FillBoundaryE(patch_type); } @@ -294,7 +294,7 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) pml[lev]->ExchangeB(patch_type, { Bfield_fp[lev][0].get(), Bfield_fp[lev][1].get(), - Bfield_fp[lev][2].get() }); + Bfield_fp[lev][2].get() }, do_pml_in_domain); pml[lev]->FillBoundaryB(patch_type); } const auto& period = Geom(lev).periodicity(); @@ -308,7 +308,7 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) pml[lev]->ExchangeB(patch_type, { Bfield_cp[lev][0].get(), Bfield_cp[lev][1].get(), - Bfield_cp[lev][2].get() }); + Bfield_cp[lev][2].get() }, do_pml_in_domain); pml[lev]->FillBoundaryB(patch_type); } const auto& cperiod = Geom(lev-1).periodicity(); @@ -331,7 +331,7 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type) { if (do_pml && pml[lev]->ok()) { - pml[lev]->ExchangeF(patch_type, F_fp[lev].get()); + pml[lev]->ExchangeF(patch_type, F_fp[lev].get(), do_pml_in_domain); pml[lev]->FillBoundaryF(patch_type); } @@ -342,7 +342,7 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type) { if (do_pml && pml[lev]->ok()) { - pml[lev]->ExchangeF(patch_type, F_cp[lev].get()); + pml[lev]->ExchangeF(patch_type, F_cp[lev].get(), do_pml_in_domain); pml[lev]->FillBoundaryF(patch_type); } -- cgit v1.2.3 From afe646b66116f7d554b6c3874f818cb2249b2922 Mon Sep 17 00:00:00 2001 From: ablelly Date: Tue, 16 Jul 2019 19:10:56 +0200 Subject: Cleaned code. Flag `do_pml_in_domain` is working. --- Source/BoundaryConditions/PML.cpp | 18 +----------------- Source/FieldSolver/WarpXPushFieldsEM.cpp | 12 ++---------- Source/Initialization/WarpXInitData.cpp | 4 ++-- 3 files changed, 5 insertions(+), 29 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index e7c6e7bb2..3888a59d0 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -144,42 +144,35 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n if (amrex::grow(grid_box, idim, ncell).intersects(box)) { direct_faces.push_back(kv.first); - amrex::Print()<<"direct_faces"< 1, Box gaps not wide enough?\n"); } } - amrex::Print()<<"pml_type = "< Date: Thu, 25 Jul 2019 01:45:29 +0200 Subject: [WIP] Adding flags do_pml_Lo, do_pml_hi --- Source/BoundaryConditions/PML.H | 14 +++---- Source/BoundaryConditions/PML.cpp | 59 ++++++++++++++++-------------- Source/BoundaryConditions/PML_routines.F90 | 2 +- Source/FieldSolver/WarpXPushFieldsEM.cpp | 2 +- Source/Parallelization/WarpXComm.cpp | 18 ++++----- Source/WarpX.H | 4 +- Source/WarpX.cpp | 15 ++++++++ 7 files changed, 67 insertions(+), 47 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index c64eb8193..3f64db558 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -149,25 +149,25 @@ public: { return *sigba_cp; } void ExchangeB (const std::array& B_fp, - const std::array& B_cp, int do_pml_in_domain); + const std::array& B_cp, int do_pml_in_domain, int ncell); void ExchangeE (const std::array& E_fp, - const std::array& E_cp, int do_pml_in_domain); + const std::array& E_cp, int do_pml_in_domain, int ncell); void CopyJinPMLs (const std::array& j_fp, const std::array& j_cp); void CopyJinReg (const std::array& j_fp, const std::array& j_cp); void ExchangeB (PatchType patch_type, - const std::array& Bp, int do_pml_in_domain); + const std::array& Bp, int do_pml_in_domain, int ncell); void ExchangeE (PatchType patch_type, - const std::array& Ep, int do_pml_in_domain); + const std::array& Ep, int do_pml_in_domain, int ncell); void CopyJinPMLs (PatchType patch_type, const std::array& jp); void CopyJinReg (PatchType patch_type, const std::array& jp); - void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp, int do_pml_in_domain); - void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp, int do_pml_in_domain); + void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp, int do_pml_in_domain, int ncell); + void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp, int do_pml_in_domain, int ncell); void FillBoundary (); void FillBoundaryE (); @@ -205,7 +205,7 @@ private: static amrex::BoxArray MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell); - static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom, int do_pml_in_domain); + static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom, int do_pml_in_domain, int ncell); static void CopyRegInPMLs (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); static void CopyPMLsInReg (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); }; diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 51107b498..0e6c24131 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -708,53 +708,53 @@ PML::GetF_cp () void PML::ExchangeB (const std::array& B_fp, - const std::array& B_cp, int do_pml_in_domain) + const std::array& B_cp, int do_pml_in_domain, int ncell) { - ExchangeB(PatchType::fine, B_fp, do_pml_in_domain); - ExchangeB(PatchType::coarse, B_cp, do_pml_in_domain); + ExchangeB(PatchType::fine, B_fp, do_pml_in_domain, ncell); + ExchangeB(PatchType::coarse, B_cp, do_pml_in_domain, ncell); } void PML::ExchangeB (PatchType patch_type, - const std::array& Bp, int do_pml_in_domain) + const std::array& Bp, int do_pml_in_domain, int ncell) { if (patch_type == PatchType::fine && pml_B_fp[0] && Bp[0]) { - Exchange(*pml_B_fp[0], *Bp[0], *m_geom, do_pml_in_domain); - Exchange(*pml_B_fp[1], *Bp[1], *m_geom, do_pml_in_domain); - Exchange(*pml_B_fp[2], *Bp[2], *m_geom, do_pml_in_domain); + Exchange(*pml_B_fp[0], *Bp[0], *m_geom, do_pml_in_domain, ncell); + Exchange(*pml_B_fp[1], *Bp[1], *m_geom, do_pml_in_domain, ncell); + Exchange(*pml_B_fp[2], *Bp[2], *m_geom, do_pml_in_domain, ncell); } else if (patch_type == PatchType::coarse && pml_B_cp[0] && Bp[0]) { - Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom, do_pml_in_domain); - Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom, do_pml_in_domain); - Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom, do_pml_in_domain); + Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom, do_pml_in_domain, ncell); + Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom, do_pml_in_domain, ncell); + Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom, do_pml_in_domain, ncell); } } void PML::ExchangeE (const std::array& E_fp, - const std::array& E_cp, int do_pml_in_domain) + const std::array& E_cp, int do_pml_in_domain, int ncell) { - ExchangeE(PatchType::fine, E_fp, do_pml_in_domain); - ExchangeE(PatchType::coarse, E_cp, do_pml_in_domain); + ExchangeE(PatchType::fine, E_fp, do_pml_in_domain, ncell); + ExchangeE(PatchType::coarse, E_cp, do_pml_in_domain, ncell); } void PML::ExchangeE (PatchType patch_type, - const std::array& Ep, int do_pml_in_domain) + const std::array& Ep, int do_pml_in_domain, int ncell) { if (patch_type == PatchType::fine && pml_E_fp[0] && Ep[0]) { - Exchange(*pml_E_fp[0], *Ep[0], *m_geom, do_pml_in_domain); - Exchange(*pml_E_fp[1], *Ep[1], *m_geom, do_pml_in_domain); - Exchange(*pml_E_fp[2], *Ep[2], *m_geom, do_pml_in_domain); + Exchange(*pml_E_fp[0], *Ep[0], *m_geom, do_pml_in_domain, ncell); + Exchange(*pml_E_fp[1], *Ep[1], *m_geom, do_pml_in_domain, ncell); + Exchange(*pml_E_fp[2], *Ep[2], *m_geom, do_pml_in_domain, ncell); } else if (patch_type == PatchType::coarse && pml_E_cp[0] && Ep[0]) { - Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom, do_pml_in_domain); - Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom, do_pml_in_domain); - Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom, do_pml_in_domain); + Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom, do_pml_in_domain, ncell); + Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom, do_pml_in_domain, ncell); + Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom, do_pml_in_domain, ncell); } } @@ -811,29 +811,32 @@ PML::CopyJinReg (const std::array& j_fp, } void -PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp, int do_pml_in_domain) +PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp, int do_pml_in_domain, int ncell) { - ExchangeF(PatchType::fine, F_fp, do_pml_in_domain); - ExchangeF(PatchType::coarse, F_cp, do_pml_in_domain); + ExchangeF(PatchType::fine, F_fp, do_pml_in_domain, ncell); + ExchangeF(PatchType::coarse, F_cp, do_pml_in_domain, ncell); } void -PML::ExchangeF (PatchType patch_type, MultiFab* Fp, int do_pml_in_domain) +PML::ExchangeF (PatchType patch_type, MultiFab* Fp, int do_pml_in_domain, int ncell) { if (patch_type == PatchType::fine && pml_F_fp && Fp) { - Exchange(*pml_F_fp, *Fp, *m_geom, do_pml_in_domain); + Exchange(*pml_F_fp, *Fp, *m_geom, do_pml_in_domain, ncell); } else if (patch_type == PatchType::coarse && pml_F_cp && Fp) { - Exchange(*pml_F_cp, *Fp, *m_cgeom, do_pml_in_domain); + Exchange(*pml_F_cp, *Fp, *m_cgeom, do_pml_in_domain, ncell); } } void -PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom, int do_pml_in_domain) +PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom, int do_pml_in_domain, int ncell) { if (do_pml_in_domain){ const IntVect& ngr = reg.nGrowVect(); const IntVect& ngp = pml.nGrowVect(); + // amrex::Print()<< "##### PML EXCHANGE #####"<ok()) { - if (F) pml[lev]->ExchangeF(patch_type, F, do_pml_in_domain); + if (F) pml[lev]->ExchangeF(patch_type, F, do_pml_in_domain, pml_ncell); const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 7c00a5297..17d87c73f 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -18,7 +18,7 @@ WarpX::ExchangeWithPmlB (int lev) Bfield_fp[lev][2].get() }, { Bfield_cp[lev][0].get(), Bfield_cp[lev][1].get(), - Bfield_cp[lev][2].get() }, do_pml_in_domain); + Bfield_cp[lev][2].get() }, do_pml_in_domain, pml_ncell); } } @@ -31,7 +31,7 @@ WarpX::ExchangeWithPmlE (int lev) Efield_fp[lev][2].get() }, { Efield_cp[lev][0].get(), Efield_cp[lev][1].get(), - Efield_cp[lev][2].get() }, do_pml_in_domain); + Efield_cp[lev][2].get() }, do_pml_in_domain, pml_ncell); } } @@ -40,7 +40,7 @@ WarpX::ExchangeWithPmlF (int lev) { if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeF(F_fp[lev].get(), - F_cp[lev].get(), do_pml_in_domain); + F_cp[lev].get(), do_pml_in_domain, pml_ncell); } } @@ -252,7 +252,7 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) pml[lev]->ExchangeE(patch_type, { Efield_fp[lev][0].get(), Efield_fp[lev][1].get(), - Efield_fp[lev][2].get() }, do_pml_in_domain); + Efield_fp[lev][2].get() }, do_pml_in_domain, pml_ncell); pml[lev]->FillBoundaryE(patch_type); } @@ -267,7 +267,7 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) pml[lev]->ExchangeE(patch_type, { Efield_cp[lev][0].get(), Efield_cp[lev][1].get(), - Efield_cp[lev][2].get() }, do_pml_in_domain); + Efield_cp[lev][2].get() }, do_pml_in_domain, pml_ncell); pml[lev]->FillBoundaryE(patch_type); } @@ -294,7 +294,7 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) pml[lev]->ExchangeB(patch_type, { Bfield_fp[lev][0].get(), Bfield_fp[lev][1].get(), - Bfield_fp[lev][2].get() }, do_pml_in_domain); + Bfield_fp[lev][2].get() }, do_pml_in_domain, pml_ncell); pml[lev]->FillBoundaryB(patch_type); } const auto& period = Geom(lev).periodicity(); @@ -308,7 +308,7 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) pml[lev]->ExchangeB(patch_type, { Bfield_cp[lev][0].get(), Bfield_cp[lev][1].get(), - Bfield_cp[lev][2].get() }, do_pml_in_domain); + Bfield_cp[lev][2].get() }, do_pml_in_domain, pml_ncell); pml[lev]->FillBoundaryB(patch_type); } const auto& cperiod = Geom(lev-1).periodicity(); @@ -331,7 +331,7 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type) { if (do_pml && pml[lev]->ok()) { - pml[lev]->ExchangeF(patch_type, F_fp[lev].get(), do_pml_in_domain); + pml[lev]->ExchangeF(patch_type, F_fp[lev].get(), do_pml_in_domain, pml_ncell); pml[lev]->FillBoundaryF(patch_type); } @@ -342,7 +342,7 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type) { if (do_pml && pml[lev]->ok()) { - pml[lev]->ExchangeF(patch_type, F_cp[lev].get(), do_pml_in_domain); + pml[lev]->ExchangeF(patch_type, F_cp[lev].get(), do_pml_in_domain, pml_ncell); pml[lev]->FillBoundaryF(patch_type); } diff --git a/Source/WarpX.H b/Source/WarpX.H index 54cad9ae0..20db5c75d 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -489,10 +489,12 @@ private: // PML int do_pml = 1; int pml_ncell = 10; - int pml_delta = 10; + int pml_delta = 5; int pml_has_particles = 0; int do_pml_j_damping = 0; int do_pml_in_domain = 0; + amrex::IntVect do_pml_Lo; + amrex::IntVect do_pml_Hi; amrex::Vector > pml; amrex::Real moving_window_x = std::numeric_limits::max(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 08a8b7b88..46b1926a0 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -384,6 +384,21 @@ WarpX::ReadParameters () pp.query("do_pml_j_damping", do_pml_j_damping); pp.query("do_pml_in_domain", do_pml_in_domain); + Vector parse_do_pml_Lo(AMREX_SPACEDIM,1); + pp.queryarr("do_pml_Lo", parse_do_pml_Lo); + do_pml_Lo[0] = parse_do_pml_Lo[0]; + do_pml_Lo[1] = parse_do_pml_Lo[1]; +#if (AMREX_SPACEDIM == 3) + do_pml_Lo[2] = parse_do_pml_Lo[2]; +#endif + Vector parse_do_pml_Hi(AMREX_SPACEDIM,1); + pp.queryarr("do_pml_Hi", parse_do_pml_Hi); + do_pml_Hi[0] = parse_do_pml_Hi[0]; + do_pml_Hi[1] = parse_do_pml_Hi[1]; +#if (AMREX_SPACEDIM == 3) + do_pml_Hi[2] = parse_do_pml_Hi[2]; +#endif + if ( (do_pml_j_damping==1)&&(do_pml_in_domain==0) ){ amrex::Abort("J-damping can only be done when PML are inside simulation domain (do_pml_in_domain=1)"); } -- cgit v1.2.3 From 7a5adaa22b1e1830b2046aba6ba38de90c96138c Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 21 Aug 2019 13:13:23 -0700 Subject: Code cleanup + add test --- Examples/Tests/particles_in_PML/analysis_2d.py | 26 ++++++ Examples/Tests/particles_in_PML/inputs2d | 56 ++++++++++++ Source/BoundaryConditions/PML.cpp | 115 ++----------------------- Source/Evolve/WarpXEvolveEM.cpp | 18 +--- Source/FieldSolver/WarpXPushFieldsEM.cpp | 1 - Source/Parallelization/WarpXComm.cpp | 12 +-- Source/Particles/PhysicalParticleContainer.cpp | 32 +++---- Source/WarpX.H | 1 - Source/WarpX.cpp | 1 - 9 files changed, 113 insertions(+), 149 deletions(-) create mode 100755 Examples/Tests/particles_in_PML/analysis_2d.py create mode 100644 Examples/Tests/particles_in_PML/inputs2d (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Examples/Tests/particles_in_PML/analysis_2d.py b/Examples/Tests/particles_in_PML/analysis_2d.py new file mode 100755 index 000000000..5b0279a9b --- /dev/null +++ b/Examples/Tests/particles_in_PML/analysis_2d.py @@ -0,0 +1,26 @@ +""" +This script tests the absorption of particles in the PML. + +The input file inputs_2d is used: it features a positive and a +negative particle, going in opposite direction and eventually +leaving the box. This script tests that the field in the box +is close to 0 once the particles have left. With regular +PML, this test fails, since the particles leave a spurious +charge, with associated fields, behind them. +""" +import sys +import yt +import numpy as np +yt.funcs.mylog.setLevel(0) + +# Open plotfile specified in command line +filename = sys.argv[1] +ds = yt.load( filename ) + +# Check that the field is low enough +ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) +Ex_array = ad0['Ex'].to_ndarray() +Ey_array = ad0['Ey'].to_ndarray() +Ez_array = ad0['Ez'].to_ndarray() +max_Efield = max(Ex_array.max(), Ey_array.max(), Ez_array.max()) +assert max_Efield < 0.0003 \ No newline at end of file diff --git a/Examples/Tests/particles_in_PML/inputs2d b/Examples/Tests/particles_in_PML/inputs2d new file mode 100644 index 000000000..733032c53 --- /dev/null +++ b/Examples/Tests/particles_in_PML/inputs2d @@ -0,0 +1,56 @@ +max_step = 200 +amr.plot_int = 30 +amr.n_cell = 224 224 + +amr.blocking_factor = 8 +amr.max_grid_size = 1024 +amr.max_level = 0 + +warpx.dump_plotfiles = 1 + +# Geometry +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 0 0 # Is periodic? +geometry.prob_lo = -32.e-6 -32.e-6 # physical domain +geometry.prob_hi = 32.e-6 32.e-6 + +# PML +warpx.do_pml = 1 +warpx.pml_ncell = 12 +warpx.pml_delta = 6 +warpx.pml_has_particles = 0 +warpx.do_pml_in_domain = 0 +warpx.do_pml_j_damping = 0 + + +# Algorithms +algo.current_deposition = esirkepov +algo.charge_deposition = standard +algo.field_gathering = vectorized +algo.particle_pusher = vay +algo.maxwell_fdtd_solver = ckc +warpx.cfl = 1.0 +warpx.use_filter = 1 + +# Particle species +particles.nspecies = 2 +particles.species_names = electron proton + +electron.charge = -q_e +electron.mass = m_e +electron.injection_style = "singleparticle" +electron.single_particle_pos = 0. 0. 0. +electron.single_particle_vel = 2. 0. 0. +electron.single_particle_weight = 1. + +proton.charge = q_e +proton.mass = m_p # Very heavy ; should not move +proton.injection_style = "singleparticle" +proton.single_particle_pos = 0. 0. 0. +proton.single_particle_vel = -2. 0. 0. +proton.single_particle_weight = 1. + +# Particle shape factor in each direction +interpolation.nox = 3 +interpolation.noy = 3 +interpolation.noz = 3 diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 0b73859b4..15f18201b 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -24,22 +24,19 @@ namespace int slo = sigma.m_lo; int shi = sigma.m_hi; int sslo = sigma_star.m_lo; - Real x = 10.0; - Real theta = 10.0; - Real coeff_damp = 1.; //std::sin(theta*MathConst::pi/180.); for (int i = olo; i <= ohi+1; ++i) { Real offset = static_cast(glo-i); sigma[i-slo] = fac*(offset*offset); - sigma_cum[i-slo] = coeff_damp*(fac*(offset*offset*offset)/3.)/(PhysConst::c*x/std::sqrt(1+x*x)); + sigma_cum[i-slo] = (fac*(offset*offset*offset)/3.)/PhysConst::c; } for (int i = olo; i <= ohi; ++i) { Real offset = static_cast(glo-i) - 0.5; sigma_star[i-sslo] = fac*(offset*offset); - sigma_star_cum[i-sslo] = coeff_damp*(fac*(offset*offset*offset)/3.)/(PhysConst::c*x/std::sqrt(1+x*x)); + sigma_star_cum[i-sslo] = (fac*(offset*offset*offset)/3.)/PhysConst::c; } } @@ -52,20 +49,17 @@ namespace int slo = sigma.m_lo; int shi = sigma.m_hi; int sslo = sigma_star.m_lo; - Real x = 10.0; - Real theta = 10.0; - Real coeff_damp = 1.;//std::sin(theta*MathConst::pi/180.); for (int i = olo; i <= ohi+1; ++i) { Real offset = static_cast(i-ghi-1); sigma[i-slo] = fac*(offset*offset); - sigma_cum[i-slo] = coeff_damp*(fac*(offset*offset*offset)/3.)/(PhysConst::c*x/std::sqrt(1+x*x)); + sigma_cum[i-slo] = coeff_damp*(fac*(offset*offset*offset)/3.)/PhysConst::c; } for (int i = olo; i <= ohi; ++i) { Real offset = static_cast(i-ghi) - 0.5; sigma_star[i-sslo] = fac*(offset*offset); - sigma_star_cum[i-sslo] = coeff_damp*(fac*(offset*offset*offset)/3.)/(PhysConst::c*x/std::sqrt(1+x*x)); + sigma_star_cum[i-sslo] = coeff_damp*(fac*(offset*offset*offset)/3.)/PhysConst::c; } } @@ -119,8 +113,6 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n sigma_star_cum_fac[idim].m_hi = hi[idim]; } - amrex::Print()< 0){pml_type[2]=idim+'0';} //std::to_string(idim) + if (idim > 0){pml_type[2]=idim+'0';} } const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell); Box hioverlap = hibox & box; if (hioverlap.ok()) { - amrex::Print()<<"["< 1, Box gaps not wide enough?\n"); } } - amrex::Print()< 0) // Copy from pml to the ghost cells of regular data { - MultiFab::Copy(tmpregmf, reg, 0, 0, 1, IntVect(0)); //ngr); + MultiFab::Copy(tmpregmf, reg, 0, 0, 1, IntVect(0)); MultiFab::Copy(totpmlmf, pml, 0, 0, ncp, ngp); tmpregmf.setVal(0.0, 1, ncp-1, 0); - // ce qui foncitonne a peu pres pour le moment tmpregmf.ParallelCopy(totpmlmf,0, 0, ncp, IntVect(0), IntVect(0), period); totpmlmf.ParallelCopy(tmpregmf,0, 0, ncp, IntVect(0), ngp, period); - MultiFab::Copy(pml, totpmlmf, 0, 0, ncp, ngp); //test -// #ifdef _OPENMP -// #pragma omp parallel -// #endif - // // amrex::Print()<<"####### EXCHANGE INFORMATIONS GHOST CELLS PML #######"<>> EXCHANGE BOXES <<<<<" <>> EXCHANGE BOXES <<<<<" < 0 && step < load_balance_max_step && (step+1) % load_balance_int == 0) + if (step > 0 && (step+1) % load_balance_int == 0) { LoadBalance(); // Reset the costs to 0 @@ -329,22 +329,6 @@ WarpX::OneStep_nosub (Real cur_time) EvolveE(dt[0]); // We now have E^{n+1} FillBoundaryE(); EvolveF(0.5*dt[0], DtType::SecondHalf); - - - // if (do_pml && pml_has_particles){ - // for (int lev = 0; lev <= finest_level; ++lev) - // { - // if (pml[lev]->ok()){ - // pml[lev]->CopyJinReg({ current_fp[lev][0].get(), - // current_fp[lev][1].get(), - // current_fp[lev][2].get() }, - // { current_cp[lev][0].get(), - // current_cp[lev][1].get(), - // current_cp[lev][2].get() }); - // } - // } - // } - EvolveB(0.5*dt[0]); // We now have B^{n+1} if (do_pml) { DampPML(); diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 18c98b8f8..011d9cfcc 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -241,7 +241,6 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); - // const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp(); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 06cbf157a..c539cd89a 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -256,7 +256,7 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeE(patch_type, - { Efield_fp[lev][0].get(), + { Efield_fp[lev][0].get(), Efield_fp[lev][1].get(), Efield_fp[lev][2].get() }, do_pml_in_domain, pml_ncell, @@ -273,7 +273,7 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeE(patch_type, - { Efield_cp[lev][0].get(), + { Efield_cp[lev][0].get(), Efield_cp[lev][1].get(), Efield_cp[lev][2].get() }, do_pml_in_domain, pml_ncell, @@ -302,7 +302,7 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeB(patch_type, - { Bfield_fp[lev][0].get(), + { Bfield_fp[lev][0].get(), Bfield_fp[lev][1].get(), Bfield_fp[lev][2].get() }, do_pml_in_domain, pml_ncell, @@ -318,9 +318,9 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeB(patch_type, - { Bfield_cp[lev][0].get(), - Bfield_cp[lev][1].get(), - Bfield_cp[lev][2].get() }, + { Bfield_cp[lev][0].get(), + Bfield_cp[lev][1].get(), + Bfield_cp[lev][2].get() }, do_pml_in_domain, pml_ncell, do_pml_Lo, do_pml_Hi); pml[lev]->FillBoundaryB(patch_type); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 20819a78c..d10390204 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -33,7 +33,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); pp.query("do_continuous_injection", do_continuous_injection); - // Whether to plot back-transformed (lab-frame) diagnostics + // Whether to plot back-transformed (lab-frame) diagnostics // for this species. pp.query("do_boosted_frame_diags", do_boosted_frame_diags); @@ -42,7 +42,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); if (not do_user_plot_vars){ // By default, all particle variables are dumped to plotfiles, - // including {x,y,z,ux,uy,uz}old variables when running in a + // including {x,y,z,ux,uy,uz}old variables when running in a // boosted frame if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ plot_flags.resize(PIdx::nattribs + 6, 1); @@ -59,9 +59,9 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // If not none, set plot_flags values to 1 for elements in plot_vars. if (plot_vars[0] != "none"){ for (const auto& var : plot_vars){ - // Return error if var not in PIdx. - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - ParticleStringNames::to_index.count(var), + // Return error if var not in PIdx. + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ParticleStringNames::to_index.count(var), "plot_vars argument not in ParticleStringNames"); plot_flags[ParticleStringNames::to_index.at(var)] = 1; } @@ -130,7 +130,7 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real void PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real x_rms, Real y_rms, Real z_rms, - Real q_tot, long npart, + Real q_tot, long npart, int do_symmetrize) { const Geometry& geom = m_gdb->Geom(0); @@ -203,7 +203,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, particle_tile.push_back_real(particle_comps["xold"], x); particle_tile.push_back_real(particle_comps["yold"], y); particle_tile.push_back_real(particle_comps["zold"], z); - + particle_tile.push_back_real(particle_comps["uxold"], u[0]); particle_tile.push_back_real(particle_comps["uyold"], u[1]); particle_tile.push_back_real(particle_comps["uzold"], u[2]); @@ -931,11 +931,11 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_NS("PICSAR::FieldGather", blp_pxr_fg); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition); - + const std::array& dx = WarpX::CellSize(lev); const std::array& cdx = WarpX::CellSize(std::max(lev-1,0)); - // Get instances of NCI Godfrey filters + // Get instances of NCI Godfrey filters const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; @@ -949,7 +949,7 @@ PhysicalParticleContainer::Evolve (int lev, bool has_buffer = cEx || cjx; #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -1227,13 +1227,13 @@ PhysicalParticleContainer::Evolve (int lev, eyeli = filtered_Ey.elixir(); nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], filtered_Ey.box()); ceyfab = &filtered_Ey; - + // Filter Bx filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); bxeli = filtered_Bx.elixir(); nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Bx, (*cBx)[pti], filtered_Bx.box()); cbxfab = &filtered_Bx; - + // Filter Bz filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); bzeli = filtered_Bz.elixir(); @@ -1453,7 +1453,7 @@ PhysicalParticleContainer::SplitParticles(int lev) } // Add local arrays psplit_x etc. to the temporary // particle container pctmp_split. Split particles - // are tagged with p.id()=NoSplitParticleID so that + // are tagged with p.id()=NoSplitParticleID so that // they are not re-split when entering a higher level // AddNParticles calls Redistribute, so that particles // in pctmp_split are in the proper grids and tiles @@ -1550,7 +1550,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int thread_num = omp_get_thread_num(); #else int thread_num = 0; -#endif +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -1694,7 +1694,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real slice_box.setHi(direction, z_max); diagnostic_particles.resize(finestLevel()+1); - + for (int lev = 0; lev < nlevs; ++lev) { const Real* dx = Geom(lev).CellSize(); @@ -1776,7 +1776,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real Real uzp = uz_old_p *weight_old + uz_new_p *weight_new; diagnostic_particles[lev][index].GetRealData(DiagIdx::w).push_back(wp[i]); - + diagnostic_particles[lev][index].GetRealData(DiagIdx::x).push_back(xp); diagnostic_particles[lev][index].GetRealData(DiagIdx::y).push_back(yp); diagnostic_particles[lev][index].GetRealData(DiagIdx::z).push_back(zp); diff --git a/Source/WarpX.H b/Source/WarpX.H index b4dfe9ae4..3e8d45f49 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -524,7 +524,6 @@ private: amrex::Real const_dt = 0.5e-11; int load_balance_int = -1; - int load_balance_max_step = std::numeric_limits::max(); amrex::Vector > costs; int load_balance_with_sfc = 0; amrex::Real load_balance_knapsack_factor = 1.24; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 6abf27216..305fedbb4 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -498,7 +498,6 @@ WarpX::ReadParameters () } pp.query("load_balance_int", load_balance_int); - pp.query("load_balance_max_step", load_balance_max_step); pp.query("load_balance_with_sfc", load_balance_with_sfc); pp.query("load_balance_knapsack_factor", load_balance_knapsack_factor); -- cgit v1.2.3 From 95ca1710f97fd47a9098f13535e5f5aa2e4e339c Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 27 Aug 2019 14:19:16 -0700 Subject: Call C++ functions --- Source/BoundaryConditions/PML_current.H | 118 ++---- Source/BoundaryConditions/PML_routines.F90 | 563 +++++++---------------------- Source/FieldSolver/WarpXPushFieldsEM.cpp | 52 ++- Source/FortranInterface/WarpX_f.H | 22 -- 4 files changed, 208 insertions(+), 547 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML_current.H b/Source/BoundaryConditions/PML_current.H index 8bc2060bf..413b068ca 100644 --- a/Source/BoundaryConditions/PML_current.H +++ b/Source/BoundaryConditions/PML_current.H @@ -8,122 +8,74 @@ using namespace amrex; AMREX_GPU_HOST_DEVICE AMREX_INLINE void warpx_push_ex_pml_current_nodal (int j, int k, int l, Array4 const& Ex, Array4 const& jx, - Array4 const& sigjy, - Array4 const& sigjz, - Real mudt, IntVect pml_type) + Real const* const& sigjy, + Real const* const& sigjz, + Real mu_c2_dt) { #if (AMREX_SPACEDIM == 3) Real alpha_xy, alpha_xz; - if (pml_type[0] == 0){ - if (pml_type[1] == 0){ - alpha_xy = 0.5; - alpha_xz = 0.5; - } - if (pml_type[1] == 1){ - alpha_xy = 1.; - alpha_xz = 0.; - } - if (pml_type[1] == 2){ - alpha_xy = 0.; - alpha_xz = 1.; - } + if (sigjy[k]+sigjz[l] == 0){ + alpha_xy = 0.5; + alpha_xz = 0.5; } else { - if ((sigjy[k] == 0.) && (sigjz[l] == 0.)){ - alpha_xy = 0.5; - alpha_xz = 0.5; - } - else { - alpha_xy = sigjy[k]/(sigjy[k]+sigjz[l]); - alpha_xz = sigjz[l]/(sigjy[k]+sigjz[l]); - } + alpha_xy = sigjy[k]/(sigjy[k]+sigjz[l]); + alpha_xz = sigjz[l]/(sigjy[k]+sigjz[l]); } - Ex(j,k,l,0) = Ex(j,k,l,0) - mudt * alpha_xy * jx(j,k,l); - Ex(j,k,l,1) = Ex(j,k,l,1) - mudt * alpha_xz * jx(j,k,l); + Ex(j,k,l,0) = Ex(j,k,l,0) - mu_c2_dt * alpha_xy * jx(j,k,l); + Ex(j,k,l,1) = Ex(j,k,l,1) - mu_c2_dt * alpha_xz * jx(j,k,l); #else - Ex(j,k,0) = Ex(j,k,0) - mudt * jx(j,k); + Ex(j,k,0) = Ex(j,k,0) - mu_c2_dt * jx(j,k); #endif } AMREX_GPU_HOST_DEVICE AMREX_INLINE void warpx_push_ey_pml_current_nodal (int j, int k, int l, Array4 const& Ey, Array4 const& jy, - Array4 const& sigjx, - Array4 const& sigjz, - Real mudt, int pml_type) + Real const* const& sigjx, + Real const* const& sigjz, + Real mu_c2_dt) { #if (AMREX_SPACEDIM == 3) - Real alpha_yx, alpha yz; - if (pml_type[0] == 0){ - if (pml_type[1] == 0){ - alpha_yx = 1.; - alpha_yz = 0.; - } - if (pml_type[1] == 1){ - alpha_yx = 0.5; - alpha_yz = 0.5; - } - if (pml_type[1] == 2){ - alpha_yx = 0.; - alpha_yz = 1.; - } +Real alpha_yx, alpha_yz; + if (sigjx[j]+sigjz[l] == 0){ + alpha_yx = 0.5; + alpha_yz = 0.5; } else { - if ((sigjy[k] == 0.) && (sigjz[l] == 0.)){ - alpha_yx = 0.5; - alpha_yz = 0.5; - } - else { - alpha_yx = sigjx[j]/(sigjx[j]+sigjz[l]); - alpha_yz = sigjz[l]/(sigjx[j]+sigjz[l]); - } + alpha_yx = sigjx[j]/(sigjx[j]+sigjz[l]); + alpha_yz = sigjz[l]/(sigjx[j]+sigjz[l]); } - Ey(j,k,l,0) = Ey(j,k,l,0) - mudt * alpha_yx * jy(j,k,l); - Ey(j,k,l,1) = Ey(j,k,l,1) - mudt * alpha_yz * jy(j,k,l); + Ey(j,k,l,0) = Ey(j,k,l,0) - mu_c2_dt * alpha_yx * jy(j,k,l); + Ey(j,k,l,1) = Ey(j,k,l,1) - mu_c2_dt * alpha_yz * jy(j,k,l); #else - Ey(j,k,0) = Ey(j,k,0) - 0.5 * mudt * jy(j,k); - Ey(j,k,1) = Ey(j,k,0) - 0.5 * mudt * jy(j,k); + Ey(j,k,0) = Ey(j,k,0) - 0.5 * mu_c2_dt * jy(j,k); + Ey(j,k,1) = Ey(j,k,0) - 0.5 * mu_c2_dt * jy(j,k); #endif } AMREX_GPU_HOST_DEVICE AMREX_INLINE void warpx_push_ez_pml_current_nodal (int j, int k, int l, Array4 const& Ez, Array4 const& jz, - Array4 const& sigjx, - Array4 const& sigjy, - Real mudt, int pml_type) + Real const* const& sigjx, + Real const* const& sigjy, + Real mu_c2_dt) { #if (AMREX_SPACEDIM == 3) Real alpha_zx, alpha_zy; - if (pml_type[0] == 0){ - if (pml_type[1] == 0){ - alpha_zx = 1.; - alpha_zy = 0.; - } - if (pml_type[1] == 1){ - alpha_zx = 0.; - alpha_zy = 1.; - } - if (pml_type[1] == 2){ - alpha_zx = 0.5; - alpha_zy = 0.5; - } + if (sigjx[j]+sigjy[k]==0){ + alpha_zx = 0.5; + alpha_zy = 0.5; } else { - if ((sigjy[k] == 0.) && (sigjz[l] == 0.)){ - alpha_zx = 0.5; - alpha_zy = 0.5; - } - else { - alpha_zx = sigjx[j]/(sigjx[j]+sigjy[k]); - alpha_zy = sigjy[k]/(sigjx[j]+sigjy[k]); - } + alpha_zx = sigjx[j]/(sigjx[j]+sigjy[k]); + alpha_zy = sigjy[k]/(sigjx[j]+sigjy[k]); } - Ez(j,k,l,0) = Ez(j,k,l,0) - mudt * alpha_zx * jz(j,k,l); - Ez(j,k,l,1) = Ez(j,k,l,1) - mudt * alpha_zy * jz(j,k,l); + Ez(j,k,l,0) = Ez(j,k,l,0) - mu_c2_dt * alpha_zx * jz(j,k,l); + Ez(j,k,l,1) = Ez(j,k,l,1) - mu_c2_dt * alpha_zy * jz(j,k,l); #else - Ez(j,k,0) = Ez(j,k,0) - mudt * jz(j,k); + Ez(j,k,0) = Ez(j,k,0) - mu_c2_dt * jz(j,k); #endif } diff --git a/Source/BoundaryConditions/PML_routines.F90 b/Source/BoundaryConditions/PML_routines.F90 index 8f617d86e..69c762dd1 100644 --- a/Source/BoundaryConditions/PML_routines.F90 +++ b/Source/BoundaryConditions/PML_routines.F90 @@ -255,295 +255,53 @@ contains & Bx, Bxlo, Bxhi, & & By, Bylo, Byhi, & & Bz, Bzlo, Bzhi, & - & jx, jxlo, jxhi, & - & jy, jylo, jyhi, & - & jz, jzlo, jzhi, & - & flag, & - & pml_type, & - & sigjx, sjxlo, sjxhi, & - & sigjy, sjylo, sjyhi, & - & sigjz, sjzlo, sjzhi, & - & sigsjx, ssjxlo, ssjxhi, & - & sigsjy, ssjylo, ssjyhi, & - & sigsjz, ssjzlo, ssjzhi, & - & mudt, & & dtsdx, dtsdy, dtsdz) & bind(c,name='warpx_push_pml_evec_3d') integer, intent(in) :: xlo(3), xhi(3), ylo(3), yhi(3), zlo(3), zhi(3), & Exlo(3), Exhi(3), Eylo(3), Eyhi(3), Ezlo(3), Ezhi(3), & - Bxlo(3), Bxhi(3), Bylo(3), Byhi(3), Bzlo(3), Bzhi(3), & - jxlo(3), jxhi(3), jylo(3), jyhi(3), jzlo(3), jzhi(3), & - sjxlo, sjxhi, sjylo, sjyhi, sjzlo, sjzhi, & - ssjxlo, ssjxhi, ssjylo, ssjyhi, ssjzlo, ssjzhi, & - flag, pml_type(5) - real(amrex_real), intent(in ) :: mudt, dtsdx, dtsdy, dtsdz + Bxlo(3), Bxhi(3), Bylo(3), Byhi(3), Bzlo(3), Bzhi(3) + real(amrex_real), intent(in) :: dtsdx, dtsdy, dtsdz real(amrex_real), intent(inout) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),Exlo(3):Exhi(3),2) real(amrex_real), intent(inout) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),Eylo(3):Eyhi(3),2) real(amrex_real), intent(inout) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),Ezlo(3):Ezhi(3),2) real(amrex_real), intent(in ) :: Bx (Bxlo(1):Bxhi(1),Bxlo(2):Bxhi(2),Bxlo(3):Bxhi(3),2) real(amrex_real), intent(in ) :: By (Bylo(1):Byhi(1),Bylo(2):Byhi(2),Bylo(3):Byhi(3),2) real(amrex_real), intent(in ) :: Bz (Bzlo(1):Bzhi(1),Bzlo(2):Bzhi(2),Bzlo(3):Bzhi(3),2) - real(amrex_real), intent(in ) :: jx (jxlo(1):jxhi(1),jxlo(2):jxhi(2),jxlo(3):jxhi(3)) - real(amrex_real), intent(in ) :: jy (jylo(1):jyhi(1),jylo(2):jyhi(2),jylo(3):jyhi(3)) - real(amrex_real), intent(in ) :: jz (jzlo(1):jzhi(1),jzlo(2):jzhi(2),jzlo(3):jzhi(3)) - real(amrex_real), intent(in ) :: sigjx (sjxlo:sjxhi) - real(amrex_real), intent(in ) :: sigjy (sjylo:sjyhi) - real(amrex_real), intent(in ) :: sigjz (sjzlo:sjzhi) - real(amrex_real), intent(in ) :: sigsjx (ssjxlo:ssjxhi) - real(amrex_real), intent(in ) :: sigsjy (ssjylo:ssjyhi) - real(amrex_real), intent(in ) :: sigsjz (ssjzlo:ssjzhi) integer :: i, j, k - real(amrex_real) :: alpha_xy, alpha_xz, alpha_yx, alpha_yz, alpha_zx, alpha_zy - - if (flag == 0) then - do k = xlo(3), xhi(3) - do j = xlo(2), xhi(2) - do i = xlo(1), xhi(1) - 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)) - Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & - & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) - end do - end do - end do - - do k = ylo(3), yhi(3) - do j = ylo(2), yhi(2) - do i = ylo(1), yhi(1) - 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)) - Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & - & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) - end do - end do - end do - - do k = zlo(3), zhi(3) - do j = zlo(2), zhi(2) - do i = zlo(1), zhi(1) - 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)) - Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & - & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) - end do - end do - end do - - else !flag = 1 - !!!!! DIRECT FACE - if (pml_type(1)==0) then - if (pml_type(2)==0) then - do k = xlo(3), xhi(3) - do j = xlo(2), xhi(2) - do i = xlo(1), xhi(1) - ! compute current coefficients alpha - alpha_xy = 0.5 - alpha_xz = 0.5 - 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) - Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & - & -By(i,j ,k-1,1)-By(i,j ,k-1,2))& - & -mudt*alpha_xz*jx(i,j,k) - end do - end do - end do - - do k = ylo(3), yhi(3) - do j = ylo(2), yhi(2) - do i = ylo(1), yhi(1) - alpha_yx = 1. - alpha_yz = 0. - 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_yz*jy(i,j,k) - Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & - & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2))& - & -mudt*alpha_yx*jy(i,j,k) - end do - end do - end do - - do k = zlo(3), zhi(3) - do j = zlo(2), zhi(2) - do i = zlo(1), zhi(1) - alpha_zx = 1. - alpha_zy = 0. - 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) - Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & - & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2))& - & -mudt*alpha_zy*jz(i,j,k) - end do - end do - end do - - else if (pml_type(2)==1) then - do k = xlo(3), xhi(3) - do j = xlo(2), xhi(2) - do i = xlo(1), xhi(1) - ! compute current coefficients alpha - alpha_xy = 1. - alpha_xz = 0. - 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) - Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & - & -By(i,j ,k-1,1)-By(i,j ,k-1,2))& - & -mudt*alpha_xz*jx(i,j,k) - end do - end do - end do - - do k = ylo(3), yhi(3) - do j = ylo(2), yhi(2) - do i = ylo(1), yhi(1) - alpha_yx = 0.5 - alpha_yz = 0.5 - 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_yz*jy(i,j,k) - Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & - & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2))& - & -mudt*alpha_yx*jy(i,j,k) - end do - end do - end do - - do k = zlo(3), zhi(3) - do j = zlo(2), zhi(2) - do i = zlo(1), zhi(1) - alpha_zx = 0. - alpha_zy = 1. - 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) - Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & - & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2))& - & -mudt*alpha_zy*jz(i,j,k) - end do - end do - end do - - else !(pml_type(2)==2) - do k = xlo(3), xhi(3) - do j = xlo(2), xhi(2) - do i = xlo(1), xhi(1) - ! compute current coefficients alpha - alpha_xy = 0. - alpha_xz = 1. - 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) - Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & - & -By(i,j ,k-1,1)-By(i,j ,k-1,2))& - & -mudt*alpha_xz*jx(i,j,k) - end do - end do - end do - - do k = ylo(3), yhi(3) - do j = ylo(2), yhi(2) - do i = ylo(1), yhi(1) - alpha_yx = 0. - alpha_yz = 1. - 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_yz*jy(i,j,k) - Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & - & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2))& - & -mudt*alpha_yx*jy(i,j,k) - end do - end do - end do - do k = zlo(3), zhi(3) - do j = zlo(2), zhi(2) - do i = zlo(1), zhi(1) - alpha_zx = 0.5 - alpha_zy = 0.5 - 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) - Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & - & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2))& - & -mudt*alpha_zy*jz(i,j,k) - end do - end do - end do - end if !DIRECT FACE - - !!!!! SIDE EDGE OR CORNER - else - do k = xlo(3), xhi(3) - do j = xlo(2), xhi(2) - do i = xlo(1), xhi(1) - ! compute current coefficients alpha - 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 - - 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) - Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & - & -By(i,j ,k-1,1)-By(i,j ,k-1,2))& - & -mudt*alpha_xz*jx(i,j,k) - end do - end do - end do - - do k = ylo(3), yhi(3) - do j = ylo(2), yhi(2) - do i = ylo(1), yhi(1) - 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 - - 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_yz*jy(i,j,k) - Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & - & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2))& - & -mudt*alpha_yx*jy(i,j,k) - end do - end do - end do - - do k = zlo(3), zhi(3) - do j = zlo(2), zhi(2) - do i = zlo(1), zhi(1) - 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 - - 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) - Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & - & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2))& - & -mudt*alpha_zy*jz(i,j,k) - end do - end do - end do + do k = xlo(3), xhi(3) + do j = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + 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)) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) + end do + end do + end do - end if !PML_TYPE + do k = ylo(3), yhi(3) + do j = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + 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)) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) + end do + end do + end do - end if ! FLAG + do k = zlo(3), zhi(3) + do j = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + 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)) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) + end do + end do + end do end subroutine warpx_push_pml_evec_3d @@ -673,83 +431,43 @@ contains & Bx, Bxlo, Bxhi, & & By, Bylo, Byhi, & & Bz, Bzlo, Bzhi, & - & jx, jxlo, jxhi, & - & jy, jylo, jyhi, & - & jz, jzlo, jzhi, & - & flag, & - & mudt, dtsdx, dtsdy, dtsdz) & + & dtsdx, dtsdy, dtsdz) & bind(c,name='warpx_push_pml_evec_2d') integer, intent(in) :: xlo(2), xhi(2), ylo(2), yhi(2), zlo(2), zhi(2), & Exlo(2), Exhi(2), Eylo(2), Eyhi(2), Ezlo(2), Ezhi(2), & - Bxlo(2), Bxhi(2), Bylo(2), Byhi(2), Bzlo(2), Bzhi(2), & - jxlo(2), jxhi(2), jylo(2), jyhi(2), jzlo(2), jzhi(2), & - flag - real(amrex_real), intent(in) :: mudt, dtsdx, dtsdy, dtsdz + Bxlo(2), Bxhi(2), Bylo(2), Byhi(2), Bzlo(2), Bzhi(2) + real(amrex_real), intent(in) :: dtsdx, dtsdy, dtsdz real(amrex_real), intent(inout) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),2) real(amrex_real), intent(inout) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),2) real(amrex_real), intent(inout) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),2) real(amrex_real), intent(in ) :: Bx (Bxlo(1):Bxhi(1),Bxlo(2):Bxhi(2),2) real(amrex_real), intent(in ) :: By (Bylo(1):Byhi(1),Bylo(2):Byhi(2),2) real(amrex_real), intent(in ) :: Bz (Bzlo(1):Bzhi(1),Bzlo(2):Bzhi(2),2) - real(amrex_real), intent(in ) :: jx (jxlo(1):jxhi(1),jxlo(2):jxhi(2)) !jx (jxlo(1):jxhi(1),jxlo(2):jxhi(2),1) - real(amrex_real), intent(in ) :: jy (jylo(1):jyhi(1),jylo(2):jyhi(2)) !jy (jylo(1):jyhi(1),jylo(2):jyhi(2),1) - real(amrex_real), intent(in ) :: jz (jzlo(1):jzhi(1),jzlo(2):jzhi(2)) !jz (jzlo(1):jzhi(1),jzlo(2):jzhi(2),1) integer :: i, k - if (flag == 1) then - do k = xlo(2), xhi(2) - do i = xlo(1), xhi(1) - Ex(i,k,2) = Ex(i,k,2) - dtsdz*(By(i,k ,1)+By(i,k ,2) & - & -By(i,k-1,1)-By(i,k-1,2))& - & - mudt * jx(i,k) - end do - end do - - do k = ylo(2), yhi(2) - do i = ylo(1), yhi(1) - Ey(i,k,1) = Ey(i,k,1) + dtsdz*(Bx(i ,k ,1)+Bx(i ,k ,2) & - & -Bx(i ,k-1,1)-Bx(i ,k-1,2))& - & - mudt * 0.5 * jy(i,k) - Ey(i,k,2) = Ey(i,k,2) - dtsdx*(Bz(i ,k ,1)+Bz(i ,k ,2) & - & -Bz(i-1,k ,1)-Bz(i-1,k ,2))& - & - mudt * 0.5 * jy(i,k) - end do - end do - - do k = zlo(2), zhi(2) - do i = zlo(1), zhi(1) - Ez(i,k,1) = Ez(i,k,1) + dtsdx*(By(i ,k,1)+By(i ,k,2) & - & -By(i-1,k,1)-By(i-1,k,2))& - & - mudt * jz(i,k) - end do - end do - - else !no particles in PML - do k = xlo(2), xhi(2) - do i = xlo(1), xhi(1) - Ex(i,k,2) = Ex(i,k,2) - dtsdz*(By(i,k ,1)+By(i,k ,2) & - & -By(i,k-1,1)-By(i,k-1,2)) - end do - end do - - do k = ylo(2), yhi(2) - do i = ylo(1), yhi(1) - Ey(i,k,1) = Ey(i,k,1) + dtsdz*(Bx(i ,k ,1)+Bx(i ,k ,2) & - & -Bx(i ,k-1,1)-Bx(i ,k-1,2)) - Ey(i,k,2) = Ey(i,k,2) - dtsdx*(Bz(i ,k ,1)+Bz(i ,k ,2) & - & -Bz(i-1,k ,1)-Bz(i-1,k ,2)) - end do - end do + do k = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + Ex(i,k,2) = Ex(i,k,2) - dtsdz*(By(i,k ,1)+By(i,k ,2) & + & -By(i,k-1,1)-By(i,k-1,2)) + end do + end do - do k = zlo(2), zhi(2) - do i = zlo(1), zhi(1) - Ez(i,k,1) = Ez(i,k,1) + dtsdx*(By(i ,k,1)+By(i ,k,2) & - & -By(i-1,k,1)-By(i-1,k,2)) - end do - end do + do k = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + Ey(i,k,1) = Ey(i,k,1) + dtsdz*(Bx(i ,k ,1)+Bx(i ,k ,2) & + & -Bx(i ,k-1,1)-Bx(i ,k-1,2)) + Ey(i,k,2) = Ey(i,k,2) - dtsdx*(Bz(i ,k ,1)+Bz(i ,k ,2) & + & -Bz(i-1,k ,1)-Bz(i-1,k ,2)) + end do + end do - end if !flag + do k = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + Ez(i,k,1) = Ez(i,k,1) + dtsdx*(By(i ,k,1)+By(i ,k,2) & + & -By(i-1,k,1)-By(i-1,k,2)) + end do + end do end subroutine warpx_push_pml_evec_2d @@ -1109,7 +827,6 @@ contains do k = texlo(2), texhi(2) do i = texlo(1), texhi(1) - ! ex(i,k,1) = ex(i,k,1) * sigez(k) !No damping for Exy ex(i,k,2) = ex(i,k,2) * sigez(k) ex(i,k,3) = ex(i,k,3) * sigcx(i) end do @@ -1150,91 +867,6 @@ contains end subroutine warpx_damp_pml_2d - subroutine warpx_dampJ_pml_2d (tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & - & jx, jxlo, jxhi, jy, jylo, jyhi, jz, jzlo, jzhi, & - & sigjx, sjxlo, sjxhi, sigjz, sjzlo, sjzhi, & - & sigsjx, ssjxlo, ssjxhi, sigsjz, ssjzlo, ssjzhi) & - bind(c,name='warpx_dampJ_pml_2d') - integer, dimension(2), intent(in) :: tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & - jxlo, jxhi, jylo, jyhi, jzlo, jzhi - integer, intent(in), value :: sjxlo, sjxhi, sjzlo, sjzhi, ssjxlo, ssjxhi, ssjzlo, ssjzhi - real(amrex_real), intent(inout) :: jx(jxlo(1):jxhi(1),jxlo(2):jxhi(2)) - real(amrex_real), intent(inout) :: jy(jylo(1):jyhi(1),jylo(2):jyhi(2)) - real(amrex_real), intent(inout) :: jz(jzlo(1):jzhi(1),jzlo(2):jzhi(2)) - real(amrex_real), intent(in) :: sigjx(sjxlo:sjxhi) - real(amrex_real), intent(in) :: sigjz(sjzlo:sjzhi) - real(amrex_real), intent(in) :: sigsjx(ssjxlo:ssjxhi) - real(amrex_real), intent(in) :: sigsjz(ssjzlo:ssjzhi) - - integer :: i,k - - do k = tjxlo(2), tjxhi(2) - do i = tjxlo(1), tjxhi(1) - jx(i,k) = jx(i,k) * sigsjx(i) * sigjz(k) - end do - end do - - ! do k = tjylo(2), tjyhi(2) - ! do i = tjylo(1), tjyhi(1) - ! jy(i,k) = jy(i,k) - ! end do - ! end do - - do k = tjzlo(2), tjzhi(2) - do i = tjzlo(1), tjzhi(1) - jz(i,k) = jz(i,k) * sigjx(i) * sigsjz(k) - end do - end do - - - end subroutine warpx_dampJ_pml_2d - - subroutine warpx_dampJ_pml_3d (tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & - & jx, jxlo, jxhi, jy, jylo, jyhi, jz, jzlo, jzhi, & - & sigjx, sjxlo, sjxhi, sigjy, sjylo, sjyhi, sigjz, sjzlo, sjzhi, & - & sigsjx, ssjxlo, ssjxhi, sigsjy, ssjylo, ssjyhi, sigsjz, ssjzlo, ssjzhi) & - bind(c,name='warpx_dampJ_pml_3d') - integer, dimension(3), intent(in) :: tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & - jxlo, jxhi, jylo, jyhi, jzlo, jzhi - integer, intent(in), value :: sjxlo, sjxhi, sjylo, sjyhi, sjzlo, sjzhi, ssjxlo, ssjxhi, ssjylo, ssjyhi, ssjzlo, ssjzhi - real(amrex_real), intent(inout) :: jx(jxlo(1):jxhi(1),jxlo(2):jxhi(2),jxlo(3):jxhi(3)) - real(amrex_real), intent(inout) :: jy(jylo(1):jyhi(1),jylo(2):jyhi(2),jylo(3):jyhi(3)) - real(amrex_real), intent(inout) :: jz(jzlo(1):jzhi(1),jzlo(2):jzhi(2),jzlo(3):jzhi(3)) - real(amrex_real), intent(in) :: sigjx(sjxlo:sjxhi) - real(amrex_real), intent(in) :: sigjy(sjylo:sjyhi) - real(amrex_real), intent(in) :: sigjz(sjzlo:sjzhi) - real(amrex_real), intent(in) :: sigsjx(ssjxlo:ssjxhi) - real(amrex_real), intent(in) :: sigsjy(ssjylo:ssjyhi) - real(amrex_real), intent(in) :: sigsjz(ssjzlo:ssjzhi) - - integer :: i,j,k - - do k = tjxlo(3), tjxhi(3) - do j = tjylo(2), tjyhi(2) - do i = tjxlo(1), tjxhi(1) - jx(i,j,k) = jx(i,j,k) * sigsjx(i) * sigjy(j) * sigjz(k) - end do - end do - end do - - do k = tjxlo(3), tjxhi(3) - do j = tjylo(2), tjyhi(2) - do i = tjxlo(1), tjxhi(1) - jy(i,j,k) = jy(i,j,k) * sigjx(i) * sigsjy(j) * sigjz(k) - end do - end do - end do - - do k = tjxlo(3), tjxhi(3) - do j = tjylo(2), tjyhi(2) - do i = tjxlo(1), tjxhi(1) - jz(i,j,k) = jz(i,j,k) * sigjx(i) * sigjy(j) * sigsjz(k) - end do - end do - end do - - end subroutine warpx_dampJ_pml_3d - subroutine warpx_damp_pml_3d (texlo, texhi, teylo, teyhi, tezlo, tezhi, & & tbxlo, tbxhi, tbylo, tbyhi, tbzlo, tbzhi, & @@ -1374,4 +1006,89 @@ contains end do end subroutine warpx_damp_pml_f_3d + subroutine warpx_dampJ_pml_2d (tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & + & jx, jxlo, jxhi, jy, jylo, jyhi, jz, jzlo, jzhi, & + & sigjx, sjxlo, sjxhi, sigjz, sjzlo, sjzhi, & + & sigsjx, ssjxlo, ssjxhi, sigsjz, ssjzlo, ssjzhi) & + bind(c,name='warpx_dampJ_pml_2d') + integer, dimension(2), intent(in) :: tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & + jxlo, jxhi, jylo, jyhi, jzlo, jzhi + integer, intent(in), value :: sjxlo, sjxhi, sjzlo, sjzhi, ssjxlo, ssjxhi, ssjzlo, ssjzhi + real(amrex_real), intent(inout) :: jx(jxlo(1):jxhi(1),jxlo(2):jxhi(2)) + real(amrex_real), intent(inout) :: jy(jylo(1):jyhi(1),jylo(2):jyhi(2)) + real(amrex_real), intent(inout) :: jz(jzlo(1):jzhi(1),jzlo(2):jzhi(2)) + real(amrex_real), intent(in) :: sigjx(sjxlo:sjxhi) + real(amrex_real), intent(in) :: sigjz(sjzlo:sjzhi) + real(amrex_real), intent(in) :: sigsjx(ssjxlo:ssjxhi) + real(amrex_real), intent(in) :: sigsjz(ssjzlo:ssjzhi) + + integer :: i,k + + do k = tjxlo(2), tjxhi(2) + do i = tjxlo(1), tjxhi(1) + jx(i,k) = jx(i,k) * sigsjx(i) * sigjz(k) + end do + end do + + ! do k = tjylo(2), tjyhi(2) + ! do i = tjylo(1), tjyhi(1) + ! jy(i,k) = jy(i,k) + ! end do + ! end do + + do k = tjzlo(2), tjzhi(2) + do i = tjzlo(1), tjzhi(1) + jz(i,k) = jz(i,k) * sigjx(i) * sigsjz(k) + end do + end do + + + end subroutine warpx_dampJ_pml_2d + + subroutine warpx_dampJ_pml_3d (tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & + & jx, jxlo, jxhi, jy, jylo, jyhi, jz, jzlo, jzhi, & + & sigjx, sjxlo, sjxhi, sigjy, sjylo, sjyhi, sigjz, sjzlo, sjzhi, & + & sigsjx, ssjxlo, ssjxhi, sigsjy, ssjylo, ssjyhi, sigsjz, ssjzlo, ssjzhi) & + bind(c,name='warpx_dampJ_pml_3d') + integer, dimension(3), intent(in) :: tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & + jxlo, jxhi, jylo, jyhi, jzlo, jzhi + integer, intent(in), value :: sjxlo, sjxhi, sjylo, sjyhi, sjzlo, sjzhi, ssjxlo, ssjxhi, ssjylo, ssjyhi, ssjzlo, ssjzhi + real(amrex_real), intent(inout) :: jx(jxlo(1):jxhi(1),jxlo(2):jxhi(2),jxlo(3):jxhi(3)) + real(amrex_real), intent(inout) :: jy(jylo(1):jyhi(1),jylo(2):jyhi(2),jylo(3):jyhi(3)) + real(amrex_real), intent(inout) :: jz(jzlo(1):jzhi(1),jzlo(2):jzhi(2),jzlo(3):jzhi(3)) + real(amrex_real), intent(in) :: sigjx(sjxlo:sjxhi) + real(amrex_real), intent(in) :: sigjy(sjylo:sjyhi) + real(amrex_real), intent(in) :: sigjz(sjzlo:sjzhi) + real(amrex_real), intent(in) :: sigsjx(ssjxlo:ssjxhi) + real(amrex_real), intent(in) :: sigsjy(ssjylo:ssjyhi) + real(amrex_real), intent(in) :: sigsjz(ssjzlo:ssjzhi) + + integer :: i,j,k + + do k = tjxlo(3), tjxhi(3) + do j = tjylo(2), tjyhi(2) + do i = tjxlo(1), tjxhi(1) + jx(i,j,k) = jx(i,j,k) * sigsjx(i) * sigjy(j) * sigjz(k) + end do + end do + end do + + do k = tjxlo(3), tjxhi(3) + do j = tjylo(2), tjyhi(2) + do i = tjxlo(1), tjxhi(1) + jy(i,j,k) = jy(i,j,k) * sigjx(i) * sigsjy(j) * sigjz(k) + end do + end do + end do + + do k = tjxlo(3), tjxhi(3) + do j = tjylo(2), tjyhi(2) + do i = tjxlo(1), tjxhi(1) + jz(i,j,k) = jz(i,j,k) * sigjx(i) * sigjy(j) * sigsjz(k) + end do + end do + end do + + end subroutine warpx_dampJ_pml_3d + end module warpx_pml_module diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 011d9cfcc..764688be0 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -11,6 +11,8 @@ #include #endif +#include + #ifdef BL_USE_SENSEI_INSITU #include #endif @@ -476,7 +478,9 @@ 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); - const auto& pml_type = sigba[mfi].pml_type_array; + auto const& pml_Exfab = pml_E[0]->array(mfi); + auto const& pml_Eyfab = pml_E[1]->array(mfi); + auto const& pml_Ezfab = pml_E[2]->array(mfi); WRPX_PUSH_PML_EVEC( tex.loVect(), tex.hiVect(), @@ -488,23 +492,33 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), -#if (AMREX_SPACEDIM==2) - BL_TO_FORTRAN_3D((*pml_j[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_j[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_j[2])[mfi]), - &pml_has_particles, - &mu_c2_dt, -#endif -#if (AMREX_SPACEDIM==3) - BL_TO_FORTRAN_3D((*pml_j[0])[mfi]), - BL_TO_FORTRAN_3D((*pml_j[1])[mfi]), - BL_TO_FORTRAN_3D((*pml_j[2])[mfi]), - &pml_has_particles, - &pml_type, - WRPX_PML_SIGMACOEFF_TO_FORTRAN(sigba[mfi]), - &mu_c2_dt, -#endif - &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); + &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); + + if (pml_has_particles) { + // Update the E field in the PML, using the current + // deposited by the particles in the PML + auto const& pml_jxfab = pml_j[0]->array(mfi); + auto const& pml_jyfab = pml_j[1]->array(mfi); + auto const& pml_jzfab = pml_j[2]->array(mfi); + const Real* sigmaj_x = sigba[mfi].sigma[0].data(); + const Real* sigmaj_y = sigba[mfi].sigma[1].data(); + const Real* sigmaj_z = sigba[mfi].sigma[2].data(); + amrex::ParallelFor( tex, tey, tez, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + warpx_push_ex_pml_current_nodal(i,j,k, + pml_Exfab, pml_jxfab, sigmaj_y, sigmaj_z, mu_c2_dt); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + warpx_push_ey_pml_current_nodal(i,j,k, + pml_Eyfab, pml_jyfab, sigmaj_x, sigmaj_z, mu_c2_dt); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + warpx_push_ez_pml_current_nodal(i,j,k, + pml_Ezfab, pml_jzfab, sigmaj_x, sigmaj_y, mu_c2_dt); + } + ); + } + if (pml_F) { WRPX_PUSH_PML_EVEC_F( @@ -515,7 +529,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), BL_TO_FORTRAN_3D((*pml_F )[mfi]), - &dtsdx_c2, &dtsdy_c2, &dtsdz_c2, + &dtsdx_c2, &dtsdy_c2, &dtsdz_c2, &WarpX::maxwell_fdtd_solver_id); } } diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index f3bce6223..28724eae6 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -1,4 +1,3 @@ - #include #ifdef __cplusplus @@ -329,27 +328,6 @@ extern "C" const BL_FORT_FAB_ARG_3D(bx), const BL_FORT_FAB_ARG_3D(by), const BL_FORT_FAB_ARG_3D(bz), -#if (AMREX_SPACEDIM == 2) - BL_FORT_FAB_ARG_3D(jx), - BL_FORT_FAB_ARG_3D(jy), - BL_FORT_FAB_ARG_3D(jz), - const int* flag, - const amrex::Real* mudt, -#endif -#if (AMREX_SPACEDIM == 3) - BL_FORT_FAB_ARG_3D(jx), - BL_FORT_FAB_ARG_3D(jy), - BL_FORT_FAB_ARG_3D(jz), - const int* flag, - const std::array* pml_type,//const std::string pml_type, - const amrex::Real* sigjx, const int* sigjx_lo, const int* sigjx_hi, - const amrex::Real* sigjy, const int* sigjy_lo, const int* sigjy_hi, - const amrex::Real* sigjz, const int* sigjz_lo, const int* sigjz_hi, - const amrex::Real* sigsjx, const int* sigsjx_lo, const int* sigsjx_hi, - const amrex::Real* sigsjy, const int* sigsjy_lo, const int* sigsjy_hi, - const amrex::Real* sigjsz, const int* sigsjz_lo, const int* sigsjz_hi, - const amrex::Real* mudt, -#endif const amrex::Real* dtsdx, const amrex::Real* dtsdy, const amrex::Real* dtsdz); -- cgit v1.2.3 From 1baaaea38f3ad1561970d53d239b93213d377435 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 28 Aug 2019 09:29:13 -0700 Subject: remove unnecessary references and clean function names --- Source/BoundaryConditions/PML_current.H | 63 ++++++++++++++++---------------- Source/Evolve/WarpXEvolveEM.cpp | 27 ++++++++------ Source/FieldSolver/WarpXPushFieldsEM.cpp | 6 +-- 3 files changed, 49 insertions(+), 47 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML_current.H b/Source/BoundaryConditions/PML_current.H index 38f0e652e..5d23b1554 100644 --- a/Source/BoundaryConditions/PML_current.H +++ b/Source/BoundaryConditions/PML_current.H @@ -6,11 +6,11 @@ using namespace amrex; AMREX_GPU_HOST_DEVICE AMREX_INLINE -void warpx_push_ex_pml_current_nodal (int j, int k, int l, Array4 const& Ex, - Array4 const& jx, - Real const* const& sigjy, - Real const* const& sigjz, - Real mu_c2_dt) +void push_ex_pml_current (int j, int k, int l, Array4 const& Ex, + Array4 const& jx, + Real const * const sigjy, + Real const * const sigjz, + Real mu_c2_dt) { #if (AMREX_SPACEDIM == 3) Real alpha_xy, alpha_xz; @@ -30,14 +30,14 @@ void warpx_push_ex_pml_current_nodal (int j, int k, int l, Array4 const& E } AMREX_GPU_HOST_DEVICE AMREX_INLINE -void warpx_push_ey_pml_current_nodal (int j, int k, int l, Array4 const& Ey, - Array4 const& jy, - Real const* const& sigjx, - Real const* const& sigjz, - Real mu_c2_dt) +void push_ey_pml_current (int j, int k, int l, Array4 const& Ey, + Array4 const& jy, + Real const * const sigjx, + Real const * const sigjz, + Real mu_c2_dt) { #if (AMREX_SPACEDIM == 3) -Real alpha_yx, alpha_yz; + Real alpha_yx, alpha_yz; if (sigjx[j]+sigjz[l] == 0){ alpha_yx = 0.5; alpha_yz = 0.5; @@ -55,11 +55,11 @@ Real alpha_yx, alpha_yz; } AMREX_GPU_HOST_DEVICE AMREX_INLINE -void warpx_push_ez_pml_current_nodal (int j, int k, int l, Array4 const& Ez, - Array4 const& jz, - Real const* const& sigjx, - Real const* const& sigjy, - Real mu_c2_dt) +void push_ez_pml_current (int j, int k, int l, Array4 const& Ez, + Array4 const& jz, + Real const * const sigjx, + Real const * const sigjy, + Real mu_c2_dt) { #if (AMREX_SPACEDIM == 3) Real alpha_zx, alpha_zy; @@ -78,13 +78,12 @@ void warpx_push_ez_pml_current_nodal (int j, int k, int l, Array4 const& E #endif } - AMREX_GPU_HOST_DEVICE AMREX_INLINE -void warpx_damp_jx_pml (int j, int k, int l, - Array4 const& jx, - Real const* const& sigsjx, - Real const* const& sigjy, - Real const* const& sigjz) +void damp_jx_pml (int j, int k, int l, + Array4 const& jx, + Real const* const sigsjx, + Real const* const sigjy, + Real const* const sigjz) { #if (AMREX_SPACEDIM == 3) jx(j,k,l) = jx(j,k,l) * sigsjx[j] * sigjy[k] * sigjz[l]; @@ -94,11 +93,11 @@ void warpx_damp_jx_pml (int j, int k, int l, } AMREX_GPU_HOST_DEVICE AMREX_INLINE -void warpx_damp_jy_pml (int j, int k, int l, - Array4 const& jy, - Real const* const& sigjx, - Real const* const& sigsjy, - Real const* const& sigjz) +void damp_jy_pml (int j, int k, int l, + Array4 const& jy, + Real const * const sigjx, + Real const * const sigsjy, + Real const * const sigjz) { #if (AMREX_SPACEDIM == 3) jy(j,k,l) = jy(j,k,l) * sigjx[j] * sigsjy[k] * sigjz[l]; @@ -108,11 +107,11 @@ void warpx_damp_jy_pml (int j, int k, int l, } AMREX_GPU_HOST_DEVICE AMREX_INLINE -void warpx_damp_jz_pml (int j, int k, int l, - Array4 const& jz, - Real const* const& sigjx, - Real const* const& sigjy, - Real const* const& sigsjz) +void damp_jz_pml (int j, int k, int l, + Array4 const& jz, + Real const * const sigjx, + Real const * const sigjy, + Real const * const sigsjz) { #if (AMREX_SPACEDIM == 3) jz(j,k,l) = jz(j,k,l) * sigjx[j] * sigjy[k] * sigsjz[l]; diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index ceccbf789..c32db95e6 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -352,16 +352,16 @@ WarpX::OneStep_nosub (Real cur_time) amrex::ParallelFor( tjx, tjy, tjz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_damp_jx_pml(i,j,k, - pml_jxfab, sigma_star_cum_fac_j_x, sigma_cum_fac_j_y, sigma_cum_fac_j_z); + damp_jx_pml(i, j, k, pml_jxfab, sigma_star_cum_fac_j_x, + sigma_cum_fac_j_y, sigma_cum_fac_j_z); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_damp_jy_pml(i,j,k, - pml_jyfab, sigma_cum_fac_j_x, sigma_star_cum_fac_j_y, sigma_cum_fac_j_z); + damp_jy_pml(i, j, k, pml_jyfab, sigma_cum_fac_j_x, + sigma_star_cum_fac_j_y, sigma_cum_fac_j_z); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_damp_jz_pml(i,j,k, - pml_jzfab, sigma_cum_fac_j_x, sigma_cum_fac_j_y, sigma_star_cum_fac_j_z); + damp_jz_pml(i, j, k, pml_jzfab, sigma_cum_fac_j_x, + sigma_cum_fac_j_y, sigma_star_cum_fac_j_z); } ); } @@ -391,16 +391,19 @@ WarpX::OneStep_nosub (Real cur_time) amrex::ParallelFor( tjx, tjy, tjz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_damp_jx_pml(i,j,k, - pml_jxfab, sigma_star_cum_fac_j_x, sigma_cum_fac_j_y, sigma_cum_fac_j_z); + damp_jx_pml(i, j, k, pml_jxfab, + sigma_star_cum_fac_j_x, + sigma_cum_fac_j_y, sigma_cum_fac_j_z); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_damp_jy_pml(i,j,k, - pml_jyfab, sigma_cum_fac_j_x, sigma_star_cum_fac_j_y, sigma_cum_fac_j_z); + damp_jy_pml(i, j, k, pml_jyfab, sigma_cum_fac_j_x, + sigma_star_cum_fac_j_y, + sigma_cum_fac_j_z); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_damp_jz_pml(i,j,k, - pml_jzfab, sigma_cum_fac_j_x, sigma_cum_fac_j_y, sigma_star_cum_fac_j_z); + damp_jz_pml(i, j, k, pml_jzfab, sigma_cum_fac_j_x, + sigma_cum_fac_j_y, + sigma_star_cum_fac_j_z); } ); } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 764688be0..612dc4e0f 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -505,15 +505,15 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) const Real* sigmaj_z = sigba[mfi].sigma[2].data(); amrex::ParallelFor( tex, tey, tez, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_ex_pml_current_nodal(i,j,k, + push_ex_pml_current(i,j,k, pml_Exfab, pml_jxfab, sigmaj_y, sigmaj_z, mu_c2_dt); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_ey_pml_current_nodal(i,j,k, + push_ey_pml_current(i,j,k, pml_Eyfab, pml_jyfab, sigmaj_x, sigmaj_z, mu_c2_dt); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - warpx_push_ez_pml_current_nodal(i,j,k, + push_ez_pml_current(i,j,k, pml_Ezfab, pml_jzfab, sigmaj_x, sigmaj_y, mu_c2_dt); } ); -- cgit v1.2.3 From a99bdf09f60dccb4077a3cccb018a719b229cab7 Mon Sep 17 00:00:00 2001 From: ablelly Date: Wed, 28 Aug 2019 23:30:24 +0200 Subject: Added shift for the sigma functions --- Source/BoundaryConditions/PML_current.H | 42 ++++++++++++++++------------ Source/BoundaryConditions/WarpXEvolvePML.cpp | 29 ++++++++++++++++--- Source/FieldSolver/WarpXPushFieldsEM.cpp | 23 +++++++++++---- 3 files changed, 67 insertions(+), 27 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML_current.H b/Source/BoundaryConditions/PML_current.H index 3e45fe3a5..910186a96 100644 --- a/Source/BoundaryConditions/PML_current.H +++ b/Source/BoundaryConditions/PML_current.H @@ -10,17 +10,18 @@ void push_ex_pml_current (int j, int k, int l, Array4 const& Ex, Array4 const& jx, Real const * const sigjy, Real const * const sigjz, + int ylo, int zlo, Real mu_c2_dt) { #if (AMREX_SPACEDIM == 3) Real alpha_xy, alpha_xz; - if (sigjy[k]+sigjz[l] == 0){ + if (sigjy[k-ylo]+sigjz[l-zlo] == 0){ alpha_xy = 0.5; alpha_xz = 0.5; } else { - alpha_xy = sigjy[k]/(sigjy[k]+sigjz[l]); - alpha_xz = sigjz[l]/(sigjy[k]+sigjz[l]); + alpha_xy = sigjy[k-ylo]/(sigjy[k-ylo]+sigjz[l-zlo]); + alpha_xz = sigjz[l-zlo]/(sigjy[k-ylo]+sigjz[l-zlo]); } Ex(j,k,l,0) = Ex(j,k,l,0) - mu_c2_dt * alpha_xy * jx(j,k,l); Ex(j,k,l,1) = Ex(j,k,l,1) - mu_c2_dt * alpha_xz * jx(j,k,l); @@ -34,17 +35,18 @@ void push_ey_pml_current (int j, int k, int l, Array4 const& Ey, Array4 const& jy, Real const * const sigjx, Real const * const sigjz, + int xlo, int zlo, Real mu_c2_dt) { #if (AMREX_SPACEDIM == 3) Real alpha_yx, alpha_yz; - if (sigjx[j]+sigjz[l] == 0){ + if (sigjx[j-xlo]+sigjz[l-zlo] == 0){ alpha_yx = 0.5; alpha_yz = 0.5; } else { - alpha_yx = sigjx[j]/(sigjx[j]+sigjz[l]); - alpha_yz = sigjz[l]/(sigjx[j]+sigjz[l]); + alpha_yx = sigjx[j-xlo]/(sigjx[j-xlo]+sigjz[l-zlo]); + alpha_yz = sigjz[l-zlo]/(sigjx[j-xlo]+sigjz[l-zlo]); } Ey(j,k,l,0) = Ey(j,k,l,0) - mu_c2_dt * alpha_yx * jy(j,k,l); Ey(j,k,l,1) = Ey(j,k,l,1) - mu_c2_dt * alpha_yz * jy(j,k,l); @@ -59,17 +61,18 @@ void push_ez_pml_current (int j, int k, int l, Array4 const& Ez, Array4 const& jz, Real const * const sigjx, Real const * const sigjy, + int xlo, int ylo, Real mu_c2_dt) { #if (AMREX_SPACEDIM == 3) Real alpha_zx, alpha_zy; - if (sigjx[j]+sigjy[k]==0){ + if (sigjx[j-xlo]+sigjy[k-ylo]==0){ alpha_zx = 0.5; alpha_zy = 0.5; } else { - alpha_zx = sigjx[j]/(sigjx[j]+sigjy[k]); - alpha_zy = sigjy[k]/(sigjx[j]+sigjy[k]); + alpha_zx = sigjx[j-xlo]/(sigjx[j-xlo]+sigjy[k-ylo]); + alpha_zy = sigjy[k-ylo]/(sigjx[j-xlo]+sigjy[k-ylo]); } Ez(j,k,l,0) = Ez(j,k,l,0) - mu_c2_dt * alpha_zx * jz(j,k,l); Ez(j,k,l,1) = Ez(j,k,l,1) - mu_c2_dt * alpha_zy * jz(j,k,l); @@ -83,12 +86,13 @@ void damp_jx_pml (int j, int k, int l, Array4 const& jx, Real const* const sigsjx, Real const* const sigjy, - Real const* const sigjz) + Real const* const sigjz, + int xlo, int ylo, int zlo) { #if (AMREX_SPACEDIM == 3) - jx(j,k,l) = jx(j,k,l) * sigsjx[j] * sigjy[k] * sigjz[l]; + jx(j,k,l) = jx(j,k,l) * sigsjx[j-xlo] * sigjy[k-ylo] * sigjz[l-zlo]; #else - jx(j,k,l) = jx(j,k,l) * sigsjx[j] * sigjz[k]; + jx(j,k,l) = jx(j,k,l) * sigsjx[j-xlo] * sigjz[k-zlo]; #endif } @@ -97,12 +101,13 @@ void damp_jy_pml (int j, int k, int l, Array4 const& jy, Real const * const sigjx, Real const * const sigsjy, - Real const * const sigjz) + Real const * const sigjz, + int xlo, int ylo, int zlo) { #if (AMREX_SPACEDIM == 3) - jy(j,k,l) = jy(j,k,l) * sigjx[j] * sigsjy[k] * sigjz[l]; + jy(j,k,l) = jy(j,k,l) * sigjx[j-xlo] * sigsjy[k-ylo] * sigjz[l-zlo]; #else - jy(j,k,l) = jy(j,k,l) * sigjx[j] * sigjz[k]; + jy(j,k,l) = jy(j,k,l) * sigjx[j-xlo] * sigjz[k-zlo]; #endif } @@ -111,12 +116,13 @@ void damp_jz_pml (int j, int k, int l, Array4 const& jz, Real const * const sigjx, Real const * const sigjy, - Real const * const sigsjz) + Real const * const sigsjz, + int xlo, int ylo, int zlo) { #if (AMREX_SPACEDIM == 3) - jz(j,k,l) = jz(j,k,l) * sigjx[j] * sigjy[k] * sigsjz[l]; + jz(j,k,l) = jz(j,k,l) * sigjx[j-xlo] * sigjy[k-ylo] * sigsjz[l-zlo]; #else - jz(j,k,l) = jz(j,k,l) * sigjx[j] * sigsjz[k]; + jz(j,k,l) = jz(j,k,l) * sigjx[j-xlo] * sigsjz[k-zlo]; #endif } diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp index 1e5519451..a801aabe4 100644 --- a/Source/BoundaryConditions/WarpXEvolvePML.cpp +++ b/Source/BoundaryConditions/WarpXEvolvePML.cpp @@ -112,7 +112,7 @@ WarpX::DampJPML (int lev, PatchType patch_type) const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp(); const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp() : pml[lev]->GetMultiSigmaBox_cp(); - + #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -138,18 +138,39 @@ WarpX::DampJPML (int lev, PatchType patch_type) const Box& tjy = mfi.tilebox(jy_nodal_flag); const Box& tjz = mfi.tilebox(jz_nodal_flag); + auto const& AMREX_RESTRICT x_lo = sigba[mfi].sigma_cum_fac[0].lo(); +#if (AMREX_SPACEDIM == 3) + auto const& AMREX_RESTRICT y_lo = sigba[mfi].sigma_cum_fac[1].lo(); + auto const& AMREX_RESTRICT z_lo = sigba[mfi].sigma_cum_fac[2].lo(); +#else + int y_lo = 0; + auto const& AMREX_RESTRICT z_lo = sigba[mfi].sigma_cum_fac[1].lo(); +#endif + + auto const& AMREX_RESTRICT xs_lo = sigba[mfi].sigma_star_cum_fac[0].lo(); +#if (AMREX_SPACEDIM == 3) + auto const& AMREX_RESTRICT ys_lo = sigba[mfi].sigma_star_cum_fac[1].lo(); + auto const& AMREX_RESTRICT zs_lo = sigba[mfi].sigma_star_cum_fac[2].lo(); +#else + int ys_lo = 0; + auto const& AMREX_RESTRICT zs_lo = sigba[mfi].sigma_star_cum_fac[1].lo(); +#endif + amrex::ParallelFor( tjx, tjy, tjz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { damp_jx_pml(i, j, k, pml_jxfab, sigma_star_cum_fac_j_x, - sigma_cum_fac_j_y, sigma_cum_fac_j_z); + sigma_cum_fac_j_y, sigma_cum_fac_j_z, + xs_lo,y_lo, z_lo); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { damp_jy_pml(i, j, k, pml_jyfab, sigma_cum_fac_j_x, - sigma_star_cum_fac_j_y, sigma_cum_fac_j_z); + sigma_star_cum_fac_j_y, sigma_cum_fac_j_z, + x_lo,ys_lo, z_lo); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { damp_jz_pml(i, j, k, pml_jzfab, sigma_cum_fac_j_x, - sigma_cum_fac_j_y, sigma_star_cum_fac_j_z); + sigma_cum_fac_j_y, sigma_star_cum_fac_j_z, + x_lo,y_lo, zs_lo); } ); } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 612dc4e0f..4dc31f994 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -493,7 +493,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); - + if (pml_has_particles) { // Update the E field in the PML, using the current // deposited by the particles in the PML @@ -503,18 +503,31 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) const Real* sigmaj_x = sigba[mfi].sigma[0].data(); const Real* sigmaj_y = sigba[mfi].sigma[1].data(); const Real* sigmaj_z = sigba[mfi].sigma[2].data(); - amrex::ParallelFor( tex, tey, tez, + + auto const& AMREX_RESTRICT x_lo = sigba[mfi].sigma[0].lo(); +#if (AMREX_SPACEDIM == 3) + auto const& AMREX_RESTRICT y_lo = sigba[mfi].sigma[1].lo(); + auto const& AMREX_RESTRICT z_lo = sigba[mfi].sigma[2].lo(); +#else + int y_lo = 0; + auto const& AMREX_RESTRICT z_lo = sigba[mfi].sigma[1].lo(); +#endif + + amrex::ParallelFor( tex, tey, tez, [=] AMREX_GPU_DEVICE (int i, int j, int k) { push_ex_pml_current(i,j,k, - pml_Exfab, pml_jxfab, sigmaj_y, sigmaj_z, mu_c2_dt); + pml_Exfab, pml_jxfab, sigmaj_y, sigmaj_z, + y_lo, z_lo, mu_c2_dt); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { push_ey_pml_current(i,j,k, - pml_Eyfab, pml_jyfab, sigmaj_x, sigmaj_z, mu_c2_dt); + pml_Eyfab, pml_jyfab, sigmaj_x, sigmaj_z, + x_lo, z_lo, mu_c2_dt); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { push_ez_pml_current(i,j,k, - pml_Ezfab, pml_jzfab, sigmaj_x, sigmaj_y, mu_c2_dt); + pml_Ezfab, pml_jzfab, sigmaj_x, sigmaj_y, + x_lo, y_lo, mu_c2_dt); } ); } -- cgit v1.2.3 From e86dc1d80ff4ada735f23e03372b8a5d2fc3499b Mon Sep 17 00:00:00 2001 From: ablelly Date: Thu, 29 Aug 2019 17:43:56 +0200 Subject: layout corrections --- Source/BoundaryConditions/WarpXEvolvePML.cpp | 2 -- Source/FieldSolver/WarpXPushFieldsEM.cpp | 4 ++-- 2 files changed, 2 insertions(+), 4 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp index 505fa83d2..813288842 100644 --- a/Source/BoundaryConditions/WarpXEvolvePML.cpp +++ b/Source/BoundaryConditions/WarpXEvolvePML.cpp @@ -12,9 +12,7 @@ #include #endif -#ifndef PML_CURRENT_H_ #include -#endif using namespace amrex; diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 4dc31f994..94e49f0bc 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -492,7 +492,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), - &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); + &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); if (pml_has_particles) { // Update the E field in the PML, using the current @@ -542,7 +542,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) BL_TO_FORTRAN_3D((*pml_E[1])[mfi]), BL_TO_FORTRAN_3D((*pml_E[2])[mfi]), BL_TO_FORTRAN_3D((*pml_F )[mfi]), - &dtsdx_c2, &dtsdy_c2, &dtsdz_c2, + &dtsdx_c2, &dtsdy_c2, &dtsdz_c2, &WarpX::maxwell_fdtd_solver_id); } } -- cgit v1.2.3 From a025c00a7436299a754879edae0d70abe5fec1d3 Mon Sep 17 00:00:00 2001 From: ablelly Date: Thu, 29 Aug 2019 22:32:01 +0200 Subject: Removed useless variables in PML Exchange functions --- Source/BoundaryConditions/PML.H | 28 ++++----------- Source/BoundaryConditions/PML.cpp | 62 ++++++++++++++------------------ Source/FieldSolver/WarpXPushFieldsEM.cpp | 2 +- Source/Parallelization/WarpXComm.cpp | 27 +++++--------- 4 files changed, 44 insertions(+), 75 deletions(-) (limited to 'Source/FieldSolver/WarpXPushFieldsEM.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index ad7326dbb..0a9b23b7f 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -134,33 +134,21 @@ public: #endif void ExchangeB (const std::array& B_fp, - const std::array& B_cp, int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), - const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); + const std::array& B_cp, int do_pml_in_domain); void ExchangeE (const std::array& E_fp, - const std::array& E_cp, int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), - const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); + const std::array& E_cp, int do_pml_in_domain); void CopyJtoPMLs (const std::array& j_fp, const std::array& j_cp); void ExchangeB (PatchType patch_type, - const std::array& Bp, int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), - const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); + const std::array& Bp, int do_pml_in_domain); void ExchangeE (PatchType patch_type, - const std::array& Ep, int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), - const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); + const std::array& Ep, int do_pml_in_domain); void CopyJtoPMLs (PatchType patch_type, const std::array& jp); - void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp, int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), - const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); - void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp, int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), - const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); + void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp, int do_pml_in_domain); + void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp, int do_pml_in_domain); void FillBoundary (); void FillBoundaryE (); @@ -205,9 +193,7 @@ private: const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); - static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom, int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), - const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); + static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom, int do_pml_in_domain); static void CopyToPML (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); }; diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 8c35e34e2..7d1906525 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -638,57 +638,53 @@ PML::GetF_cp () void PML::ExchangeB (const std::array& B_fp, - const std::array& B_cp, int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) + const std::array& B_cp, int do_pml_in_domain) { - ExchangeB(PatchType::fine, B_fp, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); - ExchangeB(PatchType::coarse, B_cp, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); + ExchangeB(PatchType::fine, B_fp, do_pml_in_domain); + ExchangeB(PatchType::coarse, B_cp, do_pml_in_domain); } void PML::ExchangeB (PatchType patch_type, - const std::array& Bp, int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) + const std::array& Bp, int do_pml_in_domain) { if (patch_type == PatchType::fine && pml_B_fp[0] && Bp[0]) { - Exchange(*pml_B_fp[0], *Bp[0], *m_geom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); - Exchange(*pml_B_fp[1], *Bp[1], *m_geom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); - Exchange(*pml_B_fp[2], *Bp[2], *m_geom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); + Exchange(*pml_B_fp[0], *Bp[0], *m_geom, do_pml_in_domain); + Exchange(*pml_B_fp[1], *Bp[1], *m_geom, do_pml_in_domain); + Exchange(*pml_B_fp[2], *Bp[2], *m_geom, do_pml_in_domain); } else if (patch_type == PatchType::coarse && pml_B_cp[0] && Bp[0]) { - Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); - Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); - Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); + Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom, do_pml_in_domain); + Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom, do_pml_in_domain); + Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom, do_pml_in_domain); } } void PML::ExchangeE (const std::array& E_fp, - const std::array& E_cp, int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) + const std::array& E_cp, int do_pml_in_domain) { - ExchangeE(PatchType::fine, E_fp, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); - ExchangeE(PatchType::coarse, E_cp, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); + ExchangeE(PatchType::fine, E_fp, do_pml_in_domain); + ExchangeE(PatchType::coarse, E_cp, do_pml_in_domain); } void PML::ExchangeE (PatchType patch_type, - const std::array& Ep, int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) + const std::array& Ep, int do_pml_in_domain) { if (patch_type == PatchType::fine && pml_E_fp[0] && Ep[0]) { - Exchange(*pml_E_fp[0], *Ep[0], *m_geom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); - Exchange(*pml_E_fp[1], *Ep[1], *m_geom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); - Exchange(*pml_E_fp[2], *Ep[2], *m_geom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); + Exchange(*pml_E_fp[0], *Ep[0], *m_geom, do_pml_in_domain); + Exchange(*pml_E_fp[1], *Ep[1], *m_geom, do_pml_in_domain); + Exchange(*pml_E_fp[2], *Ep[2], *m_geom, do_pml_in_domain); } else if (patch_type == PatchType::coarse && pml_E_cp[0] && Ep[0]) { - Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); - Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); - Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); + Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom, do_pml_in_domain); + Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom, do_pml_in_domain); + Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom, do_pml_in_domain); } } @@ -720,30 +716,26 @@ PML::CopyJtoPMLs (const std::array& j_fp, void -PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp, int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) +PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp, int do_pml_in_domain) { - ExchangeF(PatchType::fine, F_fp, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); - ExchangeF(PatchType::coarse, F_cp, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); + ExchangeF(PatchType::fine, F_fp, do_pml_in_domain); + ExchangeF(PatchType::coarse, F_cp, do_pml_in_domain); } void -PML::ExchangeF (PatchType patch_type, MultiFab* Fp, int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) +PML::ExchangeF (PatchType patch_type, MultiFab* Fp, int do_pml_in_domain) { if (patch_type == PatchType::fine && pml_F_fp && Fp) { - Exchange(*pml_F_fp, *Fp, *m_geom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); + Exchange(*pml_F_fp, *Fp, *m_geom, do_pml_in_domain); } else if (patch_type == PatchType::coarse && pml_F_cp && Fp) { - Exchange(*pml_F_cp, *Fp, *m_cgeom, do_pml_in_domain, ncell, do_pml_Lo, do_pml_Hi); + Exchange(*pml_F_cp, *Fp, *m_cgeom, do_pml_in_domain); } } void PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom, - int do_pml_in_domain, int ncell, - const amrex::IntVect do_pml_Lo, - const amrex::IntVect do_pml_Hi) + int do_pml_in_domain) { if (do_pml_in_domain){ const IntVect& ngr = reg.nGrowVect(); diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 94e49f0bc..dae18002d 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -461,7 +461,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) if (do_pml && pml[lev]->ok()) { - if (F) pml[lev]->ExchangeF(patch_type, F, do_pml_in_domain, pml_ncell); + if (F) pml[lev]->ExchangeF(patch_type, F, do_pml_in_domain); const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp(); const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp(); diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index c539cd89a..feb9d3564 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -19,8 +19,7 @@ WarpX::ExchangeWithPmlB (int lev) { Bfield_cp[lev][0].get(), Bfield_cp[lev][1].get(), Bfield_cp[lev][2].get() }, - do_pml_in_domain, pml_ncell, - do_pml_Lo, do_pml_Hi); + do_pml_in_domain); } } @@ -34,8 +33,7 @@ WarpX::ExchangeWithPmlE (int lev) { Efield_cp[lev][0].get(), Efield_cp[lev][1].get(), Efield_cp[lev][2].get() }, - do_pml_in_domain, pml_ncell, - do_pml_Lo, do_pml_Hi); + do_pml_in_domain); } } @@ -45,8 +43,7 @@ WarpX::ExchangeWithPmlF (int lev) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeF(F_fp[lev].get(), F_cp[lev].get(), - do_pml_in_domain, pml_ncell, - do_pml_Lo, do_pml_Hi); + do_pml_in_domain); } } @@ -259,8 +256,7 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) { Efield_fp[lev][0].get(), Efield_fp[lev][1].get(), Efield_fp[lev][2].get() }, - do_pml_in_domain, pml_ncell, - do_pml_Lo, do_pml_Hi); + do_pml_in_domain); pml[lev]->FillBoundaryE(patch_type); } @@ -276,8 +272,7 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) { Efield_cp[lev][0].get(), Efield_cp[lev][1].get(), Efield_cp[lev][2].get() }, - do_pml_in_domain, pml_ncell, - do_pml_Lo, do_pml_Hi); + do_pml_in_domain); pml[lev]->FillBoundaryE(patch_type); } @@ -305,8 +300,7 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) { Bfield_fp[lev][0].get(), Bfield_fp[lev][1].get(), Bfield_fp[lev][2].get() }, - do_pml_in_domain, pml_ncell, - do_pml_Lo, do_pml_Hi); + do_pml_in_domain); pml[lev]->FillBoundaryB(patch_type); } const auto& period = Geom(lev).periodicity(); @@ -321,8 +315,7 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type) { Bfield_cp[lev][0].get(), Bfield_cp[lev][1].get(), Bfield_cp[lev][2].get() }, - do_pml_in_domain, pml_ncell, - do_pml_Lo, do_pml_Hi); + do_pml_in_domain); pml[lev]->FillBoundaryB(patch_type); } const auto& cperiod = Geom(lev-1).periodicity(); @@ -346,8 +339,7 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeF(patch_type, F_fp[lev].get(), - do_pml_in_domain, pml_ncell, - do_pml_Lo, do_pml_Hi); + do_pml_in_domain); pml[lev]->FillBoundaryF(patch_type); } @@ -359,8 +351,7 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type) if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeF(patch_type, F_cp[lev].get(), - do_pml_in_domain, pml_ncell, - do_pml_Lo, do_pml_Hi); + do_pml_in_domain); pml[lev]->FillBoundaryF(patch_type); } -- cgit v1.2.3