diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index c58542db2..84a44eac3 100755 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -12,17 +12,17 @@ install: stage: install only: - develop - - develop-mat + - mergingPML script: # Force workdir cleaning in case of retried - echo "CI_PIPELINE_ID = " $CI_PIPELINE_ID - env - - if [ -d /sps/gitlab-runner/$CI_PIPELINE_ID ] ; then rm -rf /sps/gitlab-runner/$CI_PIPELINE_ID ; fi + - if [ -d /sps3/gitlab-runner/$CI_PIPELINE_ID ] ; then rm -rf /sps/gitlab-runner/$CI_PIPELINE_ID ; fi # Move in test dir - - mkdir -p /sps/gitlab-runner/$CI_PIPELINE_ID - - cp -r $CI_PROJECT_DIR /sps/gitlab-runner/$CI_PIPELINE_ID - - cd /sps/gitlab-runner/$CI_PIPELINE_ID/smilei + - mkdir -p /sps3/gitlab-runner/$CI_PIPELINE_ID + - cp -r $CI_PROJECT_DIR /sps3/gitlab-runner/$CI_PIPELINE_ID + - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei # Set the environment - make uninstall_happi - make happi @@ -31,91 +31,91 @@ compile_default: stage: compile_default only: - develop - - develop-mat + - mergingPML script: # Move in test dir - - cd /sps/gitlab-runner/$CI_PIPELINE_ID/smilei/validation + - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei/validation # Run validation script - - python validation.py -c -v -l "/sps/gitlab-runner/logs" -t "00:30:00" + - python validation.py -c -v -l "/sps3/gitlab-runner/logs" -t "00:30:00" run1D: stage: run_default only: - develop - - develop-mat + - mergingPML script: # Move in test dir - - cd /sps/gitlab-runner/$CI_PIPELINE_ID/smilei/validation + - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei/validation # Run validation script - - python validation.py -b "tst1d_*_*.py" -m 2 -o 12 -v -r 1 -l "/sps/gitlab-runner/logs" - - python validation.py -b "tst_ionization_current_1d*.py" -m 2 -o 12 -v -r 1 -l "/sps/gitlab-runner/logs" + - python validation.py -b "tst1d_*_*.py" -m 2 -o 12 -v -r 1 -l "/sps3/gitlab-runner/logs" + - python validation.py -b "tst_ionization_current_1d*.py" -m 2 -o 12 -v -r 1 -l "/sps3/gitlab-runner/logs" run2D: stage: run_default only: - develop - - develop-mat + - mergingPML script: # Move in test dir - - cd /sps/gitlab-runner/$CI_PIPELINE_ID/smilei/validation + - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei/validation # Run validation script - - python validation.py -b "tst2d_??_*.py" -m 4 -o 4 -n 1 -v -r 1 -l "/sps/gitlab-runner/logs" - - python validation.py -b "tst2d_s_??_*.py" -m 4 -o 4 -n 1 -v -r 1 -l "/sps/gitlab-runner/logs" - - python validation.py -b "tst2d_v_??_*.py" -m 4 -o 4 -n 1 -v -r 1 -l "/sps/gitlab-runner/logs" - - python validation.py -b "tst_ionization_current_2d*.py" -m 4 -o 4 -n 1 -v -r 1 -l "/sps/gitlab-runner/logs" + - python validation.py -b "tst2d_??_*.py" -m 4 -o 4 -n 1 -v -r 1 -l "/sps3/gitlab-runner/logs" + - python validation.py -b "tst2d_s_??_*.py" -m 4 -o 4 -n 1 -v -r 1 -l "/sps3/gitlab-runner/logs" + - python validation.py -b "tst2d_v_??_*.py" -m 4 -o 4 -n 1 -v -r 1 -l "/sps3/gitlab-runner/logs" + - python validation.py -b "tst_ionization_current_2d*.py" -m 4 -o 4 -n 1 -v -r 1 -l "/sps3/gitlab-runner/logs" run3D: stage: run_default only: - develop - - develop-mat + - mergingPML script: # Move in test dir - - cd /sps/gitlab-runner/$CI_PIPELINE_ID/smilei/validation + - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei/validation # Run validation script - - python validation.py -b "tst3d_??_*.py" -m 8 -o 12 --resource-file resources.json -v -r 1 -l "/sps/gitlab-runner/logs" - - python validation.py -b "tst3d_s_??_*.py" -m 8 -o 12 --resource-file resources.json -v -r 1 -l "/sps/gitlab-runner/logs" - - python validation.py -b "tst3d_v_??_*.py" -m 8 -o 12 --resource-file resources.json -v -r 1 -l "/sps/gitlab-runner/logs" - - python validation.py -b "tst3d_a_??_*.py" -m 8 -o 12 --resource-file resources.json -v -r 1 -l "/sps/gitlab-runner/logs" - - python validation.py -b "tst_ionization_current_3d*.py" -m 4 -o 4 -n 1 -v -r 1 -l "/sps/gitlab-runner/logs" + - python validation.py -b "tst3d_??_*.py" -m 8 -o 12 --resource-file resources.json -v -r 1 -l "/sps3/gitlab-runner/logs" + - python validation.py -b "tst3d_s_??_*.py" -m 8 -o 12 --resource-file resources.json -v -r 1 -l "/sps3/gitlab-runner/logs" + - python validation.py -b "tst3d_v_??_*.py" -m 8 -o 12 --resource-file resources.json -v -r 1 -l "/sps3/gitlab-runner/logs" + - python validation.py -b "tst3d_a_??_*.py" -m 8 -o 12 --resource-file resources.json -v -r 1 -l "/sps3/gitlab-runner/logs" + - python validation.py -b "tst_ionization_current_3d*.py" -m 4 -o 4 -n 1 -v -r 1 -l "/sps3/gitlab-runner/logs" runAM: stage: run_default only: - develop - - develop-mat + - mergingPML script: # Move in test dir - - cd /sps/gitlab-runner/$CI_PIPELINE_ID/smilei/validation + - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei/validation # Run validation script - - python validation.py -b "tstAM_??_*.py" -m 8 -o 12 -v -t 00:30:00 -l "/sps/gitlab-runner/logs" + - python validation.py -b "tstAM_??_*.py" -m 8 -o 12 -v -t 00:30:00 -l "/sps3/gitlab-runner/logs" runCollisions: stage: run_default only: - develop - - develop-mat + - mergingPML script: # Move in test dir - - cd /sps/gitlab-runner/$CI_PIPELINE_ID/smilei/validation + - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei/validation # Run checking script - - python validation.py -b "tst_collisions*.py" -m 4 -o 4 -n 1 -v -r 1 -l "/sps/gitlab-runner/logs" + - python validation.py -b "tst_collisions*.py" -m 4 -o 4 -n 1 -v -r 1 -l "/sps3/gitlab-runner/logs" compile_picsar: stage: compile_picsar only: - develop - - develop-mat + - mergingPML script: - - cd /sps/gitlab-runner/$CI_PIPELINE_ID/smilei + - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei - make clean - python validation/validation.py -k picsar -c -v @@ -123,20 +123,20 @@ run_picsar: stage: run_picsar only: - develop - - develop-mat + - mergingPML script: - - cd /sps/gitlab-runner/$CI_PIPELINE_ID/smilei/validation + - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei/validation - python validation.py -k picsar -b "tstAM_picsar_04_laser_wake.py" -m 4 -o 4 -n 1 -v compile_debug: stage: compile_debug only: - develop - - develop-mat + - mergingPML script: - - cd /sps/gitlab-runner/$CI_PIPELINE_ID/smilei + - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei - make clean - python validation/validation.py -k debug -c -v @@ -144,10 +144,10 @@ compile_no_mpi_threadmultiple: stage: compile_no_mpi_threadmultiple only: - develop - - develop-mat + - mergingPML script: - - cd /sps/gitlab-runner/$CI_PIPELINE_ID/smilei + - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei - make clean - python validation/validation.py -k no_mpi_tm -c -v @@ -155,10 +155,10 @@ compile_no_openmp: stage: compile_no_openmp only: - develop - - develop-mat + - mergingPML script: - - cd /sps/gitlab-runner/$CI_PIPELINE_ID/smilei + - cd /sps3/gitlab-runner/$CI_PIPELINE_ID/smilei - make clean - python validation/validation.py -k noopenmp -c -v diff --git a/benchmarks/tst2d_18_em_pml.py b/benchmarks/tst2d_18_em_pml.py new file mode 100755 index 000000000..18dc5ee09 --- /dev/null +++ b/benchmarks/tst2d_18_em_pml.py @@ -0,0 +1,46 @@ +# ---------------------------------------------------------------------------------------- +# SIMULATION PARAMETERS FOR THE PIC-CODE SMILEI +# ---------------------------------------------------------------------------------------- + +import numpy as np + +l0 = 2.*np.pi # laser wavelength +t0 = l0 # optical cycle +Lsim = [20.*l0,20.*l0] # length of the simulation +Tsim = 30.*t0 # duration of the simulation +resx = 8. # nb of cells in on laser wavelength +rest = 1.1*np.sqrt(2)*resx # time of timestep in one optical cycle + +Main( + geometry = "2Dcartesian", + maxwell_solver = "Yee", + interpolation_order = 2 , + cell_length = [l0/resx,l0/resx], + grid_length = Lsim, + number_of_patches = [ 16, 16 ], + timestep = t0/rest, + simulation_time = Tsim, + EM_boundary_conditions = [ + ['PML','PML'], + ['PML','PML'], + ], + number_of_pml_cells = [[10,10],[10,10]], + random_seed = smilei_mpi_rank +) + +Antenna( + field='Jz', + time_profile = lambda t: np.sin(2.*np.pi*t/t0)*np.sin(2.*np.pi*t/(4*t0))*(1.-np.heaviside(t-2*t0, 1)), + space_profile=gaussian(1., xfwhm=l0, yfwhm=l0, xcenter=Main.grid_length[0]*0.5, ycenter=Main.grid_length[1]*0.5) +) + +globalEvery=1. + +DiagScalar( + every=globalEvery, + vars=['Uelm'] +) + +DiagFields( + every = 100, + fields = ['Ez']) diff --git a/doc/Sphinx/_static/Coordinate_Reference_AMcylindrical.png b/doc/Sphinx/_static/Coordinate_Reference_AMcylindrical.png index fb158d8fd..b71cb5571 100644 Binary files a/doc/Sphinx/_static/Coordinate_Reference_AMcylindrical.png and b/doc/Sphinx/_static/Coordinate_Reference_AMcylindrical.png differ diff --git a/doc/Sphinx/namelist.rst b/doc/Sphinx/namelist.rst index f61342026..4dd58f659 100755 --- a/doc/Sphinx/namelist.rst +++ b/doc/Sphinx/namelist.rst @@ -132,7 +132,8 @@ The block ``Main`` is **mandatory** and has the following syntax:: Boundary conditions must be set to ``"remove"`` for particles, ``"silver-muller"`` for longitudinal EM boundaries and ``"buneman"`` for transverse EM boundaries. - Vectorization, collisions, scalar diagnostics and + You can alternatively use ``"PML"`` for all boundaries. + Vectorization, collisions and order-4 interpolation are not supported yet. .. py:data:: interpolation_order @@ -272,7 +273,7 @@ The block ``Main`` is **mandatory** and has the following syntax:: :default: ``[["periodic"]]`` The boundary conditions for the electromagnetic fields. Each boundary may have one of - the following conditions: ``"periodic"``, ``"silver-muller"``, ``"reflective"`` or ``"ramp??"``. + the following conditions: ``"periodic"``, ``"silver-muller"``, ``"reflective"``, ``"ramp??"`` or ``"PML"``. | **Syntax 1:** ``[[bc_all]]``, identical for all boundaries. | **Syntax 2:** ``[[bc_X], [bc_Y], ...]``, different depending on x, y or z. @@ -294,6 +295,14 @@ The block ``Main`` is **mandatory** and has the following syntax:: Over the first half, the fields remain untouched. Over the second half, all fields are progressively reduced down to zero. + * ``"PML"`` stands for Perfectly Matched Layer. It is an open boundary condition. + The number of cells in the layer is defined by ``"number_of_pml_cells"``. + It supports laser injection as in ``"silver-muller"``. + + .. warning:: + + In the current release, in order to use PML all ``"EM_boundary_conditions"`` of the simulation must be ``"PML"``. + .. py:data:: EM_boundary_conditions_k :type: list of lists of floats @@ -311,6 +320,13 @@ The block ``Main`` is **mandatory** and has the following syntax:: | **Syntax 1:** ``[[1,0,0]]``, identical for all boundaries. | **Syntax 2:** ``[[1,0,0],[-1,0,0], ...]``, different on each boundary. +.. py:data:: number_of_pml_cells + + :type: list of lists of integer + :default: ``[[6,6],[6,6],[6,6]]`` + + Defines the number of cells in the ``"PML"`` layers using the same alternative syntaxes as ``"EM_boundary_conditions"``. + .. py:data:: time_fields_frozen :default: 0. @@ -2272,7 +2288,7 @@ This is done by including the block ``DiagScalar``:: .. warning:: - Scalars diagnostics are not yet supported in ``"AMcylindrical"`` geometry. + Scalars diagnostics min/max cell are not yet supported in ``"AMcylindrical"`` geometry. The full list of available scalars is given in the table below. diff --git a/doc/Sphinx/post-processing.rst b/doc/Sphinx/post-processing.rst index 7a396cf79..8568cd4da 100755 --- a/doc/Sphinx/post-processing.rst +++ b/doc/Sphinx/post-processing.rst @@ -419,9 +419,9 @@ and only one mode between those three. * ``hindex`` : the starting index of each proc in the hilbert curve * ``number_of_cells`` : the number of cells in each proc - * ``number_of_particles`` : the number of particles in each proc (except frozen ones) + * ``number_of_particles`` : the total number of non-frozen macro-particles in each proc (includes all species) * ``number_of_frozen_particles`` : the number of frozen particles in each proc - * ``total_load`` : the `load` of each proc (number of particles and cells with cell_load coefficient) + * ``total_load`` : the `load` of each proc (number of macro-particles and cells weighted by cell_load coefficients) * ``timer_global`` : global simulation time (only available for proc 0) * ``timer_particles`` : time spent computing particles by each proc * ``timer_maxwell`` : time spent solving maxwell by each proc diff --git a/doc/Sphinx/releases.rst b/doc/Sphinx/releases.rst index 74f1fc52f..1ea516637 100755 --- a/doc/Sphinx/releases.rst +++ b/doc/Sphinx/releases.rst @@ -20,6 +20,7 @@ Get Smilei Changes made in the repository (not released) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +* ``Perfectly Matched Layers`` boundary conditions for EM fields are available. Documentation is updated and a 2D Cartesian benchmark is provided. * Improved performance for ARM-based processors including the Fujitsu A64FX * Improved performance for GNU, LLVM, arm-clang and Fujitsu compilers on all types of architectures * Lasers can be injected from all boundaries @@ -75,7 +76,7 @@ Projects * More spectral solvers * GPU support - * Perfectly-matched layers + * Perfectly-matched layers for the envelope model ---- diff --git a/src/ElectroMagn/ElectroMagn.cpp b/src/ElectroMagn/ElectroMagn.cpp index e9f75c6f0..82f713d12 100755 --- a/src/ElectroMagn/ElectroMagn.cpp +++ b/src/ElectroMagn/ElectroMagn.cpp @@ -444,6 +444,8 @@ void ElectroMagn::boundaryConditions( int itime, double time_dual, Patch *patch, if( emBoundCond.size()>2 ) { if( emBoundCond[2]!=NULL ) { // <=> if !periodic emBoundCond[2]->apply( this, time_dual, patch ); + } + if( emBoundCond[3]!=NULL ) { // <=> if !periodic emBoundCond[3]->apply( this, time_dual, patch ); } } diff --git a/src/ElectroMagnBC/ElectroMagnBC.cpp b/src/ElectroMagnBC/ElectroMagnBC.cpp index 712b76170..b4c9cb9a3 100755 --- a/src/ElectroMagnBC/ElectroMagnBC.cpp +++ b/src/ElectroMagnBC/ElectroMagnBC.cpp @@ -41,6 +41,8 @@ ElectroMagnBC::ElectroMagnBC( Params ¶ms, Patch *patch, unsigned int i_bound d[i] = params.cell_length[i]; dt_ov_d[i] = dt / d[i]; } + + pml_solver_ = NULL; } // Destructor for ElectromagnBC diff --git a/src/ElectroMagnBC/ElectroMagnBC.h b/src/ElectroMagnBC/ElectroMagnBC.h index 887ebab11..25c5f0dd7 100755 --- a/src/ElectroMagnBC/ElectroMagnBC.h +++ b/src/ElectroMagnBC/ElectroMagnBC.h @@ -11,6 +11,7 @@ class Patch; class ElectroMagn; class Laser; class Field; +class Solver; class ElectroMagnBC { @@ -28,6 +29,19 @@ class ElectroMagnBC //! Vector for the various lasers std::vector vecLaser; + + virtual Field* getExPML() { ERROR("Not using PML");return NULL;} + virtual Field* getEyPML() { ERROR("Not using PML");return NULL;} + virtual Field* getEzPML() { ERROR("Not using PML");return NULL;} + virtual Field* getBxPML() { ERROR("Not using PML");return NULL;} + virtual Field* getByPML() { ERROR("Not using PML");return NULL;} + virtual Field* getBzPML() { ERROR("Not using PML");return NULL;} + virtual Field* getDxPML() { ERROR("Not using PML");return NULL;} + virtual Field* getDyPML() { ERROR("Not using PML");return NULL;} + virtual Field* getDzPML() { ERROR("Not using PML");return NULL;} + virtual Field* getHxPML() { ERROR("Not using PML");return NULL;} + virtual Field* getHyPML() { ERROR("Not using PML");return NULL;} + virtual Field* getHzPML() { ERROR("Not using PML");return NULL;} protected: @@ -48,6 +62,7 @@ class ElectroMagnBC //! Ratio of the time-step by the spatial-step std::vector dt_ov_d; + Solver* pml_solver_; }; #endif diff --git a/src/ElectroMagnBC/ElectroMagnBC2D_PML.cpp b/src/ElectroMagnBC/ElectroMagnBC2D_PML.cpp new file mode 100755 index 000000000..e24e691e7 --- /dev/null +++ b/src/ElectroMagnBC/ElectroMagnBC2D_PML.cpp @@ -0,0 +1,957 @@ +#include "ElectroMagnBC2D_PML.h" + +#include + +#include +#include + +#include "Params.h" +#include "Patch.h" +#include "ElectroMagn.h" +#include "Field2D.h" +#include "Tools.h" +#include "Laser.h" + +#include "SolverFactory.h" + +using namespace std; + +ElectroMagnBC2D_PML::ElectroMagnBC2D_PML( Params ¶ms, Patch *patch, unsigned int i_boundary ) + : ElectroMagnBC2D( params, patch, i_boundary ) +{ + + std::vector n_space(params.n_space); + std::vector oversize(params.oversize); + if( params.multiple_decomposition ) { + n_space = params.n_space_region; + oversize = params.region_oversize; + } + + pml_solver_ = SolverFactory::createPML( params ); + if (params.maxwell_sol=="Yee"){ + nsolver=2; + //MESSAGE("FDTD scheme in PML region : Yee."); + } + else if (params.maxwell_sol=="Bouchard"){ + nsolver=4; + //MESSAGE("FDTD scheme in PML region : Bouchard."); + } + else { + WARNING("The solver you use in the main domain is not the same as in the PML region. FDTD scheme in PML region : Yee."); + nsolver=2; + } + + if ( ( i_boundary_ == 0 && patch->isXmin() ) + || ( i_boundary_ == 1 && patch->isXmax() ) + || ( i_boundary_ == 2 && patch->isYmin() ) + || ( i_boundary_ == 3 && patch->isYmax() ) ) { + + int iDim = 0*((i_boundary_==0)||(i_boundary_==1))+1*((i_boundary_==2)||(i_boundary_==3)); + int min_or_max = (i_boundary_)%2; + + domain_oversize_x = oversize[0] ; + domain_oversize_y = oversize[1] ; + + if (patch->isXmin() ) {//&& i_boundary_ == 0 ) { + ncells_pml_xmin = params.number_of_pml_cells[0][0]; + ncells_pml_domain_xmin = ncells_pml_xmin + 1*oversize[0] + nsolver/2; + domain_oversize_x = oversize[0] ; + } + else { + ncells_pml_xmin = 0; + ncells_pml_domain_xmin = 0; + } + if (patch->isXmax() ) {//&& i_boundary_ == 1 ) { + ncells_pml_xmax = params.number_of_pml_cells[0][1]; + ncells_pml_domain_xmax = ncells_pml_xmax + 1*oversize[0] + nsolver/2; + domain_oversize_x = oversize[0] ; + } + else { + ncells_pml_xmax = 0; + ncells_pml_domain_xmax = 0; + } + if (patch->isYmin() ) {//&& i_boundary_ == 2 ) { + ncells_pml_ymin = params.number_of_pml_cells[1][0]; + ncells_pml_domain_ymin = ncells_pml_ymin + 1*oversize[1] + nsolver/2; + domain_oversize_y = oversize[1] ; + } + else { + ncells_pml_ymin = 0; + ncells_pml_domain_ymin = 0; + } + if (patch->isYmax() ) {//&& i_boundary_ == 3 ) { + ncells_pml_ymax = params.number_of_pml_cells[1][1]; + ncells_pml_domain_ymax = ncells_pml_ymax + 1*oversize[1] + nsolver/2; + domain_oversize_y = oversize[1] ; + } + else { + ncells_pml_ymax = 0; + ncells_pml_domain_ymax = 0; + } + + ncells_pml = params.number_of_pml_cells[iDim][min_or_max]; + ncells_pml_domain = ncells_pml+1*oversize[iDim] + nsolver/2; + + // Define min and max idx to exchange + // the good data f(solver,oversize) + if (min_or_max==0){ + // if min border : Exchange of data (for domain to pml-domain) + // min2exchange <= i < max2exchange + min2exchange = 1*nsolver/2 ; + max2exchange = 2*nsolver/2 ; + // Solver + solvermin = nsolver/2 ; + solvermax = ncells_pml_domain - oversize[iDim] ; + } + else if (min_or_max==1){ + // if max border : Exchange of data (for domain to pml-domain) + // min2exchange <= i < max2exchange + min2exchange = 1*nsolver/2 ; + max2exchange = 2*nsolver/2 ; + // Solver + solvermin = oversize[iDim] + nsolver/2 - nsolver/2 + 1 ; + solvermax = ncells_pml_domain-nsolver/2 ; + } + + if (ncells_pml==0){ + ERROR("PML domain have to be >0 cells in thickness"); + } + + std::vector dimPrim( params.nDim_field ); + for( int i=0 ; isetDomainSizeAndCoefficients( iDim, min_or_max, ncells_pml_domain, startpml, ncells_pml_min, ncells_pml_max, patch ); + + Ex_ = new Field2D( dimPrim, 0, false, "Ex_pml" ); + Ey_ = new Field2D( dimPrim, 1, false, "Ey_pml" ); + Ez_ = new Field2D( dimPrim, 2, false, "Ez_pml" ); + Bx_ = new Field2D( dimPrim, 0, true, "Bx_pml" ); + By_ = new Field2D( dimPrim, 1, true, "By_pml" ); + Bz_ = new Field2D( dimPrim, 2, true, "Bz_pml" ); + Dx_ = new Field2D( dimPrim, 0, false, "Dx_pml" ); + Dy_ = new Field2D( dimPrim, 1, false, "Dy_pml" ); + Dz_ = new Field2D( dimPrim, 2, false, "Dz_pml" ); + Hx_ = new Field2D( dimPrim, 0, true, "Hx_pml" ); + Hy_ = new Field2D( dimPrim, 1, true, "Hy_pml" ); + Hz_ = new Field2D( dimPrim, 2, true, "Hz_pml" ); + + //Laser parameter + double pyKx, pyKy; //, pyKz; + double kx, ky; //, kz; + double Knorm; + double omega = 1. ; + + factor_laser_space_time = 2.*dt_ov_d[0] ; + + // Xmin boundary + pyKx = params.EM_BCs_k[0][0]; + pyKy = params.EM_BCs_k[0][1]; + Knorm = sqrt( pyKx*pyKx + pyKy*pyKy ) ; + kx = omega*pyKx/Knorm; + ky = omega*pyKy/Knorm; + + factor_laser_angle_W = kx/Knorm; + + // Xmax boundary + pyKx = params.EM_BCs_k[1][0]; + pyKy = params.EM_BCs_k[1][1]; + Knorm = sqrt( pyKx*pyKx + pyKy*pyKy ) ; + kx = omega*pyKx/Knorm; + ky = omega*pyKy/Knorm; + + factor_laser_angle_E = kx/Knorm; + + // Ymin boundary + pyKx = params.EM_BCs_k[2][0]; + pyKy = params.EM_BCs_k[2][1]; + Knorm = sqrt( pyKx*pyKx + pyKy*pyKy ) ; + kx = omega*pyKx/Knorm; + ky = omega*pyKy/Knorm; + + factor_laser_angle_S = ky/Knorm; + + // Ymax boundary + pyKx = params.EM_BCs_k[3][0]; + pyKy = params.EM_BCs_k[3][1]; + Knorm = sqrt( pyKx*pyKx + pyKy*pyKy ) ; + kx = omega*pyKx/Knorm; + ky = omega*pyKy/Knorm; + + factor_laser_angle_N = ky/Knorm; + } +} + + +ElectroMagnBC2D_PML::~ElectroMagnBC2D_PML() +{ + delete Ex_; + delete Dx_; + delete Hx_; + delete Bx_; + delete Ey_; + delete Dy_; + delete Hy_; + delete By_; + delete Ez_; + delete Dz_; + delete Hz_; + delete Bz_; + + if (pml_solver_!=NULL) { + delete pml_solver_; + } +} + + +void ElectroMagnBC2D_PML::save_fields( Field *my_field, Patch *patch ) +{ +} + + +void ElectroMagnBC2D_PML::disableExternalFields() +{ +} + +// --------------------------------------------------------------------------------------------------------------------- +// Apply Boundary Conditions +// --------------------------------------------------------------------------------------------------------------------- + +void ElectroMagnBC2D_PML::apply( ElectroMagn *EMfields, double time_dual, Patch *patch ) +{ + int iDim = 0*((i_boundary_==0)||(i_boundary_==1))+1*((i_boundary_==2)||(i_boundary_==3)); + int min_or_max = (i_boundary_)%2; + + Field2D *Ex_domain = static_cast( EMfields->Ex_ ); + Field2D *Ey_domain = static_cast( EMfields->Ey_ ); + Field2D *Ez_domain = static_cast( EMfields->Ez_ ); + Field2D *Bx_domain = static_cast( EMfields->Bx_ ); + Field2D *By_domain = static_cast( EMfields->By_ ); + Field2D *Bz_domain = static_cast( EMfields->Bz_ ); + + if( i_boundary_ == 0 && patch->isXmin() ) { + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + + // 2. Exchange field PML <- Domain + for ( int i=min2exchange ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + //Injecting a laser + vector yp( 1 ); + for( unsigned int j=patch->isYmin() ; jisYmax() ; j++ ) { + + double byW = 0.; + yp[0] = patch->getDomainLocalMin( 1 ) + ( ( int )j - ( int )EMfields->oversize[1] )*d[1]; + + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + byW += vecLaser[ilaser]->getAmplitude0( yp, time_dual, j, 0 ); + } + (*Hy_)( ncells_pml_domain-domain_oversize_x-nsolver/2, j ) += factor_laser_space_time*factor_laser_angle_W*byW ; + (*By_)( ncells_pml_domain-domain_oversize_x-nsolver/2, j ) += factor_laser_space_time*factor_laser_angle_W*byW ; + } + + vector yd( 1 ); + for( unsigned int j=patch->isYmin() ; jisYmax() ; j++ ) { + + double bzW = 0.; + yd[0] = patch->getDomainLocalMin( 1 ) + ( ( int )j - 0.5 - ( int )EMfields->oversize[1] )*d[1]; + + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bzW += vecLaser[ilaser]->getAmplitude1( yd, time_dual, j, 0 ); + } + (*Hz_)( ncells_pml_domain-domain_oversize_x-nsolver/2, j ) += factor_laser_space_time*factor_laser_angle_W*bzW ; + (*Bz_)( ncells_pml_domain-domain_oversize_x-nsolver/2, j ) += factor_laser_space_time*factor_laser_angle_W*bzW ; + } + + // 4. Exchange PML -> Domain + // Primals in x-direction + for (int i=0 ; i < nsolver/2 ; i++){ + for ( int j=0 ; jisXmax() ) { + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + // 2. Exchange field Domain -> PML + for ( int i=min2exchange ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + //Injecting a laser + vector yp( 1 ); + for( unsigned int j=patch->isYmin() ; jisYmax() ; j++ ) { + + double byE = 0.; + yp[0] = patch->getDomainLocalMin( 1 ) + ( ( int )j - ( int )EMfields->oversize[1] )*d[1]; + + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + byE += vecLaser[ilaser]->getAmplitude0( yp, time_dual, j, 0 ); + } + (*Hy_)( domain_oversize_x+nsolver/2, j ) += factor_laser_space_time*factor_laser_angle_E*byE ; + (*By_)( domain_oversize_x+nsolver/2, j ) += factor_laser_space_time*factor_laser_angle_E*byE ; + } + + vector yd( 1 ); + for( unsigned int j=patch->isYmin() ; jisYmax() ; j++ ) { + + double bzE = 0.; + yd[0] = patch->getDomainLocalMin( 1 ) + ( ( int )j - 0.5 - ( int )EMfields->oversize[1] )*d[1]; + + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bzE += vecLaser[ilaser]->getAmplitude1( yd, time_dual, j, 0 ); + } + (*Hz_)( domain_oversize_x+nsolver/2, j ) += factor_laser_space_time*factor_laser_angle_E*bzE ; + (*Bz_)( domain_oversize_x+nsolver/2, j ) += factor_laser_space_time*factor_laser_angle_E*bzE ; + } + + // 4. Exchange Domain -> PML + // Primals in x-direction + for (int i=0 ; i < nsolver/2-1 ; i++){ + for ( int j=0 ; jisYmin() ) { + + ElectroMagnBC2D_PML* pml_fields_xmin = NULL; + ElectroMagnBC2D_PML* pml_fields_xmax = NULL; + + if(ncells_pml_xmin != 0){ + pml_fields_xmin = static_cast( EMfields->emBoundCond[0] ); + } + if(ncells_pml_xmax != 0){ + pml_fields_xmax = static_cast( EMfields->emBoundCond[1] ); + } + + Field2D* Ex_pml_xmin = NULL; + Field2D* Ey_pml_xmin = NULL; + Field2D* Ez_pml_xmin = NULL; + Field2D* Hx_pml_xmin = NULL; + Field2D* Hy_pml_xmin = NULL; + Field2D* Hz_pml_xmin = NULL; + Field2D* Dx_pml_xmin = NULL; + Field2D* Dy_pml_xmin = NULL; + Field2D* Dz_pml_xmin = NULL; + Field2D* Bx_pml_xmin = NULL; + Field2D* By_pml_xmin = NULL; + Field2D* Bz_pml_xmin = NULL; + + Field2D* Ex_pml_xmax = NULL; + Field2D* Ey_pml_xmax = NULL; + Field2D* Ez_pml_xmax = NULL; + Field2D* Hx_pml_xmax = NULL; + Field2D* Hy_pml_xmax = NULL; + Field2D* Hz_pml_xmax = NULL; + Field2D* Dx_pml_xmax = NULL; + Field2D* Dy_pml_xmax = NULL; + Field2D* Dz_pml_xmax = NULL; + Field2D* Bx_pml_xmax = NULL; + Field2D* By_pml_xmax = NULL; + Field2D* Bz_pml_xmax = NULL; + + if(ncells_pml_xmin != 0){ + Ex_pml_xmin = pml_fields_xmin->Ex_; + Ey_pml_xmin = pml_fields_xmin->Ey_; + Ez_pml_xmin = pml_fields_xmin->Ez_; + Hx_pml_xmin = pml_fields_xmin->Hx_; + Hy_pml_xmin = pml_fields_xmin->Hy_; + Hz_pml_xmin = pml_fields_xmin->Hz_; + Dx_pml_xmin = pml_fields_xmin->Dx_; + Dy_pml_xmin = pml_fields_xmin->Dy_; + Dz_pml_xmin = pml_fields_xmin->Dz_; + Bx_pml_xmin = pml_fields_xmin->Bx_; + By_pml_xmin = pml_fields_xmin->By_; + Bz_pml_xmin = pml_fields_xmin->Bz_; + } + + if(ncells_pml_xmax != 0){ + Ex_pml_xmax = pml_fields_xmax->Ex_; + Ey_pml_xmax = pml_fields_xmax->Ey_; + Ez_pml_xmax = pml_fields_xmax->Ez_; + Hx_pml_xmax = pml_fields_xmax->Hx_; + Hy_pml_xmax = pml_fields_xmax->Hy_; + Hz_pml_xmax = pml_fields_xmax->Hz_; + Dx_pml_xmax = pml_fields_xmax->Dx_; + Dy_pml_xmax = pml_fields_xmax->Dy_; + Dz_pml_xmax = pml_fields_xmax->Dz_; + Bx_pml_xmax = pml_fields_xmax->Bx_; + By_pml_xmax = pml_fields_xmax->By_; + Bz_pml_xmax = pml_fields_xmax->Bz_; + } + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + // // 2. Exchange field PML <- Domain + for ( int j=min2exchange ; jisXmin()) { + if(ncells_pml_xmin != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + //Injecting a laser + vector xp( 1 ); + for( unsigned int i=patch->isXmin() ; iisXmax() ; i++ ) { + + double bxS = 0.; + xp[0] = patch->getDomainLocalMin( 0 ) + ( ( int )i - ( int )EMfields->oversize[0] )*d[0]; + + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bxS += vecLaser[ilaser]->getAmplitude0( xp, time_dual, i, 0 ); + } + (*Hx_)( ncells_pml_xmin+i, ncells_pml_domain-domain_oversize_y-nsolver/2 ) += factor_laser_space_time*factor_laser_angle_S*bxS ; + (*Bx_)( ncells_pml_xmin+i, ncells_pml_domain-domain_oversize_y-nsolver/2 ) += factor_laser_space_time*factor_laser_angle_S*bxS ; + } + + vector xd( 1 ); + for( unsigned int i=patch->isXmin() ; iisXmax() ; i++ ) { + + double bzS = 0.; + xd[0] = patch->getDomainLocalMin( 0 ) + ( ( int )i - 0.5 - ( int )EMfields->oversize[0] )*d[0]; + + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bzS += vecLaser[ilaser]->getAmplitude1( xd, time_dual, i, 0 ); + } + (*Hz_)( ncells_pml_xmin+i, ncells_pml_domain-domain_oversize_y-nsolver/2 ) += factor_laser_space_time*factor_laser_angle_S*bzS ; + (*Bz_)( ncells_pml_xmin+i, ncells_pml_domain-domain_oversize_y-nsolver/2 ) += factor_laser_space_time*factor_laser_angle_S*bzS ; + } + + // 4. Exchange PML -> Domain + // Primals in y-direction + for (int j=0 ; j < nsolver/2 ; j++){ + for ( int i=0 ; i PML x MIN + // Primal in y-direction + for (int j=0 ; j < nsolver/2 ; j++){ + if (patch->isXmin()) { + if(ncells_pml_xmin != 0){ + for ( int i=0 ; iisXmin()) { + if(ncells_pml_xmin != 0){ + for ( int i=0 ; i PML x MAX + // Primal in y-direction + for (int j=0 ; j < nsolver/2 ; j++){ + if (patch->isXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; iisYmax() ) { + + ElectroMagnBC2D_PML* pml_fields_xmin = NULL ; + ElectroMagnBC2D_PML* pml_fields_xmax = NULL ; + + if(ncells_pml_xmin != 0){ + pml_fields_xmin = static_cast( EMfields->emBoundCond[0] ); + } + if(ncells_pml_xmax != 0){ + pml_fields_xmax = static_cast( EMfields->emBoundCond[1] ); + } + + Field2D* Ex_pml_xmin = NULL; + Field2D* Ey_pml_xmin = NULL; + Field2D* Ez_pml_xmin = NULL; + Field2D* Hx_pml_xmin = NULL; + Field2D* Hy_pml_xmin = NULL; + Field2D* Hz_pml_xmin = NULL; + Field2D* Dx_pml_xmin = NULL; + Field2D* Dy_pml_xmin = NULL; + Field2D* Dz_pml_xmin = NULL; + Field2D* Bx_pml_xmin = NULL; + Field2D* By_pml_xmin = NULL; + Field2D* Bz_pml_xmin = NULL; + + Field2D* Ex_pml_xmax = NULL; + Field2D* Ey_pml_xmax = NULL; + Field2D* Ez_pml_xmax = NULL; + Field2D* Hx_pml_xmax = NULL; + Field2D* Hy_pml_xmax = NULL; + Field2D* Hz_pml_xmax = NULL; + Field2D* Dx_pml_xmax = NULL; + Field2D* Dy_pml_xmax = NULL; + Field2D* Dz_pml_xmax = NULL; + Field2D* Bx_pml_xmax = NULL; + Field2D* By_pml_xmax = NULL; + Field2D* Bz_pml_xmax = NULL; + + if(ncells_pml_xmin != 0){ + Ex_pml_xmin = pml_fields_xmin->Ex_; + Ey_pml_xmin = pml_fields_xmin->Ey_; + Ez_pml_xmin = pml_fields_xmin->Ez_; + Hx_pml_xmin = pml_fields_xmin->Hx_; + Hy_pml_xmin = pml_fields_xmin->Hy_; + Hz_pml_xmin = pml_fields_xmin->Hz_; + Dx_pml_xmin = pml_fields_xmin->Dx_; + Dy_pml_xmin = pml_fields_xmin->Dy_; + Dz_pml_xmin = pml_fields_xmin->Dz_; + Bx_pml_xmin = pml_fields_xmin->Bx_; + By_pml_xmin = pml_fields_xmin->By_; + Bz_pml_xmin = pml_fields_xmin->Bz_; + } + + if(ncells_pml_xmax != 0){ + Ex_pml_xmax = pml_fields_xmax->Ex_; + Ey_pml_xmax = pml_fields_xmax->Ey_; + Ez_pml_xmax = pml_fields_xmax->Ez_; + Hx_pml_xmax = pml_fields_xmax->Hx_; + Hy_pml_xmax = pml_fields_xmax->Hy_; + Hz_pml_xmax = pml_fields_xmax->Hz_; + Dx_pml_xmax = pml_fields_xmax->Dx_; + Dy_pml_xmax = pml_fields_xmax->Dy_; + Dz_pml_xmax = pml_fields_xmax->Dz_; + Bx_pml_xmax = pml_fields_xmax->Bx_; + By_pml_xmax = pml_fields_xmax->By_; + Bz_pml_xmax = pml_fields_xmax->Bz_; + } + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + // 2. Exchange field PML <- Domain + for ( int j=min2exchange ; jisXmin()) { + if(ncells_pml_xmin != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + //Injecting a laser + vector xp( 1 ); + for( unsigned int i=patch->isXmin() ; iisXmax() ; i++ ) { + + double bxN = 0.; + xp[0] = patch->getDomainLocalMin( 0 ) + ( ( int )i - ( int )EMfields->oversize[0] )*d[0]; + + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bxN += vecLaser[ilaser]->getAmplitude0( xp, time_dual, i, 0 ); + } + (*Hx_)( ncells_pml_xmin+i , domain_oversize_y+nsolver/2 ) += factor_laser_space_time*factor_laser_angle_N*bxN ; + (*Bx_)( ncells_pml_xmin+i , domain_oversize_y+nsolver/2 ) += factor_laser_space_time*factor_laser_angle_N*bxN ; + } + + vector xd( 1 ); + for( unsigned int i=patch->isXmin() ; iisXmax() ; i++ ) { + + double bzN = 0.; + xd[0] = patch->getDomainLocalMin( 0 ) + ( ( int )i - 0.5 - ( int )EMfields->oversize[0] )*d[0]; + + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bzN += vecLaser[ilaser]->getAmplitude1( xd, time_dual, i, 0 ); + } + (*Hz_)( ncells_pml_xmin+i, domain_oversize_y+nsolver/2 ) += factor_laser_space_time*factor_laser_angle_N*bzN ; + (*Bz_)( ncells_pml_xmin+i, domain_oversize_y+nsolver/2 ) += factor_laser_space_time*factor_laser_angle_N*bzN ; + } + + // 4. Exchange PML -> Domain + // Primals in y-direction + for (int j=0 ; j < nsolver/2-1 ; j++){ + for ( int i=0 ; i PML x MIN + // Primals in y-direction + for (int j=0 ; j < nsolver/2-1 ; j++){ + if (patch->isXmin()) { + if(ncells_pml_xmin != 0){ + for ( int i=0 ; iisXmin()) { + if(ncells_pml_xmin != 0){ + for ( int i=0 ; i PML x MAX + // Primals in y-direction + for (int j=0 ; j < nsolver/2-1 ; j++){ + if (patch->isXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; i +#include "Tools.h" +#include "ElectroMagnBC2D.h" +#include "ElectroMagn2D.h" +#include "Field2D.h" + +class Params; +class ElectroMagn; +class Field; +class PML_Solver2D_Yee; +class PML_Solver2D_Bouchard; + +class ElectroMagnBC2D_PML : public ElectroMagnBC2D +{ +public: + + ElectroMagnBC2D_PML( Params ¶ms, Patch *patch, unsigned int _min_max ); + ~ElectroMagnBC2D_PML(); + + virtual void apply( ElectroMagn *EMfields, double time_dual, Patch *patch ) override; + + void save_fields( Field *, Patch *patch ) override; + void disableExternalFields() override; + + Field2D* Ex_ = NULL; + Field2D* Ey_ = NULL; + Field2D* Ez_ = NULL; + Field2D* Bx_ = NULL; + Field2D* By_ = NULL; + Field2D* Bz_ = NULL; + Field2D* Dx_ = NULL; + Field2D* Dy_ = NULL; + Field2D* Dz_ = NULL; + Field2D* Hx_ = NULL; + Field2D* Hy_ = NULL; + Field2D* Hz_ = NULL; + + Field* getExPML() override { return Ex_; }; + Field* getEyPML() override { return Ey_; }; + Field* getEzPML() override { return Ez_; }; + Field* getBxPML() override { return Bx_; }; + Field* getByPML() override { return By_; }; + Field* getBzPML() override { return Bz_; }; + Field* getDxPML() override { return Dx_; }; + Field* getDyPML() override { return Dy_; }; + Field* getDzPML() override { return Dz_; }; + Field* getHxPML() override { return Hx_; }; + Field* getHyPML() override { return Hy_; }; + Field* getHzPML() override { return Hz_; }; + + int domain_oversize_x; + int domain_oversize_y; + + int ncells_pml_xmin ; + int ncells_pml_xmax ; + int ncells_pml_ymin ; + int ncells_pml_ymax ; + + int ncells_pml_domain_xmin; + int ncells_pml_domain_xmax; + int ncells_pml_domain_ymin; + int ncells_pml_domain_ymax; + + int ncells_pml; + int ncells_pml_domain; + + int nsolver; + + int min2exchange; + int max2exchange; + int solvermin; + int solvermax; + + int ypml_size_in_x; + int startpml; + + double factor_laser_space_time ; + double factor_laser_angle_W ; + double factor_laser_angle_E ; + double factor_laser_angle_S ; + double factor_laser_angle_N ; + +}; + +#endif + diff --git a/src/ElectroMagnBC/ElectroMagnBC3D_PML.cpp b/src/ElectroMagnBC/ElectroMagnBC3D_PML.cpp new file mode 100644 index 000000000..036d1a36a --- /dev/null +++ b/src/ElectroMagnBC/ElectroMagnBC3D_PML.cpp @@ -0,0 +1,2384 @@ +#include "ElectroMagnBC3D_PML.h" + +#include + +#include +#include + +#include "Params.h" +#include "Patch.h" +#include "ElectroMagn.h" +#include "Field3D.h" +#include "Tools.h" +#include "Laser.h" + +#include "SolverFactory.h" + +using namespace std; + +ElectroMagnBC3D_PML::ElectroMagnBC3D_PML( Params ¶ms, Patch *patch, unsigned int i_boundary ) + : ElectroMagnBC3D( params, patch, i_boundary ) +{ + + std::vector n_space(params.n_space); + std::vector oversize(params.oversize); + if( params.multiple_decomposition ) { + n_space = params.n_space_region; + oversize = params.region_oversize; + } + + pml_solver_ = SolverFactory::createPML( params ); + if (params.maxwell_sol=="Yee"){ + nsolver=2; + //MESSAGE("FDTD scheme in PML region : Yee."); + } + else if (params.maxwell_sol=="Bouchard"){ + nsolver=4; + //MESSAGE("FDTD scheme in PML region : Bouchard."); + } + else { + WARNING("The solver you use in the main domain is not the same as in the PML region. FDTD scheme in PML region : Yee."); + nsolver=2; + } + + if ( ( i_boundary_ == 0 && patch->isXmin() ) + || ( i_boundary_ == 1 && patch->isXmax() ) + || ( i_boundary_ == 2 && patch->isYmin() ) + || ( i_boundary_ == 3 && patch->isYmax() ) + || ( i_boundary_ == 4 && patch->isZmin() ) + || ( i_boundary_ == 5 && patch->isZmax() ) ) { + + int iDim = 0*((i_boundary_==0)||(i_boundary_==1))+1*((i_boundary_==2)||(i_boundary_==3))+2*((i_boundary_==4)||(i_boundary_==5)); + int min_or_max = (i_boundary_)%2; + + domain_oversize_x = oversize[0] ; + domain_oversize_y = oversize[1] ; + domain_oversize_z = oversize[2] ; + + if (patch->isXmin() ) {//&& i_boundary_ == 0 ) { + ncells_pml_xmin = params.number_of_pml_cells[0][0]; + ncells_pml_domain_xmin = ncells_pml_xmin + 1*oversize[0] + nsolver/2; + domain_oversize_x = oversize[0] ; + } + else { + ncells_pml_xmin = 0; + ncells_pml_domain_xmin = 0; + } + if (patch->isXmax() ) {//&& i_boundary_ == 1 ) { + ncells_pml_xmax = params.number_of_pml_cells[0][1]; + ncells_pml_domain_xmax = ncells_pml_xmax + 1*oversize[0] + nsolver/2; + domain_oversize_x = oversize[0] ; + } + else { + ncells_pml_xmax = 0; + ncells_pml_domain_xmax = 0; + } + if (patch->isYmin() ) {//&& i_boundary_ == 2 ) { + ncells_pml_ymin = params.number_of_pml_cells[1][0]; + ncells_pml_domain_ymin = ncells_pml_ymin + 1*oversize[1] + nsolver/2; + domain_oversize_y = oversize[1] ; + } + else { + ncells_pml_ymin = 0; + ncells_pml_domain_ymin = 0; + } + if (patch->isYmax() ) {//&& i_boundary_ == 3 ) { + ncells_pml_ymax = params.number_of_pml_cells[1][1]; + ncells_pml_domain_ymax = ncells_pml_ymax + 1*oversize[1] + nsolver/2; + domain_oversize_y = oversize[1] ; + } + else { + ncells_pml_ymax = 0; + ncells_pml_domain_ymax = 0; + } + if (patch->isZmin() ) {//&& i_boundary_ == 2 ) { + ncells_pml_zmin = params.number_of_pml_cells[2][0]; + ncells_pml_domain_zmin = ncells_pml_zmin + 1*oversize[2] + nsolver/2; + domain_oversize_z = oversize[2] ; + } + else { + ncells_pml_zmin = 0; + ncells_pml_domain_zmin = 0; + } + if (patch->isZmax() ) {//&& i_boundary_ == 3 ) { + ncells_pml_zmax = params.number_of_pml_cells[2][1]; + ncells_pml_domain_zmax = ncells_pml_zmax + 1*oversize[2] + nsolver/2; + domain_oversize_z = oversize[2] ; + } + else { + ncells_pml_zmax = 0; + ncells_pml_domain_zmax = 0; + } + + ncells_pml = params.number_of_pml_cells[iDim][min_or_max]; + ncells_pml_domain = ncells_pml+1*oversize[iDim] + nsolver/2; + + // Define min and max idx to exchange + // the good data f(solver,oversize) + if (min_or_max==0){ + // if min border : Exchange of data (for domain to pml-domain) + // min2exchange <= i < max2exchange + min2exchange = 1*nsolver/2 ; + max2exchange = 2*nsolver/2 ; + // Solver + solvermin = nsolver/2 ; + solvermax = ncells_pml_domain - oversize[iDim] ; + } + else if (min_or_max==1){ + // if max border : Exchange of data (for domain to pml-domain) + // min2exchange <= i < max2exchange + min2exchange = 1*nsolver/2 ; + max2exchange = 2*nsolver/2 ; + // Solver + solvermin = oversize[iDim] + nsolver/2 - nsolver/2 + 1 ; + solvermax = ncells_pml_domain-nsolver/2 ; + } + + if (ncells_pml==0){ + ERROR("PML domain have to be >0 cells in thickness"); + } + + std::vector dimPrim( params.nDim_field ); + for( int i=0 ; i ncells_pml + // +1 -> +nsolver/2 + // +2*oversize -> +1*oversize + if ( iDim==1 ){ + // If the PML domain in Y is in Xmin or Xmax too, add cell orthogonally + dimPrim[iDim-1] += ncells_pml_xmin + ncells_pml_xmax ; + ypml_size_in_x = dimPrim[iDim-1] ; + } + if ( iDim==2 ){ + dimPrim[0] += ncells_pml_xmin + ncells_pml_xmax ; + dimPrim[1] += ncells_pml_ymin + ncells_pml_ymax ; + zpml_size_in_x = dimPrim[0] ; + zpml_size_in_y = dimPrim[1] ; + } + + startpml = oversize[iDim]+nsolver/2; + + int ncells_pml_min[2]; + ncells_pml_min[0] = ncells_pml_xmin; + ncells_pml_min[1] = ncells_pml_ymin; + int ncells_pml_max[2]; + ncells_pml_max[0] = ncells_pml_xmax; + ncells_pml_max[1] = ncells_pml_ymax; + + pml_solver_->setDomainSizeAndCoefficients( iDim, min_or_max, ncells_pml_domain, startpml, ncells_pml_min, ncells_pml_max, patch ); + + Ex_ = new Field3D( dimPrim, 0, false, "Ex_pml" ); + Ey_ = new Field3D( dimPrim, 1, false, "Ey_pml" ); + Ez_ = new Field3D( dimPrim, 2, false, "Ez_pml" ); + Bx_ = new Field3D( dimPrim, 0, true, "Bx_pml" ); + By_ = new Field3D( dimPrim, 1, true, "By_pml" ); + Bz_ = new Field3D( dimPrim, 2, true, "Bz_pml" ); + Dx_ = new Field3D( dimPrim, 0, false, "Dx_pml" ); + Dy_ = new Field3D( dimPrim, 1, false, "Dy_pml" ); + Dz_ = new Field3D( dimPrim, 2, false, "Dz_pml" ); + Hx_ = new Field3D( dimPrim, 0, true, "Hx_pml" ); + Hy_ = new Field3D( dimPrim, 1, true, "Hy_pml" ); + Hz_ = new Field3D( dimPrim, 2, true, "Hz_pml" ); + + //Laser parameter + double pyKx, pyKy, pyKz; + double kx, ky, kz; + double Knorm; + double omega = 1. ; + + // Xmin boundary + pyKx = params.EM_BCs_k[0][0]; + pyKy = params.EM_BCs_k[0][1]; + pyKz = params.EM_BCs_k[0][2]; + Knorm = sqrt( pyKx*pyKx + pyKy*pyKy + pyKz*pyKz ) ; + kx = omega*pyKx/Knorm; + ky = omega*pyKy/Knorm; + kz = omega*pyKz/Knorm; + + factor_laser_space_time = 2.*dt_ov_d[0] ; + factor_laser_angle_W = factor_laser_space_time*kx/Knorm; + + // Xmax boundary + pyKx = params.EM_BCs_k[1][0]; + pyKy = params.EM_BCs_k[1][1]; + pyKz = params.EM_BCs_k[1][2]; + Knorm = sqrt( pyKx*pyKx + pyKy*pyKy + pyKz*pyKz ) ; + kx = omega*pyKx/Knorm; + ky = omega*pyKy/Knorm; + kz = omega*pyKz/Knorm; + + factor_laser_space_time = 2.*dt_ov_d[0] ; + factor_laser_angle_E = factor_laser_space_time*kx/Knorm; + + // Ymin boundary + pyKx = params.EM_BCs_k[2][0]; + pyKy = params.EM_BCs_k[2][1]; + pyKz = params.EM_BCs_k[2][2]; + Knorm = sqrt( pyKx*pyKx + pyKy*pyKy + pyKz*pyKz ) ; + kx = omega*pyKx/Knorm; + ky = omega*pyKy/Knorm; + kz = omega*pyKz/Knorm; + + factor_laser_space_time = 2.*dt_ov_d[1] ; + factor_laser_angle_S = factor_laser_space_time*ky/Knorm; + + // Ymax boundary + pyKx = params.EM_BCs_k[3][0]; + pyKy = params.EM_BCs_k[3][1]; + pyKz = params.EM_BCs_k[3][2]; + Knorm = sqrt( pyKx*pyKx + pyKy*pyKy + pyKz*pyKz ) ; + kx = omega*pyKx/Knorm; + ky = omega*pyKy/Knorm; + kz = omega*pyKz/Knorm; + + factor_laser_space_time = 2.*dt_ov_d[1] ; + factor_laser_angle_N = factor_laser_space_time*ky/Knorm; + + // Zmin boundary + pyKx = params.EM_BCs_k[4][0]; + pyKy = params.EM_BCs_k[4][1]; + pyKz = params.EM_BCs_k[4][2]; + Knorm = sqrt( pyKx*pyKx + pyKy*pyKy + pyKz*pyKz ) ; + kx = omega*pyKx/Knorm; + ky = omega*pyKy/Knorm; + kz = omega*pyKz/Knorm; + + factor_laser_space_time = 2.*dt_ov_d[2] ; + factor_laser_angle_B = factor_laser_space_time*kz/Knorm; + + // Zmax boundary + pyKx = params.EM_BCs_k[5][0]; + pyKy = params.EM_BCs_k[5][1]; + pyKz = params.EM_BCs_k[5][2]; + Knorm = sqrt( pyKx*pyKx + pyKy*pyKy + pyKz*pyKz ) ; + kx = omega*pyKx/Knorm; + ky = omega*pyKy/Knorm; + kz = omega*pyKz/Knorm; + + factor_laser_space_time = 2.*dt_ov_d[2] ; + factor_laser_angle_T = factor_laser_space_time*kz/Knorm; + + + } +} + + +ElectroMagnBC3D_PML::~ElectroMagnBC3D_PML() +{ + delete Ex_; + delete Dx_; + delete Hx_; + delete Bx_; + delete Ey_; + delete Dy_; + delete Hy_; + delete By_; + delete Ez_; + delete Dz_; + delete Hz_; + delete Bz_; + + if (pml_solver_!=NULL) { + delete pml_solver_; + } +} + + +void ElectroMagnBC3D_PML::save_fields( Field *my_field, Patch *patch ) +{ +} + + +void ElectroMagnBC3D_PML::disableExternalFields() +{ +} + + +// --------------------------------------------------------------------------------------------------------------------- +// Apply Boundary Conditions +// --------------------------------------------------------------------------------------------------------------------- +void ElectroMagnBC3D_PML::apply( ElectroMagn *EMfields, double time_dual, Patch *patch ) +{ + + int iDim = 0*((i_boundary_==0)||(i_boundary_==1))+1*((i_boundary_==2)||(i_boundary_==3))+2*((i_boundary_==4)||(i_boundary_==5)); + int min_or_max = (i_boundary_)%2; + + Field3D *Ex_domain = static_cast( EMfields->Ex_ ); + Field3D *Ey_domain = static_cast( EMfields->Ey_ ); + Field3D *Ez_domain = static_cast( EMfields->Ez_ ); + Field3D *Bx_domain = static_cast( EMfields->Bx_ ); + Field3D *By_domain = static_cast( EMfields->By_ ); + Field3D *Bz_domain = static_cast( EMfields->Bz_ ); + + vector pos( 2 ); + + if( i_boundary_ == 0 && patch->isXmin() ) { + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + // 2. Exchange field PML <- Domain + for ( int i=min2exchange ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + //Injecting a laser + vector by( n_p[1]*n_d[2], 0. ); + for( unsigned int j=patch->isYmin() ; jisYmax() ; j++ ) { + pos[0] = patch->getDomainLocalMin( 1 ) + ( ( int )j - ( int )EMfields->oversize[1] )*d[1]; + for( unsigned int k=patch->isZmin() ; kisZmax() ; k++ ) { + pos[1] = patch->getDomainLocalMin( 2 ) + ( ( int )k - 0.5 - ( int )EMfields->oversize[2] )*d[2]; + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + by[ j*n_d[2]+k ] += vecLaser[ilaser]->getAmplitude0( pos, time_dual, j, k ); + } + (*Hy_)( ncells_pml_domain-domain_oversize_x-nsolver/2, j , k ) += factor_laser_angle_W*by[j*n_d[2]+k] ; + (*By_)( ncells_pml_domain-domain_oversize_x-nsolver/2, j , k ) += factor_laser_angle_W*by[j*n_d[2]+k] ; + } + } + + vector bz( n_d[1]*n_p[2], 0. ); + for( unsigned int j=patch->isYmin() ; jisYmax() ; j++ ) { + pos[0] = patch->getDomainLocalMin( 1 ) + ( ( int )j - 0.5 - ( int )EMfields->oversize[1] )*d[1]; + for( unsigned int k=patch->isZmin() ; kisZmax() ; k++ ) { + pos[1] = patch->getDomainLocalMin( 2 ) + ( ( int )k - ( int )EMfields->oversize[2] )*d[2]; + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bz[ j*n_p[2]+k] += vecLaser[ilaser]->getAmplitude1( pos, time_dual, j, k ); + } + (*Hz_)( ncells_pml_domain-domain_oversize_x-nsolver/2, j , k) += factor_laser_angle_W*bz[ j*n_p[2]+k ] ; + (*Bz_)( ncells_pml_domain-domain_oversize_x-nsolver/2, j , k) += factor_laser_angle_W*bz[ j*n_p[2]+k ] ; + } + } + + // 4. Exchange PML -> Domain + // Primals in x-direction + for (int i=0 ; i < nsolver/2 ; i++){ + for ( int j=0 ; jisXmax() ) { + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + // 2. Exchange field Domain -> PML + for ( int i=min2exchange ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + //Injecting a laser + vector by( n_p[1]*n_d[2], 0. ); + for( unsigned int j=patch->isYmin() ; jisYmax() ; j++ ) { + pos[0] = patch->getDomainLocalMin( 1 ) + ( ( int )j - ( int )EMfields->oversize[1] )*d[1]; + for( unsigned int k=patch->isZmin() ; kisZmax() ; k++ ) { + pos[1] = patch->getDomainLocalMin( 2 ) + ( ( int )k - 0.5 - ( int )EMfields->oversize[2] )*d[2]; + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + by[ j*n_d[2]+k ] += vecLaser[ilaser]->getAmplitude0( pos, time_dual, j, k ); + } + (*Hy_)( domain_oversize_x+nsolver/2, j, k ) += factor_laser_angle_E*by[j*n_d[2]+k] ; + (*By_)( domain_oversize_x+nsolver/2, j, k ) += factor_laser_angle_E*by[j*n_d[2]+k] ; + } + } + + vector bz( n_d[1]*n_p[2], 0. ); + for( unsigned int j=patch->isYmin() ; jisYmax() ; j++ ) { + pos[0] = patch->getDomainLocalMin( 1 ) + ( ( int )j - 0.5 - ( int )EMfields->oversize[1] )*d[1]; + for( unsigned int k=patch->isZmin() ; kisZmax() ; k++ ) { + pos[1] = patch->getDomainLocalMin( 2 ) + ( ( int )k - ( int )EMfields->oversize[2] )*d[2]; + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bz[ j*n_p[2]+k] += vecLaser[ilaser]->getAmplitude1( pos, time_dual, j, k ); + } + (*Hz_)( domain_oversize_x+nsolver/2, j, k ) += factor_laser_angle_E*bz[ j*n_p[2]+k ] ; + (*Bz_)( domain_oversize_x+nsolver/2, j, k ) += factor_laser_angle_E*bz[ j*n_p[2]+k ] ; + } + } + + // 4. Exchange Domain -> PML + // Primals in x-direction + for (int i=0 ; i < nsolver/2-1 ; i++){ + for ( int j=0 ; jisYmin() ) { + + ElectroMagnBC3D_PML* pml_fields_xmin = NULL ; + ElectroMagnBC3D_PML* pml_fields_xmax = NULL ; + + if(ncells_pml_xmin != 0){ + pml_fields_xmin = static_cast( EMfields->emBoundCond[0] ); + } + if(ncells_pml_xmax != 0){ + pml_fields_xmax = static_cast( EMfields->emBoundCond[1] ); + } + + Field3D* Ex_pml_xmin = NULL; + Field3D* Ey_pml_xmin = NULL; + Field3D* Ez_pml_xmin = NULL; + Field3D* Hx_pml_xmin = NULL; + Field3D* Hy_pml_xmin = NULL; + Field3D* Hz_pml_xmin = NULL; + Field3D* Dx_pml_xmin = NULL; + Field3D* Dy_pml_xmin = NULL; + Field3D* Dz_pml_xmin = NULL; + Field3D* Bx_pml_xmin = NULL; + Field3D* By_pml_xmin = NULL; + + Field3D* Bz_pml_xmin = NULL; + Field3D* Ex_pml_xmax = NULL; + Field3D* Ey_pml_xmax = NULL; + Field3D* Ez_pml_xmax = NULL; + Field3D* Hx_pml_xmax = NULL; + Field3D* Hy_pml_xmax = NULL; + Field3D* Hz_pml_xmax = NULL; + Field3D* Dx_pml_xmax = NULL; + Field3D* Dy_pml_xmax = NULL; + Field3D* Dz_pml_xmax = NULL; + Field3D* Bx_pml_xmax = NULL; + Field3D* By_pml_xmax = NULL; + Field3D* Bz_pml_xmax = NULL; + + if(ncells_pml_xmin != 0){ + Ex_pml_xmin = pml_fields_xmin->Ex_; + Ey_pml_xmin = pml_fields_xmin->Ey_; + Ez_pml_xmin = pml_fields_xmin->Ez_; + Hx_pml_xmin = pml_fields_xmin->Hx_; + Hy_pml_xmin = pml_fields_xmin->Hy_; + Hz_pml_xmin = pml_fields_xmin->Hz_; + Dx_pml_xmin = pml_fields_xmin->Dx_; + Dy_pml_xmin = pml_fields_xmin->Dy_; + Dz_pml_xmin = pml_fields_xmin->Dz_; + Bx_pml_xmin = pml_fields_xmin->Bx_; + By_pml_xmin = pml_fields_xmin->By_; + Bz_pml_xmin = pml_fields_xmin->Bz_; + } + + if(ncells_pml_xmax != 0){ + Ex_pml_xmax = pml_fields_xmax->Ex_; + Ey_pml_xmax = pml_fields_xmax->Ey_; + Ez_pml_xmax = pml_fields_xmax->Ez_; + Hx_pml_xmax = pml_fields_xmax->Hx_; + Hy_pml_xmax = pml_fields_xmax->Hy_; + Hz_pml_xmax = pml_fields_xmax->Hz_; + Dx_pml_xmax = pml_fields_xmax->Dx_; + Dy_pml_xmax = pml_fields_xmax->Dy_; + Dz_pml_xmax = pml_fields_xmax->Dz_; + Bx_pml_xmax = pml_fields_xmax->Bx_; + By_pml_xmax = pml_fields_xmax->By_; + Bz_pml_xmax = pml_fields_xmax->Bz_; + } + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + // 2. Exchange field PML <- Domain + for ( int j=min2exchange ; jisXmin()) { + if(ncells_pml_xmin != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + //Injecting a laser + vector bx( n_p[0]*n_d[2], 0. ); + for( unsigned int i=patch->isXmin() ; iisXmax() ; i++ ) { + pos[0] = patch->getDomainLocalMin( 0 ) + ( ( int )i - ( int )EMfields->oversize[0] )*d[0]; + for( unsigned int k=patch->isZmin() ; kisZmax() ; k++ ) { + pos[1] = patch->getDomainLocalMin( 2 ) + ( ( int )k -0.5 - ( int )EMfields->oversize[2] )*d[2]; + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bx[ i*n_d[2]+k ] += vecLaser[ilaser]->getAmplitude0( pos, time_dual, i, k ); + } + (*Hx_)( ncells_pml_xmin+i, ncells_pml_domain-domain_oversize_y-nsolver/2, k ) += factor_laser_angle_S*bx[ i*n_d[2]+k ] ; + (*Bx_)( ncells_pml_xmin+i, ncells_pml_domain-domain_oversize_y-nsolver/2, k ) += factor_laser_angle_S*bx[ i*n_d[2]+k ] ; + } + } + + vector bz( n_d[0]*n_p[2], 0. ); + for( unsigned int i=patch->isXmin() ; iisXmax() ; i++ ) { + pos[0] = patch->getDomainLocalMin( 0 ) + ( ( int )i -0.5 - ( int )EMfields->oversize[0] )*d[0]; + for( unsigned int k=patch->isZmin() ; kisZmax() ; k++ ) { + pos[1] = patch->getDomainLocalMin( 2 ) + ( ( int )k - ( int )EMfields->oversize[2] )*d[2]; + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bz[ i*n_p[2]+k ] += vecLaser[ilaser]->getAmplitude1( pos, time_dual, i, k ); + } + (*Hz_)( ncells_pml_xmin+i, ncells_pml_domain-domain_oversize_y-nsolver/2, k ) += factor_laser_angle_S*bz[ i*n_p[2]+k ] ; + (*Bz_)( ncells_pml_xmin+i, ncells_pml_domain-domain_oversize_y-nsolver/2, k ) += factor_laser_angle_S*bz[ i*n_p[2]+k ] ; + } + } + + // 4. Exchange PML -> Domain + // Primals in y-direction + for (int j=0 ; j < nsolver/2 ; j++){ + for ( int i=0 ; i PML x MIN + // Primal in y-direction + for (int j=0 ; j < nsolver/2 ; j++){ + if (patch->isXmin()) { + if(ncells_pml_xmin != 0){ + for ( int i=0 ; iisXmin()) { + if(ncells_pml_xmin != 0){ + for ( int i=0 ; i PML x MAX + // Primal in y-direction + for (int j=0 ; j < nsolver/2 ; j++){ + if (patch->isXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; iisYmax() ) { + + ElectroMagnBC3D_PML* pml_fields_xmin = NULL ; + ElectroMagnBC3D_PML* pml_fields_xmax = NULL ; + + if(ncells_pml_xmin != 0){ + pml_fields_xmin = static_cast( EMfields->emBoundCond[0] ); + } + if(ncells_pml_xmax != 0){ + pml_fields_xmax = static_cast( EMfields->emBoundCond[1] ); + } + + Field3D* Ex_pml_xmin = NULL; + Field3D* Ey_pml_xmin = NULL; + Field3D* Ez_pml_xmin = NULL; + Field3D* Hx_pml_xmin = NULL; + Field3D* Hy_pml_xmin = NULL; + Field3D* Hz_pml_xmin = NULL; + Field3D* Dx_pml_xmin = NULL; + Field3D* Dy_pml_xmin = NULL; + Field3D* Dz_pml_xmin = NULL; + Field3D* Bx_pml_xmin = NULL; + Field3D* By_pml_xmin = NULL; + Field3D* Bz_pml_xmin = NULL; + + Field3D* Ex_pml_xmax = NULL; + Field3D* Ey_pml_xmax = NULL; + Field3D* Ez_pml_xmax = NULL; + Field3D* Hx_pml_xmax = NULL; + Field3D* Hy_pml_xmax = NULL; + Field3D* Hz_pml_xmax = NULL; + Field3D* Dx_pml_xmax = NULL; + Field3D* Dy_pml_xmax = NULL; + Field3D* Dz_pml_xmax = NULL; + Field3D* Bx_pml_xmax = NULL; + Field3D* By_pml_xmax = NULL; + Field3D* Bz_pml_xmax = NULL; + + if(ncells_pml_xmin != 0){ + Ex_pml_xmin = pml_fields_xmin->Ex_; + Ey_pml_xmin = pml_fields_xmin->Ey_; + Ez_pml_xmin = pml_fields_xmin->Ez_; + Hx_pml_xmin = pml_fields_xmin->Hx_; + Hy_pml_xmin = pml_fields_xmin->Hy_; + Hz_pml_xmin = pml_fields_xmin->Hz_; + Dx_pml_xmin = pml_fields_xmin->Dx_; + Dy_pml_xmin = pml_fields_xmin->Dy_; + Dz_pml_xmin = pml_fields_xmin->Dz_; + Bx_pml_xmin = pml_fields_xmin->Bx_; + By_pml_xmin = pml_fields_xmin->By_; + Bz_pml_xmin = pml_fields_xmin->Bz_; + } + + if(ncells_pml_xmax != 0){ + Ex_pml_xmax = pml_fields_xmax->Ex_; + Ey_pml_xmax = pml_fields_xmax->Ey_; + Ez_pml_xmax = pml_fields_xmax->Ez_; + Hx_pml_xmax = pml_fields_xmax->Hx_; + Hy_pml_xmax = pml_fields_xmax->Hy_; + Hz_pml_xmax = pml_fields_xmax->Hz_; + Dx_pml_xmax = pml_fields_xmax->Dx_; + Dy_pml_xmax = pml_fields_xmax->Dy_; + Dz_pml_xmax = pml_fields_xmax->Dz_; + Bx_pml_xmax = pml_fields_xmax->Bx_; + By_pml_xmax = pml_fields_xmax->By_; + Bz_pml_xmax = pml_fields_xmax->Bz_; + } + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + // 2. Exchange field PML <- Domain + for ( int j=min2exchange ; jisXmin()) { + if(ncells_pml_xmin != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + //Injecting a laser + vector bx( n_p[0]*n_d[2], 0. ); + for( unsigned int i=patch->isXmin() ; iisXmax() ; i++ ) { + pos[0] = patch->getDomainLocalMin( 0 ) + ( ( int )i - ( int )EMfields->oversize[0] )*d[0]; + for( unsigned int k=patch->isZmin() ; kisZmax() ; k++ ) { + pos[1] = patch->getDomainLocalMin( 2 ) + ( ( int )k -0.5 - ( int )EMfields->oversize[2] )*d[2]; + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bx[ i*n_d[2]+k ] += vecLaser[ilaser]->getAmplitude0( pos, time_dual, i, k ); + } + (*Hx_)( ncells_pml_xmin+i , domain_oversize_y+nsolver/2, k ) += factor_laser_angle_N*bx[ i*n_d[2]+k ] ; + (*Bx_)( ncells_pml_xmin+i , domain_oversize_y+nsolver/2, k ) += factor_laser_angle_N*bx[ i*n_d[2]+k ] ; + } + } + + vector bz( n_d[0]*n_p[2], 0. ); + for( unsigned int i=patch->isXmin() ; iisXmax() ; i++ ) { + pos[0] = patch->getDomainLocalMin( 0 ) + ( ( int )i -0.5 - ( int )EMfields->oversize[0] )*d[0]; + for( unsigned int k=patch->isZmin() ; kisZmax() ; k++ ) { + pos[1] = patch->getDomainLocalMin( 2 ) + ( ( int )k - ( int )EMfields->oversize[2] )*d[2]; + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bz[ i*n_p[2]+k ] += vecLaser[ilaser]->getAmplitude1( pos, time_dual, i, k ); + } + (*Hz_)( ncells_pml_xmin+i, domain_oversize_y+nsolver/2, k ) += factor_laser_angle_N*bz[ i*n_p[2]+k ] ; + (*Bz_)( ncells_pml_xmin+i, domain_oversize_y+nsolver/2, k ) += factor_laser_angle_N*bz[ i*n_p[2]+k ] ; + } + } + + // 4. Exchange PML -> Domain + // Primals in y-direction + for (int j=0 ; j < nsolver/2-1 ; j++){ + for ( int i=0 ; i PML x MIN + // Primals in y-direction + for (int j=0 ; j < nsolver/2-1 ; j++){ + if (patch->isXmin()) { + if(ncells_pml_xmin != 0){ + for ( int i=0 ; iisXmin()) { + if(ncells_pml_xmin != 0){ + for ( int i=0 ; i PML x MAX + // Primals in y-direction + for (int j=0 ; j < nsolver/2-1 ; j++){ + if (patch->isXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; iisZmin() ) { + + ElectroMagnBC3D_PML* pml_fields_xmin = NULL ; + ElectroMagnBC3D_PML* pml_fields_xmax = NULL ; + ElectroMagnBC3D_PML* pml_fields_ymin = NULL ; + ElectroMagnBC3D_PML* pml_fields_ymax = NULL ; + + if(ncells_pml_xmin != 0){ + pml_fields_xmin = static_cast( EMfields->emBoundCond[0] ); + } + if(ncells_pml_xmax != 0){ + pml_fields_xmax = static_cast( EMfields->emBoundCond[1] ); + } + if(ncells_pml_ymin != 0){ + pml_fields_ymin = static_cast( EMfields->emBoundCond[2] ); + } + if(ncells_pml_ymax != 0){ + pml_fields_ymax = static_cast( EMfields->emBoundCond[3] ); + } + + Field3D* Ex_pml_xmin = NULL; + Field3D* Ey_pml_xmin = NULL; + Field3D* Ez_pml_xmin = NULL; + Field3D* Hx_pml_xmin = NULL; + Field3D* Hy_pml_xmin = NULL; + Field3D* Hz_pml_xmin = NULL; + Field3D* Dx_pml_xmin = NULL; + Field3D* Dy_pml_xmin = NULL; + Field3D* Dz_pml_xmin = NULL; + Field3D* Bx_pml_xmin = NULL; + Field3D* By_pml_xmin = NULL; + Field3D* Bz_pml_xmin = NULL; + Field3D* Ex_pml_xmax = NULL; + Field3D* Ey_pml_xmax = NULL; + Field3D* Ez_pml_xmax = NULL; + Field3D* Hx_pml_xmax = NULL; + Field3D* Hy_pml_xmax = NULL; + Field3D* Hz_pml_xmax = NULL; + Field3D* Dx_pml_xmax = NULL; + Field3D* Dy_pml_xmax = NULL; + Field3D* Dz_pml_xmax = NULL; + Field3D* Bx_pml_xmax = NULL; + Field3D* By_pml_xmax = NULL; + Field3D* Bz_pml_xmax = NULL; + + Field3D* Ex_pml_ymin = NULL; + Field3D* Ey_pml_ymin = NULL; + Field3D* Ez_pml_ymin = NULL; + Field3D* Hx_pml_ymin = NULL; + Field3D* Hy_pml_ymin = NULL; + Field3D* Hz_pml_ymin = NULL; + Field3D* Dx_pml_ymin = NULL; + Field3D* Dy_pml_ymin = NULL; + Field3D* Dz_pml_ymin = NULL; + Field3D* Bx_pml_ymin = NULL; + Field3D* By_pml_ymin = NULL; + Field3D* Bz_pml_ymin = NULL; + Field3D* Ex_pml_ymax = NULL; + Field3D* Ey_pml_ymax = NULL; + Field3D* Ez_pml_ymax = NULL; + Field3D* Hx_pml_ymax = NULL; + Field3D* Hy_pml_ymax = NULL; + Field3D* Hz_pml_ymax = NULL; + Field3D* Dx_pml_ymax = NULL; + Field3D* Dy_pml_ymax = NULL; + Field3D* Dz_pml_ymax = NULL; + Field3D* Bx_pml_ymax = NULL; + Field3D* By_pml_ymax = NULL; + Field3D* Bz_pml_ymax = NULL; + + if(ncells_pml_xmin != 0){ + Ex_pml_xmin = pml_fields_xmin->Ex_; + Ey_pml_xmin = pml_fields_xmin->Ey_; + Ez_pml_xmin = pml_fields_xmin->Ez_; + Hx_pml_xmin = pml_fields_xmin->Hx_; + Hy_pml_xmin = pml_fields_xmin->Hy_; + Hz_pml_xmin = pml_fields_xmin->Hz_; + Dx_pml_xmin = pml_fields_xmin->Dx_; + Dy_pml_xmin = pml_fields_xmin->Dy_; + Dz_pml_xmin = pml_fields_xmin->Dz_; + Bx_pml_xmin = pml_fields_xmin->Bx_; + By_pml_xmin = pml_fields_xmin->By_; + Bz_pml_xmin = pml_fields_xmin->Bz_; + } + + if(ncells_pml_xmax != 0){ + Ex_pml_xmax = pml_fields_xmax->Ex_; + Ey_pml_xmax = pml_fields_xmax->Ey_; + Ez_pml_xmax = pml_fields_xmax->Ez_; + Hx_pml_xmax = pml_fields_xmax->Hx_; + Hy_pml_xmax = pml_fields_xmax->Hy_; + Hz_pml_xmax = pml_fields_xmax->Hz_; + Dx_pml_xmax = pml_fields_xmax->Dx_; + Dy_pml_xmax = pml_fields_xmax->Dy_; + Dz_pml_xmax = pml_fields_xmax->Dz_; + Bx_pml_xmax = pml_fields_xmax->Bx_; + By_pml_xmax = pml_fields_xmax->By_; + Bz_pml_xmax = pml_fields_xmax->Bz_; + } + + if(ncells_pml_ymin != 0){ + Ex_pml_ymin = pml_fields_ymin->Ex_; + Ey_pml_ymin = pml_fields_ymin->Ey_; + Ez_pml_ymin = pml_fields_ymin->Ez_; + Hx_pml_ymin = pml_fields_ymin->Hx_; + Hy_pml_ymin = pml_fields_ymin->Hy_; + Hz_pml_ymin = pml_fields_ymin->Hz_; + Dx_pml_ymin = pml_fields_ymin->Dx_; + Dy_pml_ymin = pml_fields_ymin->Dy_; + Dz_pml_ymin = pml_fields_ymin->Dz_; + Bx_pml_ymin = pml_fields_ymin->Bx_; + By_pml_ymin = pml_fields_ymin->By_; + Bz_pml_ymin = pml_fields_ymin->Bz_; + } + + if(ncells_pml_ymax != 0){ + Ex_pml_ymax = pml_fields_ymax->Ex_; + Ey_pml_ymax = pml_fields_ymax->Ey_; + Ez_pml_ymax = pml_fields_ymax->Ez_; + Hx_pml_ymax = pml_fields_ymax->Hx_; + Hy_pml_ymax = pml_fields_ymax->Hy_; + Hz_pml_ymax = pml_fields_ymax->Hz_; + Dx_pml_ymax = pml_fields_ymax->Dx_; + Dy_pml_ymax = pml_fields_ymax->Dy_; + Dz_pml_ymax = pml_fields_ymax->Dz_; + Bx_pml_ymax = pml_fields_ymax->Bx_; + By_pml_ymax = pml_fields_ymax->By_; + Bz_pml_ymax = pml_fields_ymax->Bz_; + } + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + /* + Dans le bloc du dessous (2.) + Peut être faire passer les echanges XMAX et YMAX apres l'echange du domaine + */ + + // 2. Exchange field PML <- Domain + for ( int k=min2exchange ; kisXmin()) { + if(ncells_pml_xmin != 0){ + int idx_start = 0; + int jdx_start = ncells_pml_ymin; + // Les qtes i-Primals commencent a 0 + // Toutes les qtes i-Duals 0 + for ( int j=0 ; jisYmin()) { + if(ncells_pml_ymin != 0){ + int jdx_start = 0; + // Les qtes j-Primals commencent a 0 + // Toutes les qtes j-Duals 0 + for ( int i=0 ; iisXmax()) { + if(ncells_pml_xmax != 0){ + int idx_start = (zpml_size_in_x-1)-(ncells_pml_xmax-1) ; + int jdx_start = ncells_pml_ymin; + // Les qtes i-Primals commencent a (zpml_size_in_x+1)-ncells_pml_xmax + // Toutes les qtes i-Duals commence a idx_start+1 + for ( int j=0 ; jisYmax()) { + if(ncells_pml_ymax != 0){ + int jdx_start = (zpml_size_in_y-1)-(ncells_pml_ymax-1) ; + // Les qtes j-Primals commencent a (zpml_size_in_y+1)-ncells_pml_ymax + // Toutes les qtes j-Duals commence a jdx_start+1 + for ( int i=0 ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + //Injecting a laser + vector bx( n_p[0]*n_d[1], 0. ); // Bx(p,d,d) + for( unsigned int i=patch->isXmin() ; iisXmax() ; i++ ) { + pos[0] = patch->getDomainLocalMin( 0 ) + ( ( int )i - ( int )EMfields->oversize[0] )*d[0]; + for( unsigned int j=patch->isYmin() ; jisYmax() ; j++ ) { + pos[1] = patch->getDomainLocalMin( 1 ) + ( ( int )j - 0.5 - ( int )EMfields->oversize[1] )*d[1]; + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bx[ i*n_d[1]+j ] += vecLaser[ilaser]->getAmplitude0( pos, time_dual, i, j ); + } + (*Hx_)( ncells_pml_xmin+i, ncells_pml_ymin+j, ncells_pml_domain-domain_oversize_z-nsolver/2 ) += factor_laser_angle_B*bx[ i*n_d[1]+j ] ; + (*Bx_)( ncells_pml_xmin+i, ncells_pml_ymin+j, ncells_pml_domain-domain_oversize_z-nsolver/2 ) += factor_laser_angle_B*bx[ i*n_d[1]+j ] ; + } + } + + vector by( n_d[0]*n_p[1], 0. ); // By(d,p,d) + for( unsigned int i=patch->isXmin() ; iisXmax() ; i++ ) { + pos[0] = patch->getDomainLocalMin( 0 ) + ( ( int )i -0.5 - ( int )EMfields->oversize[0] )*d[0]; + for( unsigned int j=patch->isYmin() ; jisYmax() ; j++ ) { + pos[1] = patch->getDomainLocalMin( 1 ) + ( ( int )j - ( int )EMfields->oversize[1] )*d[1]; + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + by[ i*n_p[1]+j ] += vecLaser[ilaser]->getAmplitude1( pos, time_dual, i, j ); + } + (*Hy_)( ncells_pml_xmin+i, ncells_pml_ymin+j, ncells_pml_domain-domain_oversize_z-nsolver/2 ) += factor_laser_angle_B*by[ i*n_p[1]+j ] ; + (*By_)( ncells_pml_xmin+i, ncells_pml_ymin+j, ncells_pml_domain-domain_oversize_z-nsolver/2 ) += factor_laser_angle_B*by[ i*n_p[1]+j ] ; + } + } + + // 4. Exchange PML -> Domain + // Ici il faut a priori remplacer y<->z et j<->k + // Primals in z-direction + for (int k=0 ; k < nsolver/2 ; k++){ + int idx_start = ncells_pml_xmin; + int jdx_start = ncells_pml_ymin; + for ( int i=0 ; i PML xMIN + // Primal in z-direction + for (int k=0 ; k < nsolver/2 ; k++){ + if (patch->isXmin()) { + if(ncells_pml_xmin != 0){ + int idx_start = 0; + int jdx_start = ncells_pml_ymin; + for ( int i=0 ; iisXmin()) { + if(ncells_pml_xmin != 0){ + int idx_start = 0; + int jdx_start = ncells_pml_ymin; + for ( int i=0 ; i PML xMAX + // Primal in z-direction + for (int k=0 ; k < nsolver/2 ; k++){ + if (patch->isXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; i PML yMin + // Primal in z-direction + for (int k=0 ; k < nsolver/2 ; k++){ + if (patch->isYmin()) { + if(ncells_pml_ymin != 0){ + int jdx_start = 0; + for ( int j=0 ; jisYmin()) { + if(ncells_pml_ymin != 0){ + int jdx_start = 0; + for ( int j=0 ; j PML yMax + // Primal in z-direction + for (int k=0 ; k < nsolver/2 ; k++){ + if (patch->isYmax()) { + if(ncells_pml_ymax != 0){ + for ( int j=0 ; jisYmax()) { + if(ncells_pml_ymax != 0){ + for ( int j=0 ; jisZmax() ) { + + ElectroMagnBC3D_PML* pml_fields_xmin = NULL ; + ElectroMagnBC3D_PML* pml_fields_xmax = NULL ; + ElectroMagnBC3D_PML* pml_fields_ymin = NULL ; + ElectroMagnBC3D_PML* pml_fields_ymax = NULL ; + + if(ncells_pml_xmin != 0){ + pml_fields_xmin = static_cast( EMfields->emBoundCond[0] ); + } + if(ncells_pml_xmax != 0){ + pml_fields_xmax = static_cast( EMfields->emBoundCond[1] ); + } + if(ncells_pml_ymin != 0){ + pml_fields_ymin = static_cast( EMfields->emBoundCond[2] ); + } + if(ncells_pml_ymax != 0){ + pml_fields_ymax = static_cast( EMfields->emBoundCond[3] ); + } + + Field3D* Ex_pml_xmin = NULL; + Field3D* Ey_pml_xmin = NULL; + Field3D* Ez_pml_xmin = NULL; + Field3D* Hx_pml_xmin = NULL; + Field3D* Hy_pml_xmin = NULL; + Field3D* Hz_pml_xmin = NULL; + Field3D* Dx_pml_xmin = NULL; + Field3D* Dy_pml_xmin = NULL; + Field3D* Dz_pml_xmin = NULL; + Field3D* Bx_pml_xmin = NULL; + Field3D* By_pml_xmin = NULL; + Field3D* Bz_pml_xmin = NULL; + Field3D* Ex_pml_xmax = NULL; + Field3D* Ey_pml_xmax = NULL; + Field3D* Ez_pml_xmax = NULL; + Field3D* Hx_pml_xmax = NULL; + Field3D* Hy_pml_xmax = NULL; + Field3D* Hz_pml_xmax = NULL; + Field3D* Dx_pml_xmax = NULL; + Field3D* Dy_pml_xmax = NULL; + Field3D* Dz_pml_xmax = NULL; + Field3D* Bx_pml_xmax = NULL; + Field3D* By_pml_xmax = NULL; + Field3D* Bz_pml_xmax = NULL; + + Field3D* Ex_pml_ymin = NULL; + Field3D* Ey_pml_ymin = NULL; + Field3D* Ez_pml_ymin = NULL; + Field3D* Hx_pml_ymin = NULL; + Field3D* Hy_pml_ymin = NULL; + Field3D* Hz_pml_ymin = NULL; + Field3D* Dx_pml_ymin = NULL; + Field3D* Dy_pml_ymin = NULL; + Field3D* Dz_pml_ymin = NULL; + Field3D* Bx_pml_ymin = NULL; + Field3D* By_pml_ymin = NULL; + Field3D* Bz_pml_ymin = NULL; + Field3D* Ex_pml_ymax = NULL; + Field3D* Ey_pml_ymax = NULL; + Field3D* Ez_pml_ymax = NULL; + Field3D* Hx_pml_ymax = NULL; + Field3D* Hy_pml_ymax = NULL; + Field3D* Hz_pml_ymax = NULL; + Field3D* Dx_pml_ymax = NULL; + Field3D* Dy_pml_ymax = NULL; + Field3D* Dz_pml_ymax = NULL; + Field3D* Bx_pml_ymax = NULL; + Field3D* By_pml_ymax = NULL; + Field3D* Bz_pml_ymax = NULL; + + if(ncells_pml_xmin != 0){ + Ex_pml_xmin = pml_fields_xmin->Ex_; + Ey_pml_xmin = pml_fields_xmin->Ey_; + Ez_pml_xmin = pml_fields_xmin->Ez_; + Hx_pml_xmin = pml_fields_xmin->Hx_; + Hy_pml_xmin = pml_fields_xmin->Hy_; + Hz_pml_xmin = pml_fields_xmin->Hz_; + Dx_pml_xmin = pml_fields_xmin->Dx_; + Dy_pml_xmin = pml_fields_xmin->Dy_; + Dz_pml_xmin = pml_fields_xmin->Dz_; + Bx_pml_xmin = pml_fields_xmin->Bx_; + By_pml_xmin = pml_fields_xmin->By_; + Bz_pml_xmin = pml_fields_xmin->Bz_; + } + + if(ncells_pml_xmax != 0){ + Ex_pml_xmax = pml_fields_xmax->Ex_; + Ey_pml_xmax = pml_fields_xmax->Ey_; + Ez_pml_xmax = pml_fields_xmax->Ez_; + Hx_pml_xmax = pml_fields_xmax->Hx_; + Hy_pml_xmax = pml_fields_xmax->Hy_; + Hz_pml_xmax = pml_fields_xmax->Hz_; + Dx_pml_xmax = pml_fields_xmax->Dx_; + Dy_pml_xmax = pml_fields_xmax->Dy_; + Dz_pml_xmax = pml_fields_xmax->Dz_; + Bx_pml_xmax = pml_fields_xmax->Bx_; + By_pml_xmax = pml_fields_xmax->By_; + Bz_pml_xmax = pml_fields_xmax->Bz_; + } + + if(ncells_pml_ymin != 0){ + Ex_pml_ymin = pml_fields_ymin->Ex_; + Ey_pml_ymin = pml_fields_ymin->Ey_; + Ez_pml_ymin = pml_fields_ymin->Ez_; + Hx_pml_ymin = pml_fields_ymin->Hx_; + Hy_pml_ymin = pml_fields_ymin->Hy_; + Hz_pml_ymin = pml_fields_ymin->Hz_; + Dx_pml_ymin = pml_fields_ymin->Dx_; + Dy_pml_ymin = pml_fields_ymin->Dy_; + Dz_pml_ymin = pml_fields_ymin->Dz_; + Bx_pml_ymin = pml_fields_ymin->Bx_; + By_pml_ymin = pml_fields_ymin->By_; + Bz_pml_ymin = pml_fields_ymin->Bz_; + } + + if(ncells_pml_ymax != 0){ + Ex_pml_ymax = pml_fields_ymax->Ex_; + Ey_pml_ymax = pml_fields_ymax->Ey_; + Ez_pml_ymax = pml_fields_ymax->Ez_; + Hx_pml_ymax = pml_fields_ymax->Hx_; + Hy_pml_ymax = pml_fields_ymax->Hy_; + Hz_pml_ymax = pml_fields_ymax->Hz_; + Dx_pml_ymax = pml_fields_ymax->Dx_; + Dy_pml_ymax = pml_fields_ymax->Dy_; + Dz_pml_ymax = pml_fields_ymax->Dz_; + Bx_pml_ymax = pml_fields_ymax->Bx_; + By_pml_ymax = pml_fields_ymax->By_; + Bz_pml_ymax = pml_fields_ymax->Bz_; + } + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + /* + Dans le bloc du dessous (2.) + Peut être faire passer les echanges XMAX et YMAX apres l'echange du domaine + */ + + // 2. Exchange field PML <- Domain + for ( int k=min2exchange ; kisXmin()) { + if(ncells_pml_xmin != 0){ + int idx_start = 0; + int jdx_start = ncells_pml_ymin; + for ( int i=0 ; iisYmin()) { + if(ncells_pml_ymin != 0){ + int jdx_start = 0; + for ( int i=0 ; iisXmax()) { + if(ncells_pml_xmax != 0){ + int idx_start = (zpml_size_in_x-1)-(ncells_pml_xmax-1) ; + int jdx_start = ncells_pml_ymin; + for ( int i=0 ; iisYmax()) { + if(ncells_pml_ymax != 0){ + int jdx_start = (zpml_size_in_y-1)-(ncells_pml_ymax-1) ; + // Les qtes j-Primals commencent a (zpml_size_in_y+1)-ncells_pml_ymax + // Toutes les qtes j-Duals commence a jdx_start+1 + for ( int i=0 ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + //Injecting a laser + vector bx( n_p[0]*n_d[1], 0. ); // Bx(p,d,d) + for( unsigned int i=patch->isXmin() ; iisXmax() ; i++ ) { + pos[0] = patch->getDomainLocalMin( 0 ) + ( ( int )i - ( int )EMfields->oversize[0] )*d[0]; + for( unsigned int j=patch->isYmin() ; jisYmax() ; j++ ) { + pos[1] = patch->getDomainLocalMin( 1 ) + ( ( int )j - 0.5 - ( int )EMfields->oversize[1] )*d[1]; + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + bx[ i*n_d[1]+j ] += vecLaser[ilaser]->getAmplitude0( pos, time_dual, i, j ); + } + (*Hx_)( ncells_pml_xmin+i, ncells_pml_ymin+j, domain_oversize_z+nsolver/2 ) += factor_laser_angle_T*bx[ i*n_d[1]+j ] ; + (*Bx_)( ncells_pml_xmin+i, ncells_pml_ymin+j, domain_oversize_z+nsolver/2 ) += factor_laser_angle_T*bx[ i*n_d[1]+j ] ; + } + } + + vector by( n_d[0]*n_p[1], 0. ); // By(d,p,d) + for( unsigned int i=patch->isXmin() ; iisXmax() ; i++ ) { + pos[0] = patch->getDomainLocalMin( 0 ) + ( ( int )i -0.5 - ( int )EMfields->oversize[0] )*d[0]; + for( unsigned int j=patch->isYmin() ; jisYmax() ; j++ ) { + pos[1] = patch->getDomainLocalMin( 1 ) + ( ( int )j - ( int )EMfields->oversize[1] )*d[1]; + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + by[ i*n_p[1]+j ] += vecLaser[ilaser]->getAmplitude1( pos, time_dual, i, j ); + } + (*Hy_)( ncells_pml_xmin+i, ncells_pml_ymin+j, domain_oversize_z+nsolver/2 ) += factor_laser_angle_T*by[ i*n_p[1]+j ] ; + (*By_)( ncells_pml_xmin+i, ncells_pml_ymin+j, domain_oversize_z+nsolver/2 ) += factor_laser_angle_T*by[ i*n_p[1]+j ] ; + } + } + + // 4. Exchange PML -> Domain + // Ici il faut a priori remplacer y<->z et j<->k + // Primals in z-direction + for (int k=0 ; k < nsolver/2-1 ; k++){ + int idx_start = ncells_pml_xmin; + int jdx_start = ncells_pml_ymin; + for ( int i=0 ; iz for PML xMIN + // 5. Exchange PML z -> PML xMIN + // Primal in z-direction + for (int k=0 ; k < nsolver/2-1 ; k++){ + if (patch->isXmin()) { + if(ncells_pml_xmin != 0){ + int idx_start = 0; + int jdx_start = ncells_pml_ymin; + for ( int i=0 ; iisXmin()) { + if(ncells_pml_xmin != 0){ + int idx_start = 0; + int jdx_start = ncells_pml_ymin; + for ( int i=0 ; iz for PML xMax + // 6. Exchange PML z -> PML xMAX + // Primal in z-direction + for (int k=0 ; k < nsolver/2-1 ; k++){ + if (patch->isXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_xmax != 0){ + for ( int i=0 ; i PML yMin + // Here Change x<->y with respect to previous block 5.for PML xMin -> PML yMin + // Primal in z-direction + for (int k=0 ; k < nsolver/2-1 ; k++){ + if (patch->isYmin()) { + if(ncells_pml_ymin != 0){ + int jdx_start = 0; + for ( int j=0 ; jisYmin()) { + if(ncells_pml_ymin != 0){ + int jdx_start = 0; + for ( int j=0 ; j PML yMax + // Here Change x<->y with respect to previous block 6.for PML xMax -> PML yMax + // Primal in z-direction + for (int k=0 ; k < nsolver/2-1 ; k++){ + if (patch->isYmax()) { + if(ncells_pml_ymax != 0){ + for ( int j=0 ; jisYmax()) { + if(ncells_pml_ymax != 0){ + for ( int j=0 ; j +#include "Tools.h" +#include "ElectroMagnBC3D.h" +#include "ElectroMagn3D.h" +#include "Field3D.h" + +class Params; +class ElectroMagn; +class Field; +class PML_Solver3D_Yee; +class PML_Solver3D_Bouchard; + +class ElectroMagnBC3D_PML : public ElectroMagnBC3D +{ +public: + + ElectroMagnBC3D_PML( Params ¶ms, Patch *patch, unsigned int _min_max ); + ~ElectroMagnBC3D_PML(); + + virtual void apply( ElectroMagn *EMfields, double time_dual, Patch *patch ) override; + + void save_fields( Field *, Patch *patch ) override; + void disableExternalFields() override; + + Field3D* Ex_ = NULL; + Field3D* Ey_ = NULL; + Field3D* Ez_ = NULL; + Field3D* Bx_ = NULL; + Field3D* By_ = NULL; + Field3D* Bz_ = NULL; + Field3D* Dx_ = NULL; + Field3D* Dy_ = NULL; + Field3D* Dz_ = NULL; + Field3D* Hx_ = NULL; + Field3D* Hy_ = NULL; + Field3D* Hz_ = NULL; + + Field* getExPML() override { return Ex_; }; + Field* getEyPML() override { return Ey_; }; + Field* getEzPML() override { return Ez_; }; + Field* getBxPML() override { return Bx_; }; + Field* getByPML() override { return By_; }; + Field* getBzPML() override { return Bz_; }; + Field* getDxPML() override { return Dx_; }; + Field* getDyPML() override { return Dy_; }; + Field* getDzPML() override { return Dz_; }; + Field* getHxPML() override { return Hx_; }; + Field* getHyPML() override { return Hy_; }; + Field* getHzPML() override { return Hz_; }; + + int domain_oversize_x; + int domain_oversize_y; + int domain_oversize_z; + + int ncells_pml_xmin ; + int ncells_pml_xmax ; + int ncells_pml_ymin ; + int ncells_pml_ymax ; + int ncells_pml_zmin ; + int ncells_pml_zmax ; + + int ncells_pml_domain_xmin; + int ncells_pml_domain_xmax; + int ncells_pml_domain_ymin; + int ncells_pml_domain_ymax; + int ncells_pml_domain_zmin; + int ncells_pml_domain_zmax; + + int ncells_pml; + int ncells_pml_domain; + + int nsolver; + + int min2exchange; + int max2exchange; + int solvermin; + int solvermax; + + int ypml_size_in_x; + int zpml_size_in_x; + int zpml_size_in_y; + int startpml; + + double factor_laser_space_time ; + double factor_laser_angle_W ; + double factor_laser_angle_E ; + double factor_laser_angle_S ; + double factor_laser_angle_N ; + double factor_laser_angle_B ; + double factor_laser_angle_T ; +}; + +#endif diff --git a/src/ElectroMagnBC/ElectroMagnBCAM_BM.cpp b/src/ElectroMagnBC/ElectroMagnBCAM_BM.cpp index 31d02a36d..a07635ff1 100755 --- a/src/ElectroMagnBC/ElectroMagnBCAM_BM.cpp +++ b/src/ElectroMagnBC/ElectroMagnBCAM_BM.cpp @@ -169,7 +169,6 @@ void ElectroMagnBCAM_BM::disableExternalFields() // --------------------------------------------------------------------------------------------------------------------- void ElectroMagnBCAM_BM::apply( ElectroMagn *EMfields, double time_dual, Patch *patch ) { - //This condition can only be applied to Rmax // Loop on imode diff --git a/src/ElectroMagnBC/ElectroMagnBCAM_PML.cpp b/src/ElectroMagnBC/ElectroMagnBCAM_PML.cpp new file mode 100755 index 000000000..e37f8d3aa --- /dev/null +++ b/src/ElectroMagnBC/ElectroMagnBCAM_PML.cpp @@ -0,0 +1,1105 @@ +#include "ElectroMagnBCAM_PML.h" + +#include + +#include +#include + +#include "Params.h" +#include "Patch.h" +#include "ElectroMagn.h" +#include "FieldFactory.h" +#include "cField2D.h" +#include "Tools.h" +#include "Laser.h" +#include +#include "dcomplex.h" + +#include "SolverFactory.h" + +using namespace std; + +ElectroMagnBCAM_PML::ElectroMagnBCAM_PML( Params ¶ms, Patch *patch, unsigned int i_boundary ) + : ElectroMagnBCAM( params, patch, i_boundary ) +{ + std::vector n_space(params.n_space); + std::vector oversize(params.oversize); + if( params.multiple_decomposition ) { + n_space = params.n_space_region; + oversize = params.region_oversize; + } + + //Number of modes + Nmode = params.nmodes; + + pml_solver_ = SolverFactory::createPML( params ); + if (params.maxwell_sol=="Yee"){ + nsolver=2; + //MESSAGE("FDTD scheme in PML region : Yee."); + } + else { + //WARNING("The solver you use in the main domain is not the same as in the PML region. FDTD scheme in PML region : Yee."); + nsolver=2; + } + + El_.resize( Nmode ); + Dl_.resize( Nmode ); + Hl_.resize( Nmode ); + Bl_.resize( Nmode ); + Er_.resize( Nmode ); + Dr_.resize( Nmode ); + Hr_.resize( Nmode ); + Br_.resize( Nmode ); + Et_.resize( Nmode ); + Dt_.resize( Nmode ); + Ht_.resize( Nmode ); + Bt_.resize( Nmode ); + + if ( ( i_boundary_ == 0 && patch->isXmin() ) + || ( i_boundary_ == 1 && patch->isXmax() ) + || ( i_boundary_ == 2 && patch->isYmin() ) + || ( i_boundary_ == 3 && patch->isYmax() ) ) { + + int iDim = 0*((i_boundary_==0)||(i_boundary_==1))+1*((i_boundary_==2)||(i_boundary_==3)); + int min_or_max = (i_boundary_)%2; + + domain_oversize_l = oversize[0] ; + domain_oversize_r = oversize[1] ; + + if (patch->isXmin() ) {//&& i_boundary_ == 0 ) { + ncells_pml_lmin = params.number_of_pml_cells[0][0]; + ncells_pml_domain_lmin = ncells_pml_lmin + 1*oversize[0] + nsolver/2; + domain_oversize_l = oversize[0] ; + } + else { + ncells_pml_lmin = 0; + ncells_pml_domain_lmin = 0; + } + if (patch->isXmax() ) {//&& i_boundary_ == 1 ) { + ncells_pml_lmax = params.number_of_pml_cells[0][1]; + ncells_pml_domain_lmax = ncells_pml_lmax + 1*oversize[0] + nsolver/2; + domain_oversize_l = oversize[0] ; + } + else { + ncells_pml_lmax = 0; + ncells_pml_domain_lmax = 0; + } + if (patch->isYmin() ) {//&& i_boundary_ == 2 ) { + ncells_pml_rmin = params.number_of_pml_cells[1][0]; + ncells_pml_domain_rmin = ncells_pml_rmin + 1*oversize[1] + nsolver/2; + domain_oversize_r = oversize[1] ; + } + else { + ncells_pml_rmin = 0; + ncells_pml_domain_rmin = 0; + } + if (patch->isYmax() ) {//&& i_boundary_ == 3 ) { + ncells_pml_rmax = params.number_of_pml_cells[1][1]; + ncells_pml_domain_rmax = ncells_pml_rmax + 1*oversize[1] + nsolver/2; + domain_oversize_r = oversize[1] ; + } + else { + ncells_pml_rmax = 0; + ncells_pml_domain_rmax = 0; + } + + ncells_pml = params.number_of_pml_cells[iDim][min_or_max]; + ncells_pml_domain = ncells_pml+1*oversize[iDim]+ nsolver/2; + + // Define min and max idx to exchange + // the good data f(solver,oversize) + if (min_or_max==0){ + // if min border : Exchange of data (for domain to pml-domain) + // min2exchange <= i < max2exchange + min2exchange = 1*nsolver/2 ; + max2exchange = 2*nsolver/2 ; + // Solver + solvermin = nsolver/2 ; + solvermax = ncells_pml_domain - oversize[iDim] ; + } + else if (min_or_max==1){ + // if max border : Exchange of data (for domain to pml-domain) + // min2exchange <= i < max2exchange + min2exchange = 1*nsolver/2 ; + max2exchange = 2*nsolver/2 ; + // Solver + solvermin = oversize[iDim] + nsolver/2 - nsolver/2 + 1 ; + solvermax = ncells_pml_domain-nsolver/2 ; + } + + if (ncells_pml==0){ + ERROR("PML domain have to be >0 cells in thickness"); + } + + std::vector dimPrim( params.nDim_field ); + for( int i=0 ; isetDomainSizeAndCoefficients( iDim, min_or_max, ncells_pml_domain, startpml, ncells_pml_min, ncells_pml_max, patch ); + + for ( unsigned int imode=0 ; imodeisXmin() ) { + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + for( unsigned int imode=0 ; imode( EMfields ) )->El_[imode]; + cField2D *Er_domain = ( static_cast( EMfields ) )->Er_[imode]; + cField2D *Et_domain = ( static_cast( EMfields ) )->Et_[imode]; + cField2D *Bl_domain = ( static_cast( EMfields ) )->Bl_[imode]; + cField2D *Br_domain = ( static_cast( EMfields ) )->Br_[imode]; + cField2D *Bt_domain = ( static_cast( EMfields ) )->Bt_[imode]; + + // 2. Exchange field PML <- Domain + for ( int i=min2exchange ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + for( unsigned int imode=0 ; imode( EMfields ) )->El_[imode]; + cField2D *Er_domain = ( static_cast( EMfields ) )->Er_[imode]; + cField2D *Et_domain = ( static_cast( EMfields ) )->Et_[imode]; + cField2D *Bl_domain = ( static_cast( EMfields ) )->Bl_[imode]; + cField2D *Br_domain = ( static_cast( EMfields ) )->Br_[imode]; + cField2D *Bt_domain = ( static_cast( EMfields ) )->Bt_[imode]; + // for Br^(d,p) + vector yp( 1 ); + for( unsigned int j=3*( patch->isYmin() ) ; jisYmax() ; j++ ) { + + std::complex byW = 0.; + yp[0] = patch->getDomainLocalMin( 1 ) +( (double)j - (double)EMfields->oversize[1] )*d[1]; + + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + if (vecLaser[ilaser]->spacetime.size() > 2){ + byW += vecLaser[ilaser]->getAmplitudecomplexN(yp, time_dual, 0, 0, 2*imode); + } else { + if( imode==1 ) { + byW += vecLaser[ilaser]->getAmplitude0( yp, time_dual, 1+2*j, 0 ) + + Icpx * vecLaser[ilaser]->getAmplitude1( yp, time_dual, 1+2*j, 0 ); + } + } + } + // unsigned int i=120; + // ( *Br_domain )( i, j ) += factor_laser_space_time*factor_laser_angle_W*byW; + unsigned int i=ncells_pml_domain-domain_oversize_l-nsolver/2; + ( *Hr_[imode] )( i, j ) += factor_laser_space_time*factor_laser_angle_W*byW; + ( *Br_[imode] )( i, j ) += factor_laser_space_time*factor_laser_angle_W*byW; + } + // for Bt^(d,d) + vector yd( 1 ); + for( unsigned int j=3*( patch->isYmin() ); jisYmax() ; j++ ) { + + std::complex bzW = 0.; + yd[0] = patch->getDomainLocalMin( 1 ) + ( (double)j - 0.5 - (double)EMfields->oversize[1] )*d[1]; + + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + if (vecLaser[ilaser]->spacetime.size() > 2){ + bzW += vecLaser[ilaser]->getAmplitudecomplexN(yd, time_dual, 0, 0, 2*imode+1); + } else { + if( imode==1 ) { + bzW += vecLaser[ilaser]->getAmplitude1( yd, time_dual, 2*j, 0 ) + - Icpx * vecLaser[ilaser]->getAmplitude0( yd, time_dual, 2*j, 0 ); + } + } + } + // unsigned int i=120; + // ( *Bt_domain )( i, j ) += factor_laser_space_time*factor_laser_angle_W*bzW; + unsigned int i=ncells_pml_domain-domain_oversize_l-nsolver/2; + ( *Ht_[imode] )( i, j ) += factor_laser_space_time*factor_laser_angle_W*bzW; + ( *Bt_[imode] )( i, j ) += factor_laser_space_time*factor_laser_angle_W*bzW; + } + //Redo condition on axis for Bt because it was modified + if( patch->isYmin() && imode != 1 ) { + unsigned int i=ncells_pml_domain-domain_oversize_l-nsolver/2; + ( *Bt_[imode] )( i, 2 ) = -( *Bt_[imode] )( i, 3 ); + ( *Ht_[imode] )( i, 2 ) = -( *Ht_[imode] )( i, 3 ); + } + if( patch->isYmin() && imode == 1 ) { + unsigned int i=ncells_pml_domain-domain_oversize_l-nsolver/2; + ( *Bt_[imode] )( i, 2 )= -2.*Icpx*( *Br_[imode] )( i, 2 )-( *Bt_[imode] )( i, 3 ); + ( *Ht_[imode] )( i, 2 )= -2.*Icpx*( *Hr_[imode] )( i, 2 )-( *Ht_[imode] )( i, 3 ); + } + } + + // 4. Exchange PML -> Domain + for( unsigned int imode=0 ; imode( EMfields ) )->El_[imode]; + cField2D *Er_domain = ( static_cast( EMfields ) )->Er_[imode]; + cField2D *Et_domain = ( static_cast( EMfields ) )->Et_[imode]; + cField2D *Bl_domain = ( static_cast( EMfields ) )->Bl_[imode]; + cField2D *Br_domain = ( static_cast( EMfields ) )->Br_[imode]; + cField2D *Bt_domain = ( static_cast( EMfields ) )->Bt_[imode]; + // Primals in x-direction + for (int i=0 ; i < nsolver/2 ; i++){ + for ( int j=0 ; jisXmax() ) { + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + for( unsigned int imode=0 ; imode( EMfields ) )->El_[imode]; + cField2D *Er_domain = ( static_cast( EMfields ) )->Er_[imode]; + cField2D *Et_domain = ( static_cast( EMfields ) )->Et_[imode]; + cField2D *Bl_domain = ( static_cast( EMfields ) )->Bl_[imode]; + cField2D *Br_domain = ( static_cast( EMfields ) )->Br_[imode]; + cField2D *Bt_domain = ( static_cast( EMfields ) )->Bt_[imode]; + + for ( int i=min2exchange ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + //Injecting a laser + for( unsigned int imode=0 ; imode( EMfields ) )->El_[imode]; + cField2D *Er_domain = ( static_cast( EMfields ) )->Er_[imode]; + cField2D *Et_domain = ( static_cast( EMfields ) )->Et_[imode]; + cField2D *Bl_domain = ( static_cast( EMfields ) )->Bl_[imode]; + cField2D *Br_domain = ( static_cast( EMfields ) )->Br_[imode]; + cField2D *Bt_domain = ( static_cast( EMfields ) )->Bt_[imode]; + // for Br^(d,p) + vector yp( 1 ); + for( unsigned int j=3*( patch->isYmin() ) ; j byE = 0.; + yp[0] = patch->getDomainLocalMin( 1 ) +( (double)j - (double)EMfields->oversize[1] )*d[1]; + + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + if (vecLaser[ilaser]->spacetime.size() > 2){ + byE += vecLaser[ilaser]->getAmplitudecomplexN(yp, time_dual, 0, 0, 2*imode); + } else { + if( imode==1 ) { + byE += vecLaser[ilaser]->getAmplitude0( yp, time_dual, 1+2*j, 0 ) + + Icpx * vecLaser[ilaser]->getAmplitude1( yp, time_dual, 1+2*j, 0 ); + } + } + } + // unsigned int i=120; + // ( *Br_domain )( i, j ) += factor_laser_space_time*factor_laser_angle_E*byE; + unsigned int i=domain_oversize_l+nsolver/2; + ( *Hr_[imode] )( i, j ) += factor_laser_space_time*factor_laser_angle_E*byE; + ( *Br_[imode] )( i, j ) += factor_laser_space_time*factor_laser_angle_E*byE; + } + // for Bt^(d,d) + vector yd( 1 ); + for( unsigned int j=3*( patch->isYmin() ); j bzE = 0.; + yd[0] = patch->getDomainLocalMin( 1 ) + ( (double)j - 0.5 - (double)EMfields->oversize[1] )*d[1]; + + // Lasers + for( unsigned int ilaser=0; ilaser< vecLaser.size(); ilaser++ ) { + if (vecLaser[ilaser]->spacetime.size() > 2){ + bzE += vecLaser[ilaser]->getAmplitudecomplexN(yd, time_dual, 0, 0, 2*imode+1); + } else { + if( imode==1 ) { + bzE += vecLaser[ilaser]->getAmplitude1( yd, time_dual, 2*j, 0 ) + - Icpx * vecLaser[ilaser]->getAmplitude0( yd, time_dual, 2*j, 0 ); + } + } + } + // unsigned int i=120; + // ( *Bt_domain )( i, j ) += factor_laser_space_time*factor_laser_angle_E*bzE; + unsigned int i=domain_oversize_l+nsolver/2; + ( *Ht_[imode] )( i, j ) += factor_laser_space_time*factor_laser_angle_E*bzE; + ( *Bt_[imode] )( i, j ) += factor_laser_space_time*factor_laser_angle_E*bzE; + } + //Redo condition on axis for Bt because it was modified + if( patch->isYmin() && imode != 1 ) { + unsigned int i=domain_oversize_l+nsolver/2; + ( *Bt_[imode] )( i, 2 ) = -( *Bt_[imode] )( i, 3 ); + ( *Ht_[imode] )( i, 2 ) = -( *Ht_[imode] )( i, 3 ); + } + if( patch->isYmin() && imode == 1 ) { + unsigned int i=domain_oversize_l+nsolver/2; + ( *Bt_[imode] )( i, 2 )= -2.*Icpx*( *Br_[imode] )( i, 2 )-( *Bt_[imode] )( i, 3 ); + ( *Ht_[imode] )( i, 2 )= -2.*Icpx*( *Hr_[imode] )( i, 2 )-( *Ht_[imode] )( i, 3 ); + } + } + + // 4. Exchange Domain -> PML + for( unsigned int imode=0 ; imode( EMfields ) )->El_[imode]; + cField2D *Er_domain = ( static_cast( EMfields ) )->Er_[imode]; + cField2D *Et_domain = ( static_cast( EMfields ) )->Et_[imode]; + cField2D *Bl_domain = ( static_cast( EMfields ) )->Bl_[imode]; + cField2D *Br_domain = ( static_cast( EMfields ) )->Br_[imode]; + cField2D *Bt_domain = ( static_cast( EMfields ) )->Bt_[imode]; + // Primals in x-direction + for (int i=0 ; i < nsolver/2-1 ; i++){ + for ( int j=0 ; jisYmin() ) { + // ERROR("PML not allow on the symetric axis") + + ElectroMagnBCAM_PML* pml_fields_lmin = NULL ; + ElectroMagnBCAM_PML* pml_fields_lmax = NULL ; + + if(ncells_pml_lmin != 0){ + pml_fields_lmin = static_cast( EMfields->emBoundCond[0] ); + } + if(ncells_pml_lmax != 0){ + pml_fields_lmax = static_cast( EMfields->emBoundCond[1] ); + } + + cField2D* El_pml_lmin = NULL; + cField2D* Er_pml_lmin = NULL; + cField2D* Et_pml_lmin = NULL; + cField2D* Hl_pml_lmin = NULL; + cField2D* Hr_pml_lmin = NULL; + cField2D* Ht_pml_lmin = NULL; + cField2D* Dl_pml_lmin = NULL; + cField2D* Dr_pml_lmin = NULL; + cField2D* Dt_pml_lmin = NULL; + cField2D* Bl_pml_lmin = NULL; + cField2D* Br_pml_lmin = NULL; + cField2D* Bt_pml_lmin = NULL; + + cField2D* El_pml_lmax = NULL; + cField2D* Er_pml_lmax = NULL; + cField2D* Et_pml_lmax = NULL; + cField2D* Hl_pml_lmax = NULL; + cField2D* Hr_pml_lmax = NULL; + cField2D* Ht_pml_lmax = NULL; + cField2D* Dl_pml_lmax = NULL; + cField2D* Dr_pml_lmax = NULL; + cField2D* Dt_pml_lmax = NULL; + cField2D* Bl_pml_lmax = NULL; + cField2D* Br_pml_lmax = NULL; + cField2D* Bt_pml_lmax = NULL; + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + for( unsigned int imode=0 ; imode( EMfields ) )->El_[imode]; + cField2D *Er_domain = ( static_cast( EMfields ) )->Er_[imode]; + cField2D *Et_domain = ( static_cast( EMfields ) )->Et_[imode]; + cField2D *Bl_domain = ( static_cast( EMfields ) )->Bl_[imode]; + cField2D *Br_domain = ( static_cast( EMfields ) )->Br_[imode]; + cField2D *Bt_domain = ( static_cast( EMfields ) )->Bt_[imode]; + + if(ncells_pml_lmin != 0){ + El_pml_lmin = pml_fields_lmin->El_[imode]; + Er_pml_lmin = pml_fields_lmin->Er_[imode]; + Et_pml_lmin = pml_fields_lmin->Et_[imode]; + Hl_pml_lmin = pml_fields_lmin->Hl_[imode]; + Hr_pml_lmin = pml_fields_lmin->Hr_[imode]; + Ht_pml_lmin = pml_fields_lmin->Ht_[imode]; + Dl_pml_lmin = pml_fields_lmin->Dl_[imode]; + Dr_pml_lmin = pml_fields_lmin->Dr_[imode]; + Dt_pml_lmin = pml_fields_lmin->Dt_[imode]; + Bl_pml_lmin = pml_fields_lmin->Bl_[imode]; + Br_pml_lmin = pml_fields_lmin->Br_[imode]; + Bt_pml_lmin = pml_fields_lmin->Bt_[imode]; + } + + if(ncells_pml_lmax != 0){ + El_pml_lmax = pml_fields_lmax->El_[imode]; + Er_pml_lmax = pml_fields_lmax->Er_[imode]; + Et_pml_lmax = pml_fields_lmax->Et_[imode]; + Hl_pml_lmax = pml_fields_lmax->Hl_[imode]; + Hr_pml_lmax = pml_fields_lmax->Hr_[imode]; + Ht_pml_lmax = pml_fields_lmax->Ht_[imode]; + Dl_pml_lmax = pml_fields_lmax->Dl_[imode]; + Dr_pml_lmax = pml_fields_lmax->Dr_[imode]; + Dt_pml_lmax = pml_fields_lmax->Dt_[imode]; + Bl_pml_lmax = pml_fields_lmax->Bl_[imode]; + Br_pml_lmax = pml_fields_lmax->Br_[imode]; + Bt_pml_lmax = pml_fields_lmax->Bt_[imode]; + } + + // 2. Exchange field PML <- Domain + for ( int j=min2exchange ; jisXmin()) { + if(ncells_pml_lmin != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_lmax != 0){ + for ( int i=0 ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + for( unsigned int imode=0 ; imode( EMfields ) )->El_[imode]; + cField2D *Er_domain = ( static_cast( EMfields ) )->Er_[imode]; + cField2D *Et_domain = ( static_cast( EMfields ) )->Et_[imode]; + cField2D *Bl_domain = ( static_cast( EMfields ) )->Bl_[imode]; + cField2D *Br_domain = ( static_cast( EMfields ) )->Br_[imode]; + cField2D *Bt_domain = ( static_cast( EMfields ) )->Bt_[imode]; + + if (patch->isXmin()) { + if(ncells_pml_lmin != 0){ + El_pml_lmin = pml_fields_lmin->El_[imode]; + Er_pml_lmin = pml_fields_lmin->Er_[imode]; + Et_pml_lmin = pml_fields_lmin->Et_[imode]; + Hl_pml_lmin = pml_fields_lmin->Hl_[imode]; + Hr_pml_lmin = pml_fields_lmin->Hr_[imode]; + Ht_pml_lmin = pml_fields_lmin->Ht_[imode]; + Dl_pml_lmin = pml_fields_lmin->Dl_[imode]; + Dr_pml_lmin = pml_fields_lmin->Dr_[imode]; + Dt_pml_lmin = pml_fields_lmin->Dt_[imode]; + Bl_pml_lmin = pml_fields_lmin->Bl_[imode]; + Br_pml_lmin = pml_fields_lmin->Br_[imode]; + Bt_pml_lmin = pml_fields_lmin->Bt_[imode]; + } + } + + if (patch->isXmax()) { + if(ncells_pml_lmax != 0){ + El_pml_lmax = pml_fields_lmax->El_[imode]; + Er_pml_lmax = pml_fields_lmax->Er_[imode]; + Et_pml_lmax = pml_fields_lmax->Et_[imode]; + Hl_pml_lmax = pml_fields_lmax->Hl_[imode]; + Hr_pml_lmax = pml_fields_lmax->Hr_[imode]; + Ht_pml_lmax = pml_fields_lmax->Ht_[imode]; + Dl_pml_lmax = pml_fields_lmax->Dl_[imode]; + Dr_pml_lmax = pml_fields_lmax->Dr_[imode]; + Dt_pml_lmax = pml_fields_lmax->Dt_[imode]; + Bl_pml_lmax = pml_fields_lmax->Bl_[imode]; + Br_pml_lmax = pml_fields_lmax->Br_[imode]; + Bt_pml_lmax = pml_fields_lmax->Bt_[imode]; + } + } + + // 4. Exchange PML -> Domain + // Primals in y-direction + for (int j=0 ; j < nsolver/2 ; j++){ + for ( int i=0 ; i PML x MIN + // Primal in y-direction + for (int j=0 ; j < nsolver/2 ; j++){ + if (patch->isXmin()) { + if(ncells_pml_lmin != 0){ + for ( int i=0 ; iisXmin()) { + if(ncells_pml_lmin != 0){ + for ( int i=0 ; i PML x MAX + // Primal in y-direction + for (int j=0 ; j < nsolver/2 ; j++){ + if (patch->isXmax()) { + if(ncells_pml_lmax != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_lmax != 0){ + for ( int i=0 ; iisYmax() ) { + + ElectroMagnBCAM_PML* pml_fields_lmin = NULL ; + ElectroMagnBCAM_PML* pml_fields_lmax = NULL ; + + if(ncells_pml_lmin != 0){ + pml_fields_lmin = static_cast( EMfields->emBoundCond[0] ); + } + if(ncells_pml_lmax != 0){ + pml_fields_lmax = static_cast( EMfields->emBoundCond[1] ); + } + + cField2D* El_pml_lmin = NULL; + cField2D* Er_pml_lmin = NULL; + cField2D* Et_pml_lmin = NULL; + cField2D* Hl_pml_lmin = NULL; + cField2D* Hr_pml_lmin = NULL; + cField2D* Ht_pml_lmin = NULL; + cField2D* Dl_pml_lmin = NULL; + cField2D* Dr_pml_lmin = NULL; + cField2D* Dt_pml_lmin = NULL; + cField2D* Bl_pml_lmin = NULL; + cField2D* Br_pml_lmin = NULL; + cField2D* Bt_pml_lmin = NULL; + + cField2D* El_pml_lmax = NULL; + cField2D* Er_pml_lmax = NULL; + cField2D* Et_pml_lmax = NULL; + cField2D* Hl_pml_lmax = NULL; + cField2D* Hr_pml_lmax = NULL; + cField2D* Ht_pml_lmax = NULL; + cField2D* Dl_pml_lmax = NULL; + cField2D* Dr_pml_lmax = NULL; + cField2D* Dt_pml_lmax = NULL; + cField2D* Bl_pml_lmax = NULL; + cField2D* Br_pml_lmax = NULL; + cField2D* Bt_pml_lmax = NULL; + + // 1. Solve Maxwell_PML for E-field : + // As if B-field isn't updated + pml_solver_->compute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + //pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + for( unsigned int imode=0 ; imode( EMfields ) )->El_[imode]; + cField2D *Er_domain = ( static_cast( EMfields ) )->Er_[imode]; + cField2D *Et_domain = ( static_cast( EMfields ) )->Et_[imode]; + cField2D *Bl_domain = ( static_cast( EMfields ) )->Bl_[imode]; + cField2D *Br_domain = ( static_cast( EMfields ) )->Br_[imode]; + cField2D *Bt_domain = ( static_cast( EMfields ) )->Bt_[imode]; + + if(ncells_pml_lmin != 0){ + El_pml_lmin = pml_fields_lmin->El_[imode]; + Er_pml_lmin = pml_fields_lmin->Er_[imode]; + Et_pml_lmin = pml_fields_lmin->Et_[imode]; + Hl_pml_lmin = pml_fields_lmin->Hl_[imode]; + Hr_pml_lmin = pml_fields_lmin->Hr_[imode]; + Ht_pml_lmin = pml_fields_lmin->Ht_[imode]; + Dl_pml_lmin = pml_fields_lmin->Dl_[imode]; + Dr_pml_lmin = pml_fields_lmin->Dr_[imode]; + Dt_pml_lmin = pml_fields_lmin->Dt_[imode]; + Bl_pml_lmin = pml_fields_lmin->Bl_[imode]; + Br_pml_lmin = pml_fields_lmin->Br_[imode]; + Bt_pml_lmin = pml_fields_lmin->Bt_[imode]; + } + + if(ncells_pml_lmax != 0){ + El_pml_lmax = pml_fields_lmax->El_[imode]; + Er_pml_lmax = pml_fields_lmax->Er_[imode]; + Et_pml_lmax = pml_fields_lmax->Et_[imode]; + Hl_pml_lmax = pml_fields_lmax->Hl_[imode]; + Hr_pml_lmax = pml_fields_lmax->Hr_[imode]; + Ht_pml_lmax = pml_fields_lmax->Ht_[imode]; + Dl_pml_lmax = pml_fields_lmax->Dl_[imode]; + Dr_pml_lmax = pml_fields_lmax->Dr_[imode]; + Dt_pml_lmax = pml_fields_lmax->Dt_[imode]; + Bl_pml_lmax = pml_fields_lmax->Bl_[imode]; + Br_pml_lmax = pml_fields_lmax->Br_[imode]; + Bt_pml_lmax = pml_fields_lmax->Bt_[imode]; + } + + for ( int j=min2exchange ; jisXmin()) { + if(ncells_pml_lmin != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_lmax != 0){ + for ( int i=0 ; icompute_E_from_D( EMfields, iDim, min_or_max, solvermin, solvermax); + pml_solver_->compute_H_from_B( EMfields, iDim, min_or_max, solvermin, solvermax); + + for( unsigned int imode=0 ; imode( EMfields ) )->El_[imode]; + cField2D *Er_domain = ( static_cast( EMfields ) )->Er_[imode]; + cField2D *Et_domain = ( static_cast( EMfields ) )->Et_[imode]; + cField2D *Bl_domain = ( static_cast( EMfields ) )->Bl_[imode]; + cField2D *Br_domain = ( static_cast( EMfields ) )->Br_[imode]; + cField2D *Bt_domain = ( static_cast( EMfields ) )->Bt_[imode]; + + if (patch->isXmin()) { + if(ncells_pml_lmin != 0){ + El_pml_lmin = pml_fields_lmin->El_[imode]; + Er_pml_lmin = pml_fields_lmin->Er_[imode]; + Et_pml_lmin = pml_fields_lmin->Et_[imode]; + Hl_pml_lmin = pml_fields_lmin->Hl_[imode]; + Hr_pml_lmin = pml_fields_lmin->Hr_[imode]; + Ht_pml_lmin = pml_fields_lmin->Ht_[imode]; + Dl_pml_lmin = pml_fields_lmin->Dl_[imode]; + Dr_pml_lmin = pml_fields_lmin->Dr_[imode]; + Dt_pml_lmin = pml_fields_lmin->Dt_[imode]; + Bl_pml_lmin = pml_fields_lmin->Bl_[imode]; + Br_pml_lmin = pml_fields_lmin->Br_[imode]; + Bt_pml_lmin = pml_fields_lmin->Bt_[imode]; + } + } + + if (patch->isXmax()) { + if(ncells_pml_lmax != 0){ + El_pml_lmax = pml_fields_lmax->El_[imode]; + Er_pml_lmax = pml_fields_lmax->Er_[imode]; + Et_pml_lmax = pml_fields_lmax->Et_[imode]; + Hl_pml_lmax = pml_fields_lmax->Hl_[imode]; + Hr_pml_lmax = pml_fields_lmax->Hr_[imode]; + Ht_pml_lmax = pml_fields_lmax->Ht_[imode]; + Dl_pml_lmax = pml_fields_lmax->Dl_[imode]; + Dr_pml_lmax = pml_fields_lmax->Dr_[imode]; + Dt_pml_lmax = pml_fields_lmax->Dt_[imode]; + Bl_pml_lmax = pml_fields_lmax->Bl_[imode]; + Br_pml_lmax = pml_fields_lmax->Br_[imode]; + Bt_pml_lmax = pml_fields_lmax->Bt_[imode]; + } + } + + // 4. Exchange PML -> Domain + // Primals in y-direction + for (int j=0 ; j < nsolver/2-1 ; j++){ + for ( int i=0 ; i PML x MIN + // Primal in y-direction + for (int j=0 ; j < nsolver/2-1 ; j++){ + if (patch->isXmin()) { + if(ncells_pml_lmin != 0){ + for ( int i=0 ; iisXmin()) { + if(ncells_pml_lmin != 0){ + for ( int i=0 ; i PML x MAX + // Primal in y-direction + for (int j=0 ; j < nsolver/2-1 ; j++){ + if (patch->isXmax()) { + if(ncells_pml_lmax != 0){ + for ( int i=0 ; iisXmax()) { + if(ncells_pml_lmax != 0){ + for ( int i=0 ; i +#include +#include "Tools.h" +#include "ElectroMagnBCAM.h" +#include "ElectroMagnAM.h" +#include "cField2D.h" +#include "dcomplex.h" + +class Params; +class ElectroMagn; +class Field; +class PML_SolverAM; + +class ElectroMagnBCAM_PML : public ElectroMagnBCAM +{ +public: + + ElectroMagnBCAM_PML( Params ¶ms, Patch *patch, unsigned int _min_max ); + ~ElectroMagnBCAM_PML(); + + virtual void apply( ElectroMagn *EMfields, double time_dual, Patch *patch ) override; + + void save_fields( Field *, Patch *patch ) override; + void disableExternalFields() override; + + std::vector El_ ;//= NULL; + std::vector Er_ ;//= NULL; + std::vector Et_ ;//= NULL; + std::vector Bl_ ;//= NULL; + std::vector Br_ ;//= NULL; + std::vector Bt_ ;//= NULL; + std::vector Dl_ ;//= NULL; + std::vector Dr_ ;//= NULL; + std::vector Dt_ ;//= NULL; + std::vector Hl_ ;//= NULL; + std::vector Hr_ ;//= NULL; + std::vector Ht_ ;//= NULL; + + int domain_oversize_l; + int domain_oversize_r; + + int ncells_pml_lmin ; + int ncells_pml_lmax ; + int ncells_pml_rmin ; + int ncells_pml_rmax ; + + int ncells_pml_domain_lmin; + int ncells_pml_domain_lmax; + int ncells_pml_domain_rmin; + int ncells_pml_domain_rmax; + + int ncells_pml; + int ncells_pml_domain; + + int nsolver; + + int min2exchange; + int max2exchange; + int solvermin; + int solvermax; + + int rpml_size_in_l; + int startpml; + + double factor_laser_space_time ; + double factor_laser_angle_W ; + double factor_laser_angle_E ; + +}; + +#endif diff --git a/src/ElectroMagnBC/ElectroMagnBCAM_SM.cpp b/src/ElectroMagnBC/ElectroMagnBCAM_SM.cpp index 1216aa12b..28f4b42b7 100755 --- a/src/ElectroMagnBC/ElectroMagnBCAM_SM.cpp +++ b/src/ElectroMagnBC/ElectroMagnBCAM_SM.cpp @@ -168,7 +168,7 @@ void ElectroMagnBCAM_SM::disableExternalFields() // --------------------------------------------------------------------------------------------------------------------- void ElectroMagnBCAM_SM::apply( ElectroMagn *EMfields, double time_dual, Patch *patch ) { - + //ERROR( "Test SM" ) // Loop on imode for( unsigned int imode=0 ; imode> dx,dy then solver become like the 2d solver + //if( (dx!=dy)||(dx!=dz)||(dy!=dz) ) { + // ERROR( "Bouchard solver requires the same cell-length in x, y and z directions" ); + //} if( dx_ov_dt!=2 ) { WARNING( "Bouchard solver requires dx/dt = 2 (Magical Timestep)" ); } @@ -166,8 +173,7 @@ void MF_Solver3D_Bouchard::operator()( ElectroMagn *fields ) ( *Bz3D )( nx_d-2, j, k ) += dt_ov_dx * ( ( *Ey3D )( nx_d-3, j, k ) - ( *Ey3D )( nx_d-2, j, k ) ) + dt_ov_dy * ( ( *Ex3D )( nx_d-2, j, k ) - ( *Ex3D )( nx_d-2, j-1, k ) ); } - } - + } } if( EM3D->isYmin ) { diff --git a/src/ElectroMagnSolver/PML_Solver2D_Bouchard.cpp b/src/ElectroMagnSolver/PML_Solver2D_Bouchard.cpp new file mode 100644 index 000000000..07a9d9f09 --- /dev/null +++ b/src/ElectroMagnSolver/PML_Solver2D_Bouchard.cpp @@ -0,0 +1,732 @@ +#include "PML_Solver2D_Bouchard.h" +#include "ElectroMagn.h" +#include "ElectroMagnBC2D_PML.h" +#include "Field2D.h" +#include "Patch.h" + +#include + +PML_Solver2D_Bouchard::PML_Solver2D_Bouchard(Params ¶ms) + : Solver2D(params) +{ + //ERROR("Under development, not yet working"); + double dt = params.timestep; + dx = params.cell_length[0]; + dy = params.cell_length[1]; + double dx_ov_dt = dx/dt; + double dy_ov_dt = dy/dt; + double dt_ov_dx = dt/dx; + double dt_ov_dy = dt/dy; + if( dx!=dy ) { + ERROR( "Bouchard solver requires the same cell-length in x and y directions" ); + } + if( dx_ov_dt!=2 ) { + WARNING( "Bouchard solver requires dx/dt = 2 (Magical Timestep)" ); + } + + // On the axes v_phi^max = 1.01c and is below c @ 0.54 kxdx/pi + // So there could existe a numerical cherenkov emission at this point + // On the diagonal v_phi^max = 1.01c and is below c @ 0.85 sqrt((kxdx)^2+(kydy)^2) + double delta = 0.1222*(1-pow(2.,2))/4. ; + double beta = -0.1727*(1-0.5*pow(2.,2)-4.*delta)/4. ; + double alpha = 1-2.*beta-3.*delta; + + beta_xy = beta; + beta_yx = beta; + delta_y = delta; + delta_x = delta; + + alpha_y = alpha; + alpha_x = alpha; + + // Ax = alpha_x*dt/dx; + // Ay = alpha_y*dt/dy; + // Bx = beta_xy*dt/dx; + // By = beta_yx*dt/dy; + // Dx = delta_x*dt/dy; + // Dy = delta_y*dt/dy; + + Ax = alpha_x/dx; + Ay = alpha_y/dy; + Bx = beta_xy/dx; + By = beta_yx/dy; + Dx = delta_x/dy; + Dy = delta_y/dy; + + //Define here the value of coefficient kappa_x_max, power_kappa_x, sigma_x_max, power_sigma_x + sigma_x_max = 20.; + kappa_x_max = 80.; + sigma_power_pml_x = 2.; + kappa_power_pml_x = 4.; + //Define here the value of coefficient kappa_y_max, power_kappa_y, sigma_y_max, power_sigma_y + sigma_y_max = 20.; + kappa_y_max = 80.; + sigma_power_pml_y = 2.; + kappa_power_pml_y = 4.; +} + +PML_Solver2D_Bouchard::~PML_Solver2D_Bouchard() +{ +} + +void PML_Solver2D_Bouchard::operator() ( ElectroMagn* fields ) +{ + ERROR( "This is not a solver for the main domain" ); +} + +void PML_Solver2D_Bouchard::setDomainSizeAndCoefficients( int iDim, int min_or_max, int ncells_pml_domain, int startpml, int* ncells_pml_min, int* ncells_pml_max, Patch* patch ) +{ + if ( iDim == 0 ) { + nx_p = ncells_pml_domain; + nx_d = ncells_pml_domain+1; + } + else if ( iDim == 1 ) { + ny_p = ncells_pml_domain; + ny_d = ncells_pml_domain+1; + nx_p += ncells_pml_min[0] + ncells_pml_max[0]; + nx_d += ncells_pml_min[0] + ncells_pml_max[0]; + } + + //PML Coeffs Kappa,Sigma ... + //Primal + kappa_x_p.resize( nx_p ); + sigma_x_p.resize( nx_p ); + kappa_y_p.resize( ny_p ); + sigma_y_p.resize( ny_p ); + kappa_z_p.resize( 1 ); + sigma_z_p.resize( 1 ); + //Dual + kappa_x_d.resize( nx_d ); + sigma_x_d.resize( nx_d ); + kappa_y_d.resize( ny_d ); + sigma_y_d.resize( ny_d ); + kappa_z_d.resize( 1 ); + sigma_z_d.resize( 1 ); + + // While min and max of each dim have got the same ncells_pml, coefficient are the same + // for min and max PML + + //Coeffs for El,Dl,Hl,Bl + //Primal + c1_p_xfield.resize( ny_p ); // j-dependent + c2_p_xfield.resize( ny_p ); // j-dependent + c3_p_xfield.resize( 1 ); // k-dependent + c4_p_xfield.resize( 1 ); // k-dependent + c5_p_xfield.resize( nx_p ); // i-dependent + c6_p_xfield.resize( nx_p ); // i-dependent + //Dual + c1_d_xfield.resize( ny_d ); // j-dependent + c2_d_xfield.resize( ny_d ); // j-dependent + c3_d_xfield.resize( 1 ); // j-dependent + c4_d_xfield.resize( 1 ); // j-dependent + c5_d_xfield.resize( nx_d ); // i-dependent + c6_d_xfield.resize( nx_d ); // i-dependent + //Coefs for Er,Dr,Hr,Br + //Primal + c1_p_yfield.resize( 1 ); // k-dependent + c2_p_yfield.resize( 1 ); // k-dependent + c3_p_yfield.resize( nx_p ); // i-dependent + c4_p_yfield.resize( nx_p ); // i-dependent + c5_p_yfield.resize( ny_p ); // j-dependent + c6_p_yfield.resize( ny_p ); // j-dependent + //Dual + c1_d_yfield.resize( 1 ); // k-dependent + c2_d_yfield.resize( 1 ); // k-dependent + c3_d_yfield.resize( nx_d ); // i-dependent + c4_d_yfield.resize( nx_d ); // i-dependent + c5_d_yfield.resize( ny_d ); // k-dependent + c6_d_yfield.resize( ny_d ); // k-dependent + //Coefs for Et,Dt,Ht,Bt + //Primal + c1_p_zfield.resize( nx_p ); // j-dependent + c2_p_zfield.resize( nx_p ); // j-dependent + c3_p_zfield.resize( ny_p ); // i-dependent + c4_p_zfield.resize( ny_p ); // i-dependent + c5_p_zfield.resize( 1 ); // j-dependent + c6_p_zfield.resize( 1 ); // j-dependent + //Dual + c1_d_zfield.resize( nx_d ); // i-dependent + c2_d_zfield.resize( nx_d ); // i-dependent + c3_d_zfield.resize( ny_d ); // j-dependent + c4_d_zfield.resize( ny_d ); // j-dependent + c5_d_zfield.resize( 1 ); // k-dependent + c6_d_zfield.resize( 1 ); // k-dependent + + // Quote of the first primal grid-point where PML solver will be apply + // The first dual grid-point is at getDomainLocalMax( 0 ) - 0.5*dx and getDomainLocalMax( 1 ) - 0.5*dr + // xmax = patch->getDomainLocalMax( 0 ); + // ymax = patch->getDomainLocalMax( 1 ); + + if ( iDim == 0 ) { + // 3 cells (oversize) are vaccum so the PML media begin at r0 which is : + // Eventually the size of PML media is : + length_y_pml = 0 ; + length_x_pml = (ncells_pml_domain-startpml+0.5)*dx ; + // Primal grid + // X-direction + // Params for first cell of PML-patch (vacuum) i = 0,1,2 + for ( int i=0 ; i=3 + for ( int i=startpml; i=4 + for ( int i=startpml+1 ; i=3 + for ( int j=startpml ; j=4 + for ( int j=startpml+1 ; j( fields->emBoundCond[iDim*2+min_or_max] ); + Field2D* Ex_pml = NULL; + Field2D* Ey_pml = NULL; + Field2D* Ez_pml = NULL; + Field2D* Hx_pml = NULL; + Field2D* Hy_pml = NULL; + Field2D* Hz_pml = NULL; + Field2D* Dx_pml = NULL; + Field2D* Dy_pml = NULL; + Field2D* Dz_pml = NULL; + + Ex_pml = pml_fields->Ex_; + Ey_pml = pml_fields->Ey_; + Ez_pml = pml_fields->Ez_; + Hx_pml = pml_fields->Hx_; + Hy_pml = pml_fields->Hy_; + Hz_pml = pml_fields->Hz_; + Dx_pml = pml_fields->Dx_; + Dy_pml = pml_fields->Dy_; + Dz_pml = pml_fields->Dz_; + + if (min_or_max==0){ + isMin=true; + isMax=false; + } + else if (min_or_max==1){ + isMin=false; + isMax=true; + } + + if (iDim == 0) { + //Electric field Ex^(d,p,p) Remind that in PML, there no current + for( unsigned int k=0 ; k<1 ; k++ ) { + for( unsigned int i=solvermin ; i( fields->emBoundCond[iDim*2+min_or_max] ); + Field2D* Ex_pml = NULL; + Field2D* Ey_pml = NULL; + Field2D* Ez_pml = NULL; + Field2D* Hx_pml = NULL; + Field2D* Hy_pml = NULL; + Field2D* Hz_pml = NULL; + Field2D* Bx_pml = NULL; + Field2D* By_pml = NULL; + Field2D* Bz_pml = NULL; + + Ex_pml = pml_fields->Ex_; + Ey_pml = pml_fields->Ey_; + Ez_pml = pml_fields->Ez_; + Hx_pml = pml_fields->Hx_; + Hy_pml = pml_fields->Hy_; + Hz_pml = pml_fields->Hz_; + Bx_pml = pml_fields->Bx_; + By_pml = pml_fields->By_; + Bz_pml = pml_fields->Bz_; + + if (min_or_max==0){ + isMin=true; + isMax=false; + } + else if (min_or_max==1){ + isMin=false; + isMax=true; + } + + if (iDim==0){ + //Magnetic field Bx^(p,d,d) Remind that in PML, there no current + for( unsigned int k=0 ; k<1 ; k++ ) { + for( unsigned int i=solvermin ; i kappa_x_p; + std::vector sigma_x_p; + std::vector kappa_x_d; + std::vector sigma_x_d; + std::vector kappa_y_p; + std::vector sigma_y_p; + std::vector kappa_y_d; + std::vector sigma_y_d; + std::vector kappa_z_p; + std::vector sigma_z_p; + std::vector kappa_z_d; + std::vector sigma_z_d; + + std::vector c1_p_xfield; + std::vector c2_p_xfield; + std::vector c3_p_xfield; + std::vector c4_p_xfield; + std::vector c5_p_xfield; + std::vector c6_p_xfield; + std::vector c1_d_xfield; + std::vector c2_d_xfield; + std::vector c3_d_xfield; + std::vector c4_d_xfield; + std::vector c5_d_xfield; + std::vector c6_d_xfield; + std::vector c1_p_yfield; + std::vector c2_p_yfield; + std::vector c3_p_yfield; + std::vector c4_p_yfield; + std::vector c5_p_yfield; + std::vector c6_p_yfield; + std::vector c1_d_yfield; + std::vector c2_d_yfield; + std::vector c3_d_yfield; + std::vector c4_d_yfield; + std::vector c5_d_yfield; + std::vector c6_d_yfield; + std::vector c1_p_zfield; + std::vector c2_p_zfield; + std::vector c3_p_zfield; + std::vector c4_p_zfield; + std::vector c5_p_zfield; + std::vector c6_p_zfield; + std::vector c1_d_zfield; + std::vector c2_d_zfield; + std::vector c3_d_zfield; + std::vector c4_d_zfield; + std::vector c5_d_zfield; + std::vector c6_d_zfield; + + //double xmin; + //double ymin; + + double length_x_pml; + double length_y_pml; + double length_x_pml_xmin; + double length_x_pml_xmax; + + bool isMin ; + bool isMax ; + double Dx_pml_old ; + double Dy_pml_old ; + double Dz_pml_old ; + double Bx_pml_old ; + double By_pml_old ; + double Bz_pml_old ; + +};//END class + +#endif + diff --git a/src/ElectroMagnSolver/PML_Solver2D_Yee.cpp b/src/ElectroMagnSolver/PML_Solver2D_Yee.cpp new file mode 100755 index 000000000..ab53beb62 --- /dev/null +++ b/src/ElectroMagnSolver/PML_Solver2D_Yee.cpp @@ -0,0 +1,649 @@ +#include "PML_Solver2D_Yee.h" +#include "ElectroMagn.h" +#include "ElectroMagnBC2D_PML.h" +#include "Field2D.h" +#include "Patch.h" + +PML_Solver2D_Yee::PML_Solver2D_Yee( Params ¶ms ) + : Solver2D( params ) +{ + //Define here the value of coefficient kappa_x_max, power_kappa_x, sigma_x_max, power_sigma_x + sigma_x_max = 20.; + kappa_x_max = 80.; + sigma_power_pml_x = 2.; + kappa_power_pml_x = 4.; + //Define here the value of coefficient kappa_y_max, power_kappa_y, sigma_y_max, power_sigma_y + sigma_y_max = 20.; + kappa_y_max = 80.; + sigma_power_pml_y = 2.; + kappa_power_pml_y = 4.; +} + +PML_Solver2D_Yee::~PML_Solver2D_Yee() +{ +} + +void PML_Solver2D_Yee::operator()( ElectroMagn *fields ) +{ + ERROR( "This is not a solver for the main domain" ); + + Field2D *Ex2D = static_cast( fields->Ex_ ); + Field2D *Ey2D = static_cast( fields->Ey_ ); + Field2D *Ez2D = static_cast( fields->Ez_ ); + Field2D *Bx2D = static_cast( fields->Bx_ ); + Field2D *By2D = static_cast( fields->By_ ); + Field2D *Bz2D = static_cast( fields->Bz_ ); + +} + +void PML_Solver2D_Yee::setDomainSizeAndCoefficients( int iDim, int min_or_max, int ncells_pml_domain, int startpml, int* ncells_pml_min, int* ncells_pml_max, Patch* patch ) +{ + if ( iDim == 0 ) { + nx_p = ncells_pml_domain; + nx_d = ncells_pml_domain+1; + } + else if ( iDim == 1 ) { + ny_p = ncells_pml_domain; + ny_d = ncells_pml_domain+1; + nx_p += ncells_pml_min[0] + ncells_pml_max[0]; + nx_d += ncells_pml_min[0] + ncells_pml_max[0]; + } + + //PML Coeffs Kappa,Sigma ... + //Primal + kappa_x_p.resize( nx_p ); + sigma_x_p.resize( nx_p ); + kappa_y_p.resize( ny_p ); + sigma_y_p.resize( ny_p ); + kappa_z_p.resize( 1 ); + sigma_z_p.resize( 1 ); + //Dual + kappa_x_d.resize( nx_d ); + sigma_x_d.resize( nx_d ); + kappa_y_d.resize( ny_d ); + sigma_y_d.resize( ny_d ); + kappa_z_d.resize( 1 ); + sigma_z_d.resize( 1 ); + + // While min and max of each dim have got the same ncells_pml, coefficient are the same + // for min and max PML + + //Coeffs for El,Dl,Hl,Bl + //Primal + c1_p_xfield.resize( ny_p ); // j-dependent + c2_p_xfield.resize( ny_p ); // j-dependent + c3_p_xfield.resize( 1 ); // k-dependent + c4_p_xfield.resize( 1 ); // k-dependent + c5_p_xfield.resize( nx_p ); // i-dependent + c6_p_xfield.resize( nx_p ); // i-dependent + //Dual + c1_d_xfield.resize( ny_d ); // j-dependent + c2_d_xfield.resize( ny_d ); // j-dependent + c3_d_xfield.resize( 1 ); // j-dependent + c4_d_xfield.resize( 1 ); // j-dependent + c5_d_xfield.resize( nx_d ); // i-dependent + c6_d_xfield.resize( nx_d ); // i-dependent + //Coefs for Er,Dr,Hr,Br + //Primal + c1_p_yfield.resize( 1 ); // k-dependent + c2_p_yfield.resize( 1 ); // k-dependent + c3_p_yfield.resize( nx_p ); // i-dependent + c4_p_yfield.resize( nx_p ); // i-dependent + c5_p_yfield.resize( ny_p ); // j-dependent + c6_p_yfield.resize( ny_p ); // j-dependent + //Dual + c1_d_yfield.resize( 1 ); // k-dependent + c2_d_yfield.resize( 1 ); // k-dependent + c3_d_yfield.resize( nx_d ); // i-dependent + c4_d_yfield.resize( nx_d ); // i-dependent + c5_d_yfield.resize( ny_d ); // k-dependent + c6_d_yfield.resize( ny_d ); // k-dependent + //Coefs for Et,Dt,Ht,Bt + //Primal + c1_p_zfield.resize( nx_p ); // j-dependent + c2_p_zfield.resize( nx_p ); // j-dependent + c3_p_zfield.resize( ny_p ); // i-dependent + c4_p_zfield.resize( ny_p ); // i-dependent + c5_p_zfield.resize( 1 ); // j-dependent + c6_p_zfield.resize( 1 ); // j-dependent + //Dual + c1_d_zfield.resize( nx_d ); // i-dependent + c2_d_zfield.resize( nx_d ); // i-dependent + c3_d_zfield.resize( ny_d ); // j-dependent + c4_d_zfield.resize( ny_d ); // j-dependent + c5_d_zfield.resize( 1 ); // k-dependent + c6_d_zfield.resize( 1 ); // k-dependent + + // Quote of the first primal grid-point where PML solver will be apply + // The first dual grid-point is at getDomainLocalMax( 0 ) - 0.5*dx and getDomainLocalMax( 1 ) - 0.5*dr + // xmax = patch->getDomainLocalMax( 0 ); + // ymax = patch->getDomainLocalMax( 1 ); + + if ( iDim == 0 ) { + // 3 cells (oversize) are vaccum so the PML media begin at r0 which is : + // Eventually the size of PML media is : + length_y_pml = 0 ; + length_x_pml = (ncells_pml_domain-startpml+0.5)*dx ; + // Primal grid + // X-direction + // Params for first cell of PML-patch (vacuum) i = 0,1,2 + for ( int i=0 ; i=3 + for ( int i=startpml; i=4 + for ( int i=startpml+1 ; i=3 + for ( int j=startpml ; j=4 + for ( int j=startpml+1 ; j( fields->emBoundCond[iDim*2+min_or_max] ); + Field2D* Ex_pml = NULL; + Field2D* Ey_pml = NULL; + Field2D* Ez_pml = NULL; + Field2D* Hx_pml = NULL; + Field2D* Hy_pml = NULL; + Field2D* Hz_pml = NULL; + Field2D* Dx_pml = NULL; + Field2D* Dy_pml = NULL; + Field2D* Dz_pml = NULL; + + Ex_pml = pml_fields->Ex_; + Ey_pml = pml_fields->Ey_; + Ez_pml = pml_fields->Ez_; + Hx_pml = pml_fields->Hx_; + Hy_pml = pml_fields->Hy_; + Hz_pml = pml_fields->Hz_; + Dx_pml = pml_fields->Dx_; + Dy_pml = pml_fields->Dy_; + Dz_pml = pml_fields->Dz_; + + if (min_or_max==0){ + isMin=true; + isMax=false; + } + else if (min_or_max==1){ + isMin=false; + isMax=true; + } + + if (iDim == 0) { + //Electric field Ex^(d,p,p) Remind that in PML, there no current + for( unsigned int k=0 ; k<1 ; k++ ) { + for( unsigned int i=solvermin ; i( fields->emBoundCond[iDim*2+min_or_max] ); + Field2D* Ex_pml = NULL; + Field2D* Ey_pml = NULL; + Field2D* Ez_pml = NULL; + Field2D* Hx_pml = NULL; + Field2D* Hy_pml = NULL; + Field2D* Hz_pml = NULL; + Field2D* Bx_pml = NULL; + Field2D* By_pml = NULL; + Field2D* Bz_pml = NULL; + + Ex_pml = pml_fields->Ex_; + Ey_pml = pml_fields->Ey_; + Ez_pml = pml_fields->Ez_; + Hx_pml = pml_fields->Hx_; + Hy_pml = pml_fields->Hy_; + Hz_pml = pml_fields->Hz_; + Bx_pml = pml_fields->Bx_; + By_pml = pml_fields->By_; + Bz_pml = pml_fields->Bz_; + + if (min_or_max==0){ + isMin=true; + isMax=false; + } + else if (min_or_max==1){ + isMin=false; + isMax=true; + } + + if (iDim==0){ + //Magnetic field Bx^(p,d,d) Remind that in PML, there no current + for( unsigned int k=0 ; k<1 ; k++ ) { + for( unsigned int i=solvermin ; i kappa_x_p; + std::vector sigma_x_p; + std::vector kappa_x_d; + std::vector sigma_x_d; + std::vector kappa_y_p; + std::vector sigma_y_p; + std::vector kappa_y_d; + std::vector sigma_y_d; + std::vector kappa_z_p; + std::vector sigma_z_p; + std::vector kappa_z_d; + std::vector sigma_z_d; + + std::vector c1_p_xfield; + std::vector c2_p_xfield; + std::vector c3_p_xfield; + std::vector c4_p_xfield; + std::vector c5_p_xfield; + std::vector c6_p_xfield; + std::vector c1_d_xfield; + std::vector c2_d_xfield; + std::vector c3_d_xfield; + std::vector c4_d_xfield; + std::vector c5_d_xfield; + std::vector c6_d_xfield; + std::vector c1_p_yfield; + std::vector c2_p_yfield; + std::vector c3_p_yfield; + std::vector c4_p_yfield; + std::vector c5_p_yfield; + std::vector c6_p_yfield; + std::vector c1_d_yfield; + std::vector c2_d_yfield; + std::vector c3_d_yfield; + std::vector c4_d_yfield; + std::vector c5_d_yfield; + std::vector c6_d_yfield; + std::vector c1_p_zfield; + std::vector c2_p_zfield; + std::vector c3_p_zfield; + std::vector c4_p_zfield; + std::vector c5_p_zfield; + std::vector c6_p_zfield; + std::vector c1_d_zfield; + std::vector c2_d_zfield; + std::vector c3_d_zfield; + std::vector c4_d_zfield; + std::vector c5_d_zfield; + std::vector c6_d_zfield; + + //double xmin; + //double ymin; + + double length_x_pml; + double length_y_pml; + double length_x_pml_xmin; + double length_x_pml_xmax; + + bool isMin ; + bool isMax ; + double Dx_pml_old ; + double Dy_pml_old ; + double Dz_pml_old ; + double Bx_pml_old ; + double By_pml_old ; + double Bz_pml_old ; + +};//END class + +#endif + diff --git a/src/ElectroMagnSolver/PML_Solver3D_Bouchard.cpp b/src/ElectroMagnSolver/PML_Solver3D_Bouchard.cpp new file mode 100644 index 000000000..fa1c5efb2 --- /dev/null +++ b/src/ElectroMagnSolver/PML_Solver3D_Bouchard.cpp @@ -0,0 +1,1169 @@ +#include "PML_Solver3D_Bouchard.h" +#include "ElectroMagn.h" +#include "ElectroMagnBC3D_PML.h" +#include "Field3D.h" +#include "Patch.h" + +PML_Solver3D_Bouchard::PML_Solver3D_Bouchard( Params ¶ms ) + : Solver3D( params ) +{ + //ERROR("Under development, not yet working"); + double dt = params.timestep; + dx = params.cell_length[0]; + dy = params.cell_length[1]; + dz = params.cell_length[2]; + double dx_ov_dt = dx/dt; + double dy_ov_dt = dy/dt; + double dz_ov_dt = dz/dt; + double dt_ov_dx = dt/dx; + double dt_ov_dy = dt/dy; + double dt_ov_dz = dt/dz; + //Not necessary to have dx=dy=dz, but dispersion law are modify + //In particular if dz >> dx,dy then solver become like the 2d solver + //if( (dx!=dy)||(dx!=dz)||(dy!=dz) ) { + // ERROR( "Bouchard solver requires the same cell-length in x, y and z directions" ); + //} + if( dx_ov_dt!=2 ) { + WARNING( "Bouchard solver requires dx/dt = 2 (Magical Timestep)" ); + } + + double delta = 0.1222*(1-pow(2.,2))/4. ; + double beta = -0.1727*(1-0.5*pow(2.,2)-4.*delta)/4. ; + double alpha = 1-4.*beta-3.*delta ; + + delta_x = delta ; + delta_y = delta ; + delta_z = delta ; + beta_xy = beta ; + beta_yx = beta ; + beta_xz = beta ; + beta_zx = beta ; + beta_yz = beta ; + beta_zy = beta ; + alpha_x = alpha ; + alpha_y = alpha ; + alpha_z = alpha ; + + // Ax = alpha_x*dt/dx; + // Ay = alpha_y*dt/dy; + // Az = alpha_z*dt/dz; + // Bxy = beta_xy*dt/dx; + // Byx = beta_yx*dt/dy; + // Bxz = beta_xz*dt/dx; + // Bzx = beta_zx*dt/dz; + // Byz = beta_yz*dt/dy; + // Bzy = beta_zy*dt/dz; + // Dx = delta_x*dt/dx; + // Dy = delta_y*dt/dy; + // Dz = delta_z*dt/dz; + + Ax = alpha_x/dx; + Ay = alpha_y/dy; + Az = alpha_z/dz; + Bxy = beta_xy/dx; + Byx = beta_yx/dy; + Bxz = beta_xz/dx; + Bzx = beta_zx/dz; + Byz = beta_yz/dy; + Bzy = beta_zy/dz; + Dx = delta_x/dx; + Dy = delta_y/dy; + Dz = delta_z/dz; + + //Define here the value of coefficient kappa_x_max, power_kappa_x, sigma_x_max, power_sigma_x + sigma_x_max = 20.; + kappa_x_max = 80.; + sigma_power_pml_x = 2.; + kappa_power_pml_x = 4.; + //Define here the value of coefficient kappa_y_max, power_kappa_y, sigma_y_max, power_sigma_y + sigma_y_max = 20.; + kappa_y_max = 80.; + sigma_power_pml_y = 2.; + kappa_power_pml_y = 4.; + //Define here the value of coefficient kappa_z_max, power_kappa_z, sigma_z_max, power_sigma_z + sigma_z_max = 20.; + kappa_z_max = 80.; + sigma_power_pml_z = 2.; + kappa_power_pml_z = 4.; +} + +PML_Solver3D_Bouchard::~PML_Solver3D_Bouchard() +{ +} + +void PML_Solver3D_Bouchard::operator()( ElectroMagn *fields ) +{ + ERROR( "This is not a solver for the main domain" ); +} + +void PML_Solver3D_Bouchard::setDomainSizeAndCoefficients( int iDim, int min_or_max, int ncells_pml_domain, int startpml, int* ncells_pml_min, int* ncells_pml_max, Patch* patch ) +{ + if ( iDim == 0 ) { + nx_p = ncells_pml_domain; + nx_d = ncells_pml_domain+1; + } + else if ( iDim == 1 ) { + ny_p = ncells_pml_domain; + ny_d = ncells_pml_domain+1; + nx_p += ncells_pml_min[0] + ncells_pml_max[0]; + nx_d += ncells_pml_min[0] + ncells_pml_max[0]; + } + + else if ( iDim == 2 ) { + nz_p = ncells_pml_domain; + nz_d = ncells_pml_domain+1; + nx_p += ncells_pml_min[0] + ncells_pml_max[0]; + nx_d += ncells_pml_min[0] + ncells_pml_max[0]; + ny_p += ncells_pml_min[1] + ncells_pml_max[1]; + ny_d += ncells_pml_min[1] + ncells_pml_max[1]; + } + + //PML Coeffs Kappa,Sigma ... + //Primal + kappa_x_p.resize( nx_p ); + sigma_x_p.resize( nx_p ); + kappa_y_p.resize( ny_p ); + sigma_y_p.resize( ny_p ); + kappa_z_p.resize( nz_p ); + sigma_z_p.resize( nz_p ); + //Dual + kappa_x_d.resize( nx_d ); + sigma_x_d.resize( nx_d ); + kappa_y_d.resize( ny_d ); + sigma_y_d.resize( ny_d ); + kappa_z_d.resize( nz_d ); + sigma_z_d.resize( nz_d ); + + // While min and max of each dim have got the same ncells_pml, coefficient are the same + // for min and max PML + + //Coeffs for Ex,Dx,Hx,Bx + //Primal + c1_p_xfield.resize( ny_p ); // j-dependent + c2_p_xfield.resize( ny_p ); // j-dependent + c3_p_xfield.resize( nz_p ); // k-dependent + c4_p_xfield.resize( nz_p ); // k-dependent + c5_p_xfield.resize( nx_p ); // i-dependent + c6_p_xfield.resize( nx_p ); // i-dependent + //Dual + c1_d_xfield.resize( ny_d ); // j-dependent + c2_d_xfield.resize( ny_d ); // j-dependent + c3_d_xfield.resize( nz_d ); // j-dependent + c4_d_xfield.resize( nz_d ); // j-dependent + c5_d_xfield.resize( nx_d ); // i-dependent + c6_d_xfield.resize( nx_d ); // i-dependent + //Coefs for Ey,Dy,Hy,By + //Primal + c1_p_yfield.resize( nz_p ); // k-dependent + c2_p_yfield.resize( nz_p ); // k-dependent + c3_p_yfield.resize( nx_p ); // i-dependent + c4_p_yfield.resize( nx_p ); // i-dependent + c5_p_yfield.resize( ny_p ); // j-dependent + c6_p_yfield.resize( ny_p ); // j-dependent + //Dual + c1_d_yfield.resize( nz_d ); // k-dependent + c2_d_yfield.resize( nz_d ); // k-dependent + c3_d_yfield.resize( nx_d ); // i-dependent + c4_d_yfield.resize( nx_d ); // i-dependent + c5_d_yfield.resize( ny_d ); // k-dependent + c6_d_yfield.resize( ny_d ); // k-dependent + //Coefs for Ez,Dz,Hz,Bz + //Primal + c1_p_zfield.resize( nx_p ); // j-dependent + c2_p_zfield.resize( nx_p ); // j-dependent + c3_p_zfield.resize( ny_p ); // i-dependent + c4_p_zfield.resize( ny_p ); // i-dependent + c5_p_zfield.resize( nz_p ); // j-dependent + c6_p_zfield.resize( nz_p ); // j-dependent + //Dual + c1_d_zfield.resize( nx_d ); // i-dependent + c2_d_zfield.resize( nx_d ); // i-dependent + c3_d_zfield.resize( ny_d ); // j-dependent + c4_d_zfield.resize( ny_d ); // j-dependent + c5_d_zfield.resize( nz_d ); // k-dependent + c6_d_zfield.resize( nz_d ); // k-dependent + + // Quote of the first primal grid-point where PML solver will be apply + // The first dual grid-point is at getDomainLocalMax( 0 ) - 0.5*dx and getDomainLocalMax( 1 ) - 0.5*dy + // xmax = patch->getDomainLocalMax( 0 ); + // ymax = patch->getDomainLocalMax( 1 ); + // zmax = patch->getDomainLocalMax( 2 ); + + if ( iDim == 0 ) { + // 3 cells (oversize) are vaccum so the PML media begin at y0 which is : + // Eventually the size of PML media is : + length_x_pml = (ncells_pml_domain-startpml+0.5)*dx ; + length_y_pml = 0. ; + length_z_pml = 0. ; + // Primal grid + // X-direction + // Params for first cell of PML-patch (vacuum) i = 0,1,2 + for ( int i=0 ; i=3 + // sigma_x_max = 80.; + // kappa_x_max = 80.; + // sigma_power_pml_x = 4.; + // kappa_power_pml_x = 4.; + for ( int i=startpml ; i=4 + // sigma_x_max = 80.; + // kappa_x_max = 80.; + // sigma_power_pml_x = 4.; + // kappa_power_pml_x = 4.; + for ( int i=startpml+1 ; i=3 + // sigma_y_max = 80.; + // kappa_y_max = 80.; + // sigma_power_pml_y = 4.; + // kappa_power_pml_y = 4.; + for ( int j=startpml ; j=4 + // sigma_y_max = 80.; + // kappa_y_max = 80.; + // sigma_power_pml_y = 4.; + // kappa_power_pml_y = 4.; + for ( int j=startpml+1 ; j=3 + // sigma_z_max = 80.; + // kappa_z_max = 80.; + // sigma_power_pml_z = 4.; + // kappa_power_pml_z = 4.; + for ( int k=startpml ; k=4 + // sigma_z_max = 80.; + // kappa_z_max = 80.; + // sigma_power_pml_z = 4.; + // kappa_power_pml_z = 4.; + for ( int k=startpml+1 ; k( fields->emBoundCond[iDim*2+min_or_max] ); + Field3D* Ex_pml = NULL; + Field3D* Ey_pml = NULL; + Field3D* Ez_pml = NULL; + Field3D* Hx_pml = NULL; + Field3D* Hy_pml = NULL; + Field3D* Hz_pml = NULL; + Field3D* Dx_pml = NULL; + Field3D* Dy_pml = NULL; + Field3D* Dz_pml = NULL; + + Ex_pml = pml_fields->Ex_; + Ey_pml = pml_fields->Ey_; + Ez_pml = pml_fields->Ez_; + Hx_pml = pml_fields->Hx_; + Hy_pml = pml_fields->Hy_; + Hz_pml = pml_fields->Hz_; + Dx_pml = pml_fields->Dx_; + Dy_pml = pml_fields->Dy_; + Dz_pml = pml_fields->Dz_; + + if (min_or_max==0){ + isMin=true; + isMax=false; + } + else if (min_or_max==1){ + isMin=false; + isMax=true; + } + + if (iDim == 0) { + //Electric field Ex^(d,p,p) Remind that in PML, there no current + for( unsigned int i=solvermin ; i( fields->emBoundCond[iDim*2+min_or_max] ); + Field3D* Ex_pml = NULL; + Field3D* Ey_pml = NULL; + Field3D* Ez_pml = NULL; + Field3D* Hx_pml = NULL; + Field3D* Hy_pml = NULL; + Field3D* Hz_pml = NULL; + Field3D* Bx_pml = NULL; + Field3D* By_pml = NULL; + Field3D* Bz_pml = NULL; + + Ex_pml = pml_fields->Ex_; + Ey_pml = pml_fields->Ey_; + Ez_pml = pml_fields->Ez_; + Hx_pml = pml_fields->Hx_; + Hy_pml = pml_fields->Hy_; + Hz_pml = pml_fields->Hz_; + Bx_pml = pml_fields->Bx_; + By_pml = pml_fields->By_; + Bz_pml = pml_fields->Bz_; + + if (min_or_max==0){ + isMin=true; + isMax=false; + } + else if (min_or_max==1){ + isMin=false; + isMax=true; + } + + if (iDim==0){ + //Magnetic field Bx^(p,d,d) Remind that in PML, there no current + for( unsigned int i=solvermin ; i kappa_x_p; + std::vector sigma_x_p; + std::vector kappa_x_d; + std::vector sigma_x_d; + std::vector kappa_y_p; + std::vector sigma_y_p; + std::vector kappa_y_d; + std::vector sigma_y_d; + std::vector kappa_z_p; + std::vector sigma_z_p; + std::vector kappa_z_d; + std::vector sigma_z_d; + + std::vector c1_p_xfield; + std::vector c2_p_xfield; + std::vector c3_p_xfield; + std::vector c4_p_xfield; + std::vector c5_p_xfield; + std::vector c6_p_xfield; + std::vector c1_d_xfield; + std::vector c2_d_xfield; + std::vector c3_d_xfield; + std::vector c4_d_xfield; + std::vector c5_d_xfield; + std::vector c6_d_xfield; + std::vector c1_p_yfield; + std::vector c2_p_yfield; + std::vector c3_p_yfield; + std::vector c4_p_yfield; + std::vector c5_p_yfield; + std::vector c6_p_yfield; + std::vector c1_d_yfield; + std::vector c2_d_yfield; + std::vector c3_d_yfield; + std::vector c4_d_yfield; + std::vector c5_d_yfield; + std::vector c6_d_yfield; + std::vector c1_p_zfield; + std::vector c2_p_zfield; + std::vector c3_p_zfield; + std::vector c4_p_zfield; + std::vector c5_p_zfield; + std::vector c6_p_zfield; + std::vector c1_d_zfield; + std::vector c2_d_zfield; + std::vector c3_d_zfield; + std::vector c4_d_zfield; + std::vector c5_d_zfield; + std::vector c6_d_zfield; + + //double xmax; + //double ymax; + //double zmax; + + double length_x_pml; + double length_y_pml; + double length_z_pml; + double length_x_pml_xmin; + double length_x_pml_xmax; + double length_y_pml_ymin; + double length_y_pml_ymax; + + bool isMin ; + bool isMax ; + double Dx_pml_old ; + double Dy_pml_old ; + double Dz_pml_old ; + double Bx_pml_old ; + double By_pml_old ; + double Bz_pml_old ; + +};//END class + +#endif diff --git a/src/ElectroMagnSolver/PML_Solver3D_Yee.cpp b/src/ElectroMagnSolver/PML_Solver3D_Yee.cpp new file mode 100755 index 000000000..bf45f855a --- /dev/null +++ b/src/ElectroMagnSolver/PML_Solver3D_Yee.cpp @@ -0,0 +1,1009 @@ +#include "PML_Solver3D_Yee.h" +#include "ElectroMagn.h" +#include "ElectroMagnBC3D_PML.h" +#include "Field3D.h" +#include "Patch.h" + +PML_Solver3D_Yee::PML_Solver3D_Yee( Params ¶ms ) + : Solver3D( params ) +{ + //Define here the value of coefficient kappa_x_max, power_kappa_x, sigma_x_max, power_sigma_x + sigma_x_max = 20.; + kappa_x_max = 80.; + sigma_power_pml_x = 2.; + kappa_power_pml_x = 4.; + //Define here the value of coefficient kappa_y_max, power_kappa_y, sigma_y_max, power_sigma_y + sigma_y_max = 20.; + kappa_y_max = 80.; + sigma_power_pml_y = 2.; + kappa_power_pml_y = 4.; + //Define here the value of coefficient kappa_z_max, power_kappa_z, sigma_z_max, power_sigma_z + sigma_z_max = 20.; + kappa_z_max = 80.; + sigma_power_pml_z = 2.; + kappa_power_pml_z = 4.; +} + +PML_Solver3D_Yee::~PML_Solver3D_Yee() +{ +} + +void PML_Solver3D_Yee::operator()( ElectroMagn *fields ) +{ + ERROR( "This is not a solver for the main domain" ); + + Field3D *Ex3D = static_cast( fields->Ex_ ); + Field3D *Ey3D = static_cast( fields->Ey_ ); + Field3D *Ez3D = static_cast( fields->Ez_ ); + Field3D *Bx3D = static_cast( fields->Bx_ ); + Field3D *By3D = static_cast( fields->By_ ); + Field3D *Bz3D = static_cast( fields->Bz_ ); + +} + +void PML_Solver3D_Yee::setDomainSizeAndCoefficients( int iDim, int min_or_max, int ncells_pml_domain, int startpml, int* ncells_pml_min, int* ncells_pml_max, Patch* patch ) +{ + if ( iDim == 0 ) { + nx_p = ncells_pml_domain; + nx_d = ncells_pml_domain+1; + } + else if ( iDim == 1 ) { + ny_p = ncells_pml_domain; + ny_d = ncells_pml_domain+1; + nx_p += ncells_pml_min[0] + ncells_pml_max[0]; + nx_d += ncells_pml_min[0] + ncells_pml_max[0]; + } + + else if ( iDim == 2 ) { + nz_p = ncells_pml_domain; + nz_d = ncells_pml_domain+1; + nx_p += ncells_pml_min[0] + ncells_pml_max[0]; + nx_d += ncells_pml_min[0] + ncells_pml_max[0]; + ny_p += ncells_pml_min[1] + ncells_pml_max[1]; + ny_d += ncells_pml_min[1] + ncells_pml_max[1]; + } + + //PML Coeffs Kappa,Sigma ... + //Primal + kappa_x_p.resize( nx_p ); + sigma_x_p.resize( nx_p ); + kappa_y_p.resize( ny_p ); + sigma_y_p.resize( ny_p ); + kappa_z_p.resize( nz_p ); + sigma_z_p.resize( nz_p ); + //Dual + kappa_x_d.resize( nx_d ); + sigma_x_d.resize( nx_d ); + kappa_y_d.resize( ny_d ); + sigma_y_d.resize( ny_d ); + kappa_z_d.resize( nz_d ); + sigma_z_d.resize( nz_d ); + + // While min and max of each dim have got the same ncells_pml, coefficient are the same + // for min and max PML + + //Coeffs for Ex,Dx,Hx,Bx + //Primal + c1_p_xfield.resize( ny_p ); // j-dependent + c2_p_xfield.resize( ny_p ); // j-dependent + c3_p_xfield.resize( nz_p ); // k-dependent + c4_p_xfield.resize( nz_p ); // k-dependent + c5_p_xfield.resize( nx_p ); // i-dependent + c6_p_xfield.resize( nx_p ); // i-dependent + //Dual + c1_d_xfield.resize( ny_d ); // j-dependent + c2_d_xfield.resize( ny_d ); // j-dependent + c3_d_xfield.resize( nz_d ); // j-dependent + c4_d_xfield.resize( nz_d ); // j-dependent + c5_d_xfield.resize( nx_d ); // i-dependent + c6_d_xfield.resize( nx_d ); // i-dependent + //Coefs for Ey,Dy,Hy,By + //Primal + c1_p_yfield.resize( nz_p ); // k-dependent + c2_p_yfield.resize( nz_p ); // k-dependent + c3_p_yfield.resize( nx_p ); // i-dependent + c4_p_yfield.resize( nx_p ); // i-dependent + c5_p_yfield.resize( ny_p ); // j-dependent + c6_p_yfield.resize( ny_p ); // j-dependent + //Dual + c1_d_yfield.resize( nz_d ); // k-dependent + c2_d_yfield.resize( nz_d ); // k-dependent + c3_d_yfield.resize( nx_d ); // i-dependent + c4_d_yfield.resize( nx_d ); // i-dependent + c5_d_yfield.resize( ny_d ); // k-dependent + c6_d_yfield.resize( ny_d ); // k-dependent + //Coefs for Ez,Dz,Hz,Bz + //Primal + c1_p_zfield.resize( nx_p ); // j-dependent + c2_p_zfield.resize( nx_p ); // j-dependent + c3_p_zfield.resize( ny_p ); // i-dependent + c4_p_zfield.resize( ny_p ); // i-dependent + c5_p_zfield.resize( nz_p ); // j-dependent + c6_p_zfield.resize( nz_p ); // j-dependent + //Dual + c1_d_zfield.resize( nx_d ); // i-dependent + c2_d_zfield.resize( nx_d ); // i-dependent + c3_d_zfield.resize( ny_d ); // j-dependent + c4_d_zfield.resize( ny_d ); // j-dependent + c5_d_zfield.resize( nz_d ); // k-dependent + c6_d_zfield.resize( nz_d ); // k-dependent + + // Quote of the first primal grid-point where PML solver will be apply + // The first dual grid-point is at getDomainLocalMax( 0 ) - 0.5*dx and getDomainLocalMax( 1 ) - 0.5*dy + // xmax = patch->getDomainLocalMax( 0 ); + // ymax = patch->getDomainLocalMax( 1 ); + // zmax = patch->getDomainLocalMax( 2 ); + + if ( iDim == 0 ) { + // 3 cells (oversize) are vaccum so the PML media begin at y0 which is : + // Eventually the size of PML media is : + length_x_pml = (ncells_pml_domain-startpml+0.5)*dx ; + length_y_pml = 0. ; + length_z_pml = 0. ; + // Primal grid + // X-direction + // Params for first cell of PML-patch (vacuum) i = 0,1,2 + for ( int i=0 ; i=3 + // sigma_x_max = 80.; + // kappa_x_max = 80.; + // sigma_power_pml_x = 4.; + // kappa_power_pml_x = 4.; + for ( int i=startpml ; i=4 + // sigma_x_max = 80.; + // kappa_x_max = 80.; + // sigma_power_pml_x = 4.; + // kappa_power_pml_x = 4.; + for ( int i=startpml+1 ; i=3 + // sigma_y_max = 80.; + // kappa_y_max = 80.; + // sigma_power_pml_y = 4.; + // kappa_power_pml_y = 4.; + for ( int j=startpml ; j=4 + // sigma_y_max = 80.; + // kappa_y_max = 80.; + // sigma_power_pml_y = 4.; + // kappa_power_pml_y = 4.; + for ( int j=startpml+1 ; j=3 + // sigma_z_max = 80.; + // kappa_z_max = 80.; + // sigma_power_pml_z = 4.; + // kappa_power_pml_z = 4.; + for ( int k=startpml ; k=4 + // sigma_z_max = 80.; + // kappa_z_max = 80.; + // sigma_power_pml_z = 4.; + // kappa_power_pml_z = 4.; + for ( int k=startpml+1 ; k( fields->emBoundCond[iDim*2+min_or_max] ); + Field3D* Ex_pml = NULL; + Field3D* Ey_pml = NULL; + Field3D* Ez_pml = NULL; + Field3D* Hx_pml = NULL; + Field3D* Hy_pml = NULL; + Field3D* Hz_pml = NULL; + Field3D* Dx_pml = NULL; + Field3D* Dy_pml = NULL; + Field3D* Dz_pml = NULL; + + Ex_pml = pml_fields->Ex_; + Ey_pml = pml_fields->Ey_; + Ez_pml = pml_fields->Ez_; + Hx_pml = pml_fields->Hx_; + Hy_pml = pml_fields->Hy_; + Hz_pml = pml_fields->Hz_; + Dx_pml = pml_fields->Dx_; + Dy_pml = pml_fields->Dy_; + Dz_pml = pml_fields->Dz_; + + if (min_or_max==0){ + isMin=true; + isMax=false; + } + else if (min_or_max==1){ + isMin=false; + isMax=true; + } + + if (iDim == 0) { + //Electric field Ex^(d,p,p) Remind that in PML, there no current + for( unsigned int i=solvermin ; i( fields->emBoundCond[iDim*2+min_or_max] ); + Field3D* Ex_pml = NULL; + Field3D* Ey_pml = NULL; + Field3D* Ez_pml = NULL; + Field3D* Hx_pml = NULL; + Field3D* Hy_pml = NULL; + Field3D* Hz_pml = NULL; + Field3D* Bx_pml = NULL; + Field3D* By_pml = NULL; + Field3D* Bz_pml = NULL; + + Ex_pml = pml_fields->Ex_; + Ey_pml = pml_fields->Ey_; + Ez_pml = pml_fields->Ez_; + Hx_pml = pml_fields->Hx_; + Hy_pml = pml_fields->Hy_; + Hz_pml = pml_fields->Hz_; + Bx_pml = pml_fields->Bx_; + By_pml = pml_fields->By_; + Bz_pml = pml_fields->Bz_; + + if (min_or_max==0){ + isMin=true; + isMax=false; + } + else if (min_or_max==1){ + isMin=false; + isMax=true; + } + + if (iDim==0){ + //Magnetic field Bx^(p,d,d) Remind that in PML, there no current + for( unsigned int i=solvermin ; i kappa_x_p; + std::vector sigma_x_p; + std::vector kappa_x_d; + std::vector sigma_x_d; + std::vector kappa_y_p; + std::vector sigma_y_p; + std::vector kappa_y_d; + std::vector sigma_y_d; + std::vector kappa_z_p; + std::vector sigma_z_p; + std::vector kappa_z_d; + std::vector sigma_z_d; + + std::vector c1_p_xfield; + std::vector c2_p_xfield; + std::vector c3_p_xfield; + std::vector c4_p_xfield; + std::vector c5_p_xfield; + std::vector c6_p_xfield; + std::vector c1_d_xfield; + std::vector c2_d_xfield; + std::vector c3_d_xfield; + std::vector c4_d_xfield; + std::vector c5_d_xfield; + std::vector c6_d_xfield; + std::vector c1_p_yfield; + std::vector c2_p_yfield; + std::vector c3_p_yfield; + std::vector c4_p_yfield; + std::vector c5_p_yfield; + std::vector c6_p_yfield; + std::vector c1_d_yfield; + std::vector c2_d_yfield; + std::vector c3_d_yfield; + std::vector c4_d_yfield; + std::vector c5_d_yfield; + std::vector c6_d_yfield; + std::vector c1_p_zfield; + std::vector c2_p_zfield; + std::vector c3_p_zfield; + std::vector c4_p_zfield; + std::vector c5_p_zfield; + std::vector c6_p_zfield; + std::vector c1_d_zfield; + std::vector c2_d_zfield; + std::vector c3_d_zfield; + std::vector c4_d_zfield; + std::vector c5_d_zfield; + std::vector c6_d_zfield; + + //double xmax; + //double ymax; + //double zmax; + + double length_x_pml; + double length_y_pml; + double length_z_pml; + double length_x_pml_xmin; + double length_x_pml_xmax; + double length_y_pml_ymin; + double length_y_pml_ymax; + + bool isMin ; + bool isMax ; + double Dx_pml_old ; + double Dy_pml_old ; + double Dz_pml_old ; + double Bx_pml_old ; + double By_pml_old ; + double Bz_pml_old ; + +};//END class + +#endif + diff --git a/src/ElectroMagnSolver/PML_SolverAM.cpp b/src/ElectroMagnSolver/PML_SolverAM.cpp new file mode 100755 index 000000000..1357373a9 --- /dev/null +++ b/src/ElectroMagnSolver/PML_SolverAM.cpp @@ -0,0 +1,1049 @@ +#include "PML_SolverAM.h" +#include "ElectroMagnAM.h" +#include "ElectroMagnBCAM_PML.h" +#include "cField2D.h" +#include +#include "dcomplex.h" +#include "Patch.h" + +PML_SolverAM::PML_SolverAM( Params ¶ms ) + : SolverAM( params ) +{ +} + +PML_SolverAM::~PML_SolverAM() +{ +} + +void PML_SolverAM::operator()( ElectroMagn *fields ) +{ + ERROR( "This is not a solver for the main domain" ); + + for( unsigned int imode=0 ; imode( fields ) )->Dl_[imode]; + //cField2D *Dr = ( static_cast( fields ) )->Dr_[imode]; + //cField2D *Dt = ( static_cast( fields ) )->Dt_[imode]; + cField2D *El = ( static_cast( fields ) )->El_[imode]; + cField2D *Er = ( static_cast( fields ) )->Er_[imode]; + cField2D *Et = ( static_cast( fields ) )->Et_[imode]; + cField2D *Bl = ( static_cast( fields ) )->Bl_[imode]; + cField2D *Br = ( static_cast( fields ) )->Br_[imode]; + cField2D *Bt = ( static_cast( fields ) )->Bt_[imode]; + //cField2D *Hl = ( static_cast( fields ) )->Hl_[imode]; + //cField2D *Hr = ( static_cast( fields ) )->Hr_[imode]; + //cField2D *Ht = ( static_cast( fields ) )->Ht_[imode]; + cField2D *Jl = ( static_cast( fields ) )->Jl_[imode]; + cField2D *Jr = ( static_cast( fields ) )->Jr_[imode]; + cField2D *Jt = ( static_cast( fields ) )->Jt_[imode]; + int j_glob = ( static_cast( fields ) )->j_glob_; + bool isYmin = ( static_cast( fields ) )->isYmin; + } +} + +void PML_SolverAM::setDomainSizeAndCoefficients( int iDim, int min_or_max, int ncells_pml_domain, int startpml, int* ncells_pml_min, int* ncells_pml_max, Patch* patch ) +{ + if ( iDim == 0 ) { + // Global radial index where begin the PML domain + // because j_glob is not define for this region + if (min_or_max==0) { + j_glob_pml = patch->getCellStartingGlobalIndex( 1 ); + } + else if (min_or_max==1) { + j_glob_pml = patch->getCellStartingGlobalIndex( 1 ); + } + nl_p = ncells_pml_domain; + nl_d = ncells_pml_domain+1; + } + else if ( iDim == 1 ) { + // Global radial index where begin the PML domain + // because j_glob is not define for this region + if (min_or_max==0) { + j_glob_pml = patch->getCellStartingGlobalIndex( 1 )-ncells_pml_domain+oversize[iDim]+1; + } + else if (min_or_max==1) { + j_glob_pml = patch->getCellStartingGlobalIndex( 1 )+nr_p-oversize[iDim]-1; + } + // Redifine length of pml region + nr_p = ncells_pml_domain; + nr_d = ncells_pml_domain+1; + nl_p += ncells_pml_min[0] + ncells_pml_max[0]; + nl_d += ncells_pml_min[0] + ncells_pml_max[0]; + } + + //PML Coeffs Kappa,Sigma ... + //Primal + kappa_r_p.resize( nr_p ); + sigma_r_p.resize( nr_p ); + kappa_l_p.resize( nl_p ); + sigma_l_p.resize( nl_p ); + integrate_kappa_r_p.resize( nr_p ) ; + integrate_sigma_r_p.resize( nr_p ) ; + //Dual + kappa_r_d.resize( nr_d ); + sigma_r_d.resize( nr_d ); + kappa_l_d.resize( nl_d ); + sigma_l_d.resize( nl_d ); + integrate_kappa_r_d.resize( nr_d ) ; + integrate_sigma_r_d.resize( nr_d ) ; + + //Coeffs for El,Dl,Hl,Bl + //Primal + c1_p_lfield.resize( nr_p, 1. ); // j-dependent + c2_p_lfield.resize( nr_p, dt ); // j-dependent + c3_p_lfield.resize( nr_p, 0. ); // j-dependent + c4_p_lfield.resize( nr_p, 1. ); // j-dependent + c5_p_lfield.resize( nl_p, 1. ); // i-dependent + c6_p_lfield.resize( nl_p, 0. ); // i-dependent + //Dual + c1_d_lfield.resize( nr_d, 1. ); // j-dependent + c2_d_lfield.resize( nr_d, dt ); // j-dependent + c3_d_lfield.resize( nr_d, 0 ); // j-dependent + c4_d_lfield.resize( nr_d, 1. ); // j-dependent + c5_d_lfield.resize( nl_d, 1. ); // i-dependent + c6_d_lfield.resize( nl_d, 0. ); // i-dependent + //Coefs for Er,Dr,Hr,Br + //Primal + c1_p_rfield.resize( nl_p, 1. ); // i-dependent + c2_p_rfield.resize( nl_p, dt ); // i-dependent + c3_p_rfield.resize( nr_p, 0. ); // j-dependent + c4_p_rfield.resize( nr_p, 1. ); // j-dependent + c5_p_rfield.resize( nr_p, 1. ); // j-dependent + c6_p_rfield.resize( nr_p, 0. ); // j-dependent + //Dual + c1_d_rfield.resize( nl_d, 1. ); // i-dependent + c2_d_rfield.resize( nl_d, dt ); // i-dependent + c3_d_rfield.resize( nr_d, 0. ); // j-dependent + c4_d_rfield.resize( nr_d, 1. ); // j-dependent + c5_d_rfield.resize( nr_d, 1. ); // j-dependent + c6_d_rfield.resize( nr_d, 0. ); // j-dependent + //Coefs for Et,Dt,Ht,Bt + //Primal + c1_p_tfield.resize( nr_p, 1. ); // j-dependent + c2_p_tfield.resize( nr_p, dt ); // j-dependent + c3_p_tfield.resize( nl_p, 0. ); // i-dependent + c4_p_tfield.resize( nl_p, 1. ); // i-dependent + c5_p_tfield.resize( nr_p, 1. ); // j-dependent + c6_p_tfield.resize( nr_p, 0. ); // j-dependent + //Dual + c1_d_tfield.resize( nr_d, 1. ); // j-dependent + c2_d_tfield.resize( nr_d, dt ); // j-dependent + c3_d_tfield.resize( nl_d, 0. ); // i-dependent + c4_d_tfield.resize( nl_d, 1. ); // i-dependent + c5_d_tfield.resize( nr_d, 1. ); // j-dependent + c6_d_tfield.resize( nr_d, 0. ); // j-dependent + + // if ( iDim == 0 ) { + // rmax = 0. ; + // r0 = 0. ; + // } + // else if( iDim == 1 ) { + // rmax = patch->getDomainLocalMax( 1 ) ; + // r0 = rmax + 3.*dr ; // Have to be oversize+1 + // } + + rmax = patch->getDomainLocalMax( 1 ) ; + r0 = rmax + (oversize[1] + 1 )*dr ; + + if ( iDim == 0 ) { + // 2 cells are vaccum so the PML media begin at r0 which is : + // Eventually the size of PML media is : + length_l_pml = (ncells_pml_domain-startpml+0.5)*dl ; + length_r_pml = 0 ; + // Primal grid + // Longitudinal + // Params for first cell of PML-patch (vacuum) i = 0,1,2 + for ( int i=0 ; i=3 + sigma_l_max = 80.; + kappa_l_max = 80.; + power_pml_l = 4.; + // sigma_l_max = 0.; + // kappa_l_max = 1.; + // power_pml_l = 0.; + for ( int i=startpml ; i=4 + sigma_l_max = 80.; + kappa_l_max = 80.; + power_pml_r = 4.; + // sigma_l_max = 0.; + // kappa_l_max = 1.; + // power_pml_l = 0.; + for ( int i=startpml+1 ; i=3 + sigma_r_max = 80.; + kappa_r_max = 80.; + power_pml_r = 4.; + // sigma_r_max = 0. ; + // kappa_r_max = 1. ; + // power_pml_r = 0. ; + for ( int j=startpml ; j=4 + sigma_r_max = 80.; + kappa_r_max = 80.; + power_pml_r = 4.; + // sigma_r_max = 0.; + // kappa_r_max = 1.; + // power_pml_r = 0.; + for ( int j=startpml+1 ; j( fields->emBoundCond[iDim*2+min_or_max] ); + cField2D* El_pml = NULL; + cField2D* Er_pml = NULL; + cField2D* Et_pml = NULL; + cField2D* Hl_pml = NULL; + cField2D* Hr_pml = NULL; + cField2D* Ht_pml = NULL; + cField2D* Dl_pml = NULL; + cField2D* Dr_pml = NULL; + cField2D* Dt_pml = NULL; + + bool isYmin = ( static_cast( fields ) )->isYmin; + + if (min_or_max==0){ + isMin=true; + isMax=false; + } + else if (min_or_max==1){ + isMin=false; + isMax=true; + } + + if (iDim==0){ + for( unsigned int imode=0 ; imodeEl_[imode]; + Er_pml = pml_fields->Er_[imode]; + Et_pml = pml_fields->Et_[imode]; + Hl_pml = pml_fields->Hl_[imode]; + Hr_pml = pml_fields->Hr_[imode]; + Ht_pml = pml_fields->Ht_[imode]; + Dl_pml = pml_fields->Dl_[imode]; + Dr_pml = pml_fields->Dr_[imode]; + Dt_pml = pml_fields->Dt_[imode]; + //Electric field El^(d,p) Remind that in PML, there no current + for( unsigned int i=solvermin ; i(1,0)*( *Dl_pml )( i, j ) ; + ( *Dl_pml )( i, j ) = + c1_p_lfield[j] * ( *Dl_pml )( i, j ) + + c2_p_lfield[j] / ( ( j_glob_pml+j )*dr ) * ( ( j+j_glob_pml+0.5 ) * ( *Ht_pml )( i, j+1 ) - ( j+j_glob_pml-0.5 )*( *Ht_pml )( i, j ) ) + + c2_p_lfield[j] / ( ( j_glob_pml+j )*dr ) * Icpx*( double )imode*( *Hr_pml )( i, j ) ; + ( *El_pml )( i, j ) = + c3_p_lfield[j] * ( *El_pml )( i, j ) + + c4_p_lfield[j] * (c5_d_lfield[i]*( *Dl_pml )( i, j ) - c6_d_lfield[i]*Dl_pml_old ) ; + } + } + //Electric field Er^(p,d) Remind that in PML, there no current + for( unsigned int i=solvermin ; i(1,0)*( *Dr_pml )( i, j ) ; + ( *Dr_pml )( i, j ) = + c1_p_rfield[i] * ( *Dr_pml )( i, j ) + - c2_p_rfield[i]/dl * ( ( *Ht_pml )( i+1, j ) - ( *Ht_pml )( i, j ) ) + - c2_p_rfield[i]*Icpx*( double )imode/( ( j+j_glob_pml-0.5 )*dr)*( *Hl_pml )( i, j ) ; + ( *Er_pml )( i, j ) = + c3_d_rfield[j] * ( *Er_pml )( i, j ) + + c4_d_rfield[j] * (c5_d_rfield[j]*( *Dr_pml )( i, j ) - c6_d_rfield[j]*Dr_pml_old ); + } + } + //Electric field Et^(p,p) Remind that in PML, there no current + for( unsigned int i=solvermin ; i(1,0)*( *Dt_pml )( i, j ) ; + ( *Dt_pml )( i, j ) = + c1_p_tfield[j] * ( *Dt_pml )( i, j ) + - c2_p_tfield[j]/dr * ( ( *Hl_pml )( i, j+1 ) - ( *Hl_pml )( i, j ) ) + + c2_p_tfield[j]/dl * ( ( *Hr_pml )( i+1, j ) - ( *Hr_pml )( i, j ) ) ; + ( *Et_pml )( i, j ) = + c3_p_tfield[i] * ( *Et_pml )( i, j ) + + c4_p_tfield[i] * (c5_p_tfield[j]*( *Dt_pml )( i, j ) - c6_p_tfield[j]*Dt_pml_old ); + } + } + if( isYmin ) { + // Conditions on axis + unsigned int j=2; + if( imode==0 ) { + for( unsigned int i=0 ; i 1 + for( unsigned int i=0 ; i 1 + // for( unsigned int i=3*isMax ; iEl_[imode]; + Er_pml = pml_fields->Er_[imode]; + Et_pml = pml_fields->Et_[imode]; + Hl_pml = pml_fields->Hl_[imode]; + Hr_pml = pml_fields->Hr_[imode]; + Ht_pml = pml_fields->Ht_[imode]; + Dl_pml = pml_fields->Dl_[imode]; + Dr_pml = pml_fields->Dr_[imode]; + Dt_pml = pml_fields->Dt_[imode]; + //Electric field El^(d,p) Remind that in PML, there no current + for( unsigned int i=0 ; i(1,0)*( *Dl_pml )( i, j ) ; + ( *Dl_pml )( i, j ) = + c1_p_lfield[j] * ( *Dl_pml )( i, j ) + + c2_p_lfield[j] / ( ( j_glob_pml+j )*dr ) * ( ( j+j_glob_pml+0.5 ) * ( *Ht_pml )( i, j+1 ) - ( j+j_glob_pml-0.5 )*( *Ht_pml )( i, j ) ) + + c2_p_lfield[j] / ( ( j_glob_pml+j )*dr ) * Icpx*( double )imode*( *Hr_pml )( i, j ) ; + ( *El_pml )( i, j ) = + c3_p_lfield[j] * ( *El_pml )( i, j ) + + c4_p_lfield[j] * (c5_d_lfield[i]*( *Dl_pml )( i, j ) - c6_d_lfield[i]*Dl_pml_old ) ; + } + } + //Electric field Er^(p,d) Remind that in PML, there no current + for( unsigned int i=0 ; i(1,0)*( *Dr_pml )( i, j ) ; + ( *Dr_pml )( i, j ) = + c1_p_rfield[i] * ( *Dr_pml )( i, j ) + - c2_p_rfield[i]/dl * ( ( *Ht_pml )( i+1, j ) - ( *Ht_pml )( i, j ) ) + - c2_p_rfield[i]*Icpx*( double )imode/( ( j+j_glob_pml-0.5 )*dr)*( *Hl_pml )( i, j ) ; + ( *Er_pml )( i, j ) = + c3_d_rfield[j] * ( *Er_pml )( i, j ) + + c4_d_rfield[j] * (c5_d_rfield[j]*( *Dr_pml )( i, j ) - c6_d_rfield[j]*Dr_pml_old ); + } + } + //Electric field Et^(p,p) Remind that in PML, there no current + for( unsigned int i=0 ; i(1,0)*( *Dt_pml )( i, j ) ; + ( *Dt_pml )( i, j ) = + c1_p_tfield[j] * ( *Dt_pml )( i, j ) + - c2_p_tfield[j]/dr * ( ( *Hl_pml )( i, j+1 ) - ( *Hl_pml )( i, j ) ) + + c2_p_tfield[j]/dl * ( ( *Hr_pml )( i+1, j ) - ( *Hr_pml )( i, j ) ) ; + ( *Et_pml )( i, j ) = + c3_p_tfield[i] * ( *Et_pml )( i, j ) + + c4_p_tfield[i] * (c5_p_tfield[j]*( *Dt_pml )( i, j ) - c6_p_tfield[j]*Dt_pml_old ); + } + } + } + } +} + +void PML_SolverAM::compute_H_from_B( ElectroMagn *fields, int iDim, int min_or_max, int solvermin, int solvermax ) +{ + ElectroMagnBCAM_PML* pml_fields = static_cast( fields->emBoundCond[iDim*2+min_or_max] ); + cField2D* El_pml = NULL; + cField2D* Er_pml = NULL; + cField2D* Et_pml = NULL; + cField2D* Hl_pml = NULL; + cField2D* Hr_pml = NULL; + cField2D* Ht_pml = NULL; + cField2D* Bl_pml = NULL; + cField2D* Br_pml = NULL; + cField2D* Bt_pml = NULL; + + bool isYmin = ( static_cast( fields ) )->isYmin; + + if (min_or_max==0){ + isMin=true; + isMax=false; + } + else if (min_or_max==1){ + isMin=false; + isMax=true; + } + + if (iDim==0){ + for( unsigned int imode=0 ; imodeEl_[imode]; + Er_pml = pml_fields->Er_[imode]; + Et_pml = pml_fields->Et_[imode]; + Hl_pml = pml_fields->Hl_[imode]; + Hr_pml = pml_fields->Hr_[imode]; + Ht_pml = pml_fields->Ht_[imode]; + Bl_pml = pml_fields->Bl_[imode]; + Br_pml = pml_fields->Br_[imode]; + Bt_pml = pml_fields->Bt_[imode]; + //Magnetic field Bl^(p,d) + for( unsigned int i=solvermin ; i(1,0)*( *Bl_pml )( i, j ) ; + ( *Bl_pml )( i, j ) = + c1_d_lfield[j] * ( *Bl_pml )( i, j ) + - c2_d_lfield[j] / ( ( j_glob_pml+j-0.5 )*dr ) * ( ( double )( j+j_glob_pml )*( *Et_pml )( i, j ) - ( double )( j+j_glob_pml-1. )*( *Et_pml )( i, j-1 ) ) + - c2_d_lfield[j] / ( ( j_glob_pml+j-0.5 )*dr ) * ( Icpx*( double )imode*( *Er_pml )( i, j ) ) ; + ( *Hl_pml )( i, j ) = + c3_d_lfield[j] * ( *Hl_pml )( i, j ) + + c4_d_lfield[j] * (c5_p_lfield[i]*( *Bl_pml )( i, j ) - c6_p_lfield[i]*Bl_pml_old ) ; + } + } + //Magnetic field Br^(d,p) + for( unsigned int i=solvermin ; i(1,0)*( *Br_pml )( i, j ) ; + ( *Br_pml )( i, j ) = + c1_d_rfield[i] * ( *Br_pml )( i, j ) + + c2_d_rfield[i]/dl * ( ( *Et_pml )( i, j ) - ( *Et_pml )( i-1, j ) ) + + c2_d_rfield[i]*Icpx*( double )imode/( ( double )( j_glob_pml+j )*dr )*( *El_pml )( i, j ) ; + ( *Hr_pml )( i, j ) = + c3_p_rfield[j] * ( *Hr_pml )( i, j ) + + c4_p_rfield[j] * (c5_p_rfield[j]*( *Br_pml )( i, j ) - c6_p_rfield[j]*Br_pml_old ); + } + } + //Magnetic field Bt^(d,d) + for( unsigned int i=solvermin ; i(1,0)*( *Bt_pml )( i, j ) ; + ( *Bt_pml )( i, j ) = + c1_d_tfield[j] * ( *Bt_pml )( i, j ) + + c2_d_tfield[j]/dr * ( ( *El_pml )( i, j ) - ( *El_pml )( i, j-1 ) ) + - c2_d_tfield[j]/dl * ( ( *Er_pml )( i, j ) - ( *Er_pml )( i-1, j ) ) ; + ( *Ht_pml )( i, j ) = + c3_d_tfield[i] * ( *Ht_pml )( i, j ) + + c4_d_tfield[i] * (c5_d_tfield[j]*( *Bt_pml )( i, j ) - c6_d_tfield[j]*Bt_pml_old ); + } + } + if( isYmin ) { + unsigned int j=2; + if( imode==0 ) { + for( unsigned int i=0 ; i 1 + for( unsigned int i=0 ; i 1 + // for( unsigned int i=3*isMax ; iEl_[imode]; + Er_pml = pml_fields->Er_[imode]; + Et_pml = pml_fields->Et_[imode]; + Hl_pml = pml_fields->Hl_[imode]; + Hr_pml = pml_fields->Hr_[imode]; + Ht_pml = pml_fields->Ht_[imode]; + Bl_pml = pml_fields->Bl_[imode]; + Br_pml = pml_fields->Br_[imode]; + Bt_pml = pml_fields->Bt_[imode]; + //Magnetic field Bl^(p,d) + for( unsigned int i=0 ; i(1,0)*( *Bl_pml )( i, j ) ; + ( *Bl_pml )( i, j ) = + c1_d_lfield[j] * ( *Bl_pml )( i, j ) + - c2_d_lfield[j] / ( ( j_glob_pml+j-0.5 )*dr ) * ( ( double )( j+j_glob_pml )*( *Et_pml )( i, j ) - ( double )( j+j_glob_pml-1. )*( *Et_pml )( i, j-1 ) ) + - c2_d_lfield[j] / ( ( j_glob_pml+j-0.5 )*dr ) * ( Icpx*( double )imode*( *Er_pml )( i, j ) ) ; + ( *Hl_pml )( i, j ) = + c3_d_lfield[j] * ( *Hl_pml )( i, j ) + + c4_d_lfield[j] * (c5_p_lfield[i]*( *Bl_pml )( i, j ) - c6_p_lfield[i]*Bl_pml_old ) ; + } + } + //Magnetic field Br^(d,p) + for( unsigned int i=1 ; i(1,0)*( *Br_pml )( i, j ) ; + ( *Br_pml )( i, j ) = + c1_d_rfield[i] * ( *Br_pml )( i, j ) + + c2_d_rfield[i]/dl * ( ( *Et_pml )( i, j ) - ( *Et_pml )( i-1, j ) ) + + c2_d_rfield[i]*Icpx*( double )imode/( ( double )( j_glob_pml+j )*dr )*( *El_pml )( i, j ) ; + ( *Hr_pml )( i, j ) = + c3_p_rfield[j] * ( *Hr_pml )( i, j ) + + c4_p_rfield[j] * (c5_p_rfield[j]*( *Br_pml )( i, j ) - c6_p_rfield[j]*Br_pml_old ); + } + } + //Magnetic field Bt^(d,d) + for( unsigned int i=1 ; i(1,0)*( *Bt_pml )( i, j ) ; + ( *Bt_pml )( i, j ) = + c1_d_tfield[j] * ( *Bt_pml )( i, j ) + + c2_d_tfield[j]/dr * ( ( *El_pml )( i, j ) - ( *El_pml )( i, j-1 ) ) + - c2_d_tfield[j]/dl * ( ( *Er_pml )( i, j ) - ( *Er_pml )( i-1, j ) ) ; + ( *Ht_pml )( i, j ) = + c3_d_tfield[i] * ( *Ht_pml )( i, j ) + + c4_d_tfield[i] * (c5_d_tfield[j]*( *Bt_pml )( i, j ) - c6_d_tfield[j]*Bt_pml_old ); + } + } + } + } +} diff --git a/src/ElectroMagnSolver/PML_SolverAM.h b/src/ElectroMagnSolver/PML_SolverAM.h new file mode 100755 index 000000000..ffb1ee450 --- /dev/null +++ b/src/ElectroMagnSolver/PML_SolverAM.h @@ -0,0 +1,105 @@ +#ifndef PML_SOLVERAM_H +#define PML_SOLVERAM_H + +#include "SolverAM.h" +class ElectroMagn; + +// -------------------------------------------------------------------------------------------------------------------- +//! Class Pusher +// -------------------------------------------------------------------------------------------------------------------- +class PML_SolverAM : public SolverAM +{ + +public: + PML_SolverAM( Params ¶ms ); + virtual ~PML_SolverAM(); + + //! Overloading of () operator + virtual void operator()( ElectroMagn *fields ); + + void setDomainSizeAndCoefficients( int iDim, int min_or_max, int ncells_pml_domain, int startpml, int* ncells_pml_lmin, int* ncells_pml_lmax, Patch* patch ); + + void compute_E_from_D( ElectroMagn *fields, int iDim, int min_or_max, int solvermin, int solvermax ); + void compute_H_from_B( ElectroMagn *fields, int iDim, int min_or_max, int solvermin, int solvermax ); + +protected: + double sigma_r_max; + double kappa_r_max; + double power_pml_r; + double sigma_l_max; + double kappa_l_max; + double power_pml_l; + + std::vector kappa_r_p; + std::vector sigma_r_p; + std::vector kappa_r_d; + std::vector sigma_r_d; + std::vector integrate_kappa_r_p; + std::vector integrate_sigma_r_p; + std::vector integrate_kappa_r_d; + std::vector integrate_sigma_r_d; + std::vector kappa_l_p; + std::vector sigma_l_p; + std::vector kappa_l_d; + std::vector sigma_l_d; + + std::vector c1_p_lfield; + std::vector c2_p_lfield; + std::vector c3_p_lfield; + std::vector c4_p_lfield; + std::vector c5_p_lfield; + std::vector c6_p_lfield; + std::vector c1_d_lfield; + std::vector c2_d_lfield; + std::vector c3_d_lfield; + std::vector c4_d_lfield; + std::vector c5_d_lfield; + std::vector c6_d_lfield; + std::vector c1_p_rfield; + std::vector c2_p_rfield; + std::vector c3_p_rfield; + std::vector c4_p_rfield; + std::vector c5_p_rfield; + std::vector c6_p_rfield; + std::vector c1_d_rfield; + std::vector c2_d_rfield; + std::vector c3_d_rfield; + std::vector c4_d_rfield; + std::vector c5_d_rfield; + std::vector c6_d_rfield; + std::vector c1_p_tfield; + std::vector c2_p_tfield; + std::vector c3_p_tfield; + std::vector c4_p_tfield; + std::vector c5_p_tfield; + std::vector c6_p_tfield; + std::vector c1_d_tfield; + std::vector c2_d_tfield; + std::vector c3_d_tfield; + std::vector c4_d_tfield; + std::vector c5_d_tfield; + std::vector c6_d_tfield; + + // double lmax; + // double l0; + int j_glob_pml; + double rmax; + double r0; + double length_r_pml; + double length_l_pml; + double length_l_pml_lmin; + double length_l_pml_lmax; + + bool isMin; + bool isMax; + + std::complex Dl_pml_old; + std::complex Dr_pml_old; + std::complex Dt_pml_old; + std::complex Bl_pml_old; + std::complex Br_pml_old; + std::complex Bt_pml_old; + +};//END class + +#endif diff --git a/src/ElectroMagnSolver/Solver.h b/src/ElectroMagnSolver/Solver.h index 7f0b772e2..8d52ac3e7 100755 --- a/src/ElectroMagnSolver/Solver.h +++ b/src/ElectroMagnSolver/Solver.h @@ -22,7 +22,11 @@ class Solver virtual void densities_correction( ElectroMagn *EMfields ) {}; //! Overloading of () operator virtual void operator()( ElectroMagn *fields ) = 0; - + + virtual void setDomainSizeAndCoefficients( int iDim, int min_or_max, int ncells_pml, int startpml, int* ncells_pml_min, int* ncells_pml_max, Patch* patch ) {ERROR("Not using PML");}; + virtual void compute_E_from_D( ElectroMagn *fields, int iDim, int min_or_max, int solver_min, int solver_max ) {ERROR("Not using PML");}; + virtual void compute_H_from_B( ElectroMagn *fields, int iDim, int min_or_max, int solver_min, int solver_max ) {ERROR("Not using PML");}; + protected: };//END class diff --git a/src/ElectroMagnSolver/Solver2D.h b/src/ElectroMagnSolver/Solver2D.h index 8b25616c7..020b53851 100755 --- a/src/ElectroMagnSolver/Solver2D.h +++ b/src/ElectroMagnSolver/Solver2D.h @@ -27,6 +27,8 @@ class Solver2D : public Solver ny_d = n_space[1] +2+2*oversize[1]; dt = params.timestep; + dx = params.cell_length[0]; + dy = params.cell_length[1]; dt_ov_dx = params.timestep / params.cell_length[0]; dt_ov_dy = params.timestep / params.cell_length[1]; }; @@ -41,6 +43,8 @@ class Solver2D : public Solver unsigned int ny_p; unsigned int ny_d; double dt; + double dx; + double dy; double dt_ov_dy; double dt_ov_dx; diff --git a/src/ElectroMagnSolver/Solver3D.h b/src/ElectroMagnSolver/Solver3D.h index 5389ccd69..33989833a 100755 --- a/src/ElectroMagnSolver/Solver3D.h +++ b/src/ElectroMagnSolver/Solver3D.h @@ -29,6 +29,9 @@ class Solver3D : public Solver nz_d = n_space[2] +2+2*oversize[2]-(params.is_pxr); dt = params.timestep; + dx = params.cell_length[0]; + dy = params.cell_length[1]; + dz = params.cell_length[2]; dt_ov_dx = params.timestep / params.cell_length[0]; dt_ov_dy = params.timestep / params.cell_length[1]; dt_ov_dz = params.timestep / params.cell_length[2]; @@ -47,6 +50,9 @@ class Solver3D : public Solver unsigned int nz_p; unsigned int nz_d; double dt; + double dx; + double dy; + double dz; double dt_ov_dx; double dt_ov_dy; double dt_ov_dz; diff --git a/src/ElectroMagnSolver/SolverAM.h b/src/ElectroMagnSolver/SolverAM.h index ec0762599..92c42bc4b 100755 --- a/src/ElectroMagnSolver/SolverAM.h +++ b/src/ElectroMagnSolver/SolverAM.h @@ -13,7 +13,7 @@ class SolverAM : public Solver //! Creator for Solver SolverAM( Params ¶ms ) : Solver( params ) { - std::vector oversize(params.oversize); + oversize = params.oversize; if (params.multiple_decomposition) oversize = params.region_oversize; nl_p = params.n_space[0]+1+2*oversize[0]; @@ -33,6 +33,7 @@ class SolverAM : public Solver Nmode= params.nmodes; dt = params.timestep; + dl = params.cell_length[0]; dr = params.cell_length[1]; dt_ov_dl = params.timestep / params.cell_length[0]; dt_ov_dr = params.timestep / params.cell_length[1]; @@ -50,9 +51,11 @@ class SolverAM : public Solver unsigned int nr_d; unsigned int Nmode; double dt; + double dl; double dr; double dt_ov_dl; double dt_ov_dr; + std::vector oversize; };//END class diff --git a/src/ElectroMagnSolver/SolverFactory.h b/src/ElectroMagnSolver/SolverFactory.h index 3cc3ee3d1..ea1ef5e3f 100755 --- a/src/ElectroMagnSolver/SolverFactory.h +++ b/src/ElectroMagnSolver/SolverFactory.h @@ -26,6 +26,12 @@ #include "PXR_Solver3D_GPSTD.h" #include "PXR_SolverAM_GPSTD.h" +#include "PML_Solver2D_Bouchard.h" +#include "PML_Solver2D_Yee.h" +#include "PML_Solver3D_Bouchard.h" +#include "PML_Solver3D_Yee.h" +#include "PML_SolverAM.h" + #include "Params.h" #include "Tools.h" @@ -153,7 +159,39 @@ class SolverFactory return solver; }; - + + + // Create PML solver + // ----------------------------- + static Solver *createPML( Params ¶ms ) + { + Solver *solver = NULL; + if( params.geometry == "2Dcartesian" ) { + if (params.maxwell_sol=="Bouchard"){ + solver = new PML_Solver2D_Bouchard( params ); + } + else {//(params.maxwell_sol=="Yee") + solver = new PML_Solver2D_Yee( params ); + } + } + else if( params.geometry == "3Dcartesian" ) { + if (params.maxwell_sol=="Bouchard"){ + solver = new PML_Solver3D_Bouchard( params ); + } + else { + solver = new PML_Solver3D_Yee( params ); + } + } + else if( params.geometry == "AMcylindrical" ) { + solver = new PML_SolverAM( params ); + } + else { + ERROR( "PML configuration not implemented" ); + } + + return solver; + } + }; #endif diff --git a/src/MovWindow/SimWindow.cpp b/src/MovWindow/SimWindow.cpp index 6e45ea498..e79b6dd1a 100755 --- a/src/MovWindow/SimWindow.cpp +++ b/src/MovWindow/SimWindow.cpp @@ -172,7 +172,7 @@ void SimWindow::shift( VectorPatch &vecPatches, SmileiMPI *smpi, Params ¶ms, int Href_receiver = 0; for (int irk = 0; irk < mypatch->MPI_neighbor_[0][0]; irk++) Href_receiver += smpi->patch_count[irk]; // The tag is the patch number in the receiver vector of patches in order to avoid too large tags not supported by some MPI versions. - smpi->isend( vecPatches_old[ipatch], vecPatches_old[ipatch]->MPI_neighbor_[0][0], ( vecPatches_old[ipatch]->neighbor_[0][0] - Href_receiver ) * nmessage, params ); + smpi->isend( vecPatches_old[ipatch], vecPatches_old[ipatch]->MPI_neighbor_[0][0], ( vecPatches_old[ipatch]->neighbor_[0][0] - Href_receiver ) * nmessage, params, false ); } } } else { //In case my left neighbor belongs to me: @@ -226,7 +226,7 @@ void SimWindow::shift( VectorPatch &vecPatches, SmileiMPI *smpi, Params ¶ms, if( mypatch->MPI_neighbor_[0][1] != MPI_PROC_NULL ) { if ( mypatch->Pcoordinates[0]!=params.number_of_patches[0]-1 ) { // The tag is the patch number in the receiver vector of patches in order to avoid too large tags not supported by some MPI versions. - smpi->recv( mypatch, mypatch->MPI_neighbor_[0][1], ( mypatch->hindex - vecPatches.refHindex_ )*nmessage, params ); + smpi->recv( mypatch, mypatch->MPI_neighbor_[0][1], ( mypatch->hindex - vecPatches.refHindex_ )*nmessage, params, false ); patch_particle_created[my_thread][j] = false ; //Mark no needs of particles } } diff --git a/src/Params/OpenPMDparams.cpp b/src/Params/OpenPMDparams.cpp index 4f655d989..1e313f21d 100755 --- a/src/Params/OpenPMDparams.cpp +++ b/src/Params/OpenPMDparams.cpp @@ -108,7 +108,7 @@ OpenPMDparams::OpenPMDparams( Params &p ): fieldBoundary .addString( "open" ); fieldBoundaryParameters.addString( params->EM_BCs[i][j] ); } else { - ERROR( " impossible boundary condition " ); + //ERROR( " impossible boundary condition " ); } particleBoundary .addString( "" ); particleBoundaryParameters.addString( "" ); diff --git a/src/Params/Params.cpp b/src/Params/Params.cpp index 32eeb3180..4c77e18ec 100755 --- a/src/Params/Params.cpp +++ b/src/Params/Params.cpp @@ -393,6 +393,10 @@ Params::Params( SmileiMPI *smpi, std::vector namelistsFiles ) : ERROR_NAMELIST( "EM_boundary_conditions along "<<"xyz"[iDim]<<" must be periodic for spectral solver in cartesian geometry.", LINK_NAMELIST + std::string("#main-variables") ); } + //if ( ( (EM_BCs[0][0] == "PML") || (EM_BCs[0][1] == "PML") ) + // && ( (EM_BCs[iDim][0] != "PML") || (EM_BCs[iDim][1] != "PML") ) ) { + // ERROR( "Either all PML, either none" ); + //} } int n_envlaser = PyTools::nComponents( "LaserEnvelope" ); @@ -506,6 +510,8 @@ Params::Params( SmileiMPI *smpi, std::vector namelistsFiles ) : save_magnectic_fields_for_SM = true; PyTools::extract( "save_magnectic_fields_for_SM", save_magnectic_fields_for_SM, "Main" ); + PyTools::extractVV( "number_of_pml_cells", number_of_pml_cells, "Main" ); + // ----------------------------------- // POISSON & FILTERING OPTIONS // ----------------------------------- diff --git a/src/Params/Params.h b/src/Params/Params.h index d0cf6baeb..534939f57 100755 --- a/src/Params/Params.h +++ b/src/Params/Params.h @@ -131,6 +131,7 @@ class Params //! Are open boundaries used ? std::vector< std::vector > open_boundaries; bool save_magnectic_fields_for_SM; + std::vector< std::vector > number_of_pml_cells; //! Boundary conditions for Envelope Field std::vector< std::vector > Env_BCs; diff --git a/src/Patch/Patch.cpp b/src/Patch/Patch.cpp index add6068b1..293132ee7 100755 --- a/src/Patch/Patch.cpp +++ b/src/Patch/Patch.cpp @@ -249,6 +249,15 @@ void Patch::finalizeMPIenvironment( Params ¶ms ) nb_comms += 18; } } + if ( (dynamic_cast( EMfields->emBoundCond[bcId] ))// && dynamic_cast( EMfields->emBoundCond[bcId] )->Hx_ ) + || + (dynamic_cast( EMfields->emBoundCond[bcId] ))) //&& dynamic_cast( EMfields->emBoundCond[bcId] )->Hx_ )) + { + nb_comms += 12; + } else if (dynamic_cast( EMfields->emBoundCond[bcId] ))// && dynamic_cast( EMfields->emBoundCond[bcId] )->Hl_[0] ){ + { + nb_comms += 12*( params.nmodes ); + } } requests_.resize( nb_comms, MPI_REQUEST_NULL ); @@ -361,7 +370,6 @@ void Patch::setLocationAndAllocateFields( Params ¶ms, DomainDecomposition *d Pcoordinates[0] =xDom ; Pcoordinates[1] =yDom; - //cout << "coords = " << Pcoordinates[0] << " " << Pcoordinates[1] < xcall( 2, 0 ); // 1st direction xcall[0] = Pcoordinates[0]-1; @@ -448,14 +455,6 @@ void Patch::setLocationAndAllocateFields( Params ¶ms, DomainDecomposition *d } neighbor_[1][1] = domain_decomposition->getDomainId( xcall ); - //cout << "\t"<< neighbor_[1][1] << endl; - //cout << neighbor_[0][0] << "\t h_me \t" << neighbor_[0][1] << endl; - //cout << "\t"<< neighbor_[1][0] << endl; - // - //cout << "\t"<< MPI_neighbor_[1][1] << endl; - //cout << MPI_neighbor_[0][0] << "\t mpi_me \t" << MPI_neighbor_[0][1] << endl; - //cout << "\t"<< MPI_neighbor_[1][0] << endl; - cell_starting_global_index[0] -= oversize[0]; cell_starting_global_index[1] -= oversize[1]; } // Fin 2D @@ -471,7 +470,6 @@ void Patch::setLocationAndAllocateFields( Params ¶ms, DomainDecomposition *d Pcoordinates[0] = xDom; Pcoordinates[1] = yDom; Pcoordinates[2] = zDom; - //cout << hindex << " - coords = " << xDom << " " << yDom << " " << zDom << endl; min_local_[0] = params.offset_map[0][xDom] * params.cell_length[0]; max_local_[0] = (params.offset_map[0][xDom]+params.n_space_region[0]) * params.cell_length[0]; @@ -509,7 +507,6 @@ void Patch::setLocationAndAllocateFields( Params ¶ms, DomainDecomposition *d MPI_neighbor_[0][1] = MPI_PROC_NULL; - //cout << MPI_neighbor_[0][0] << " " << " me " << " " << MPI_neighbor_[0][1] << endl; // ---------------- Y ---------------- if (yDom>0) @@ -664,11 +661,6 @@ void Patch::setLocationAndAllocateFields( Params ¶ms, DomainDecomposition *d } }*/ - - //cout << "MPI Nei\t" << "\t" << MPI_neighbor_[1][1] << endl; - //cout << "MPI Nei\t" << MPI_neighbor_[0][0] << "\t" << MPI_me_ << "\t" << MPI_neighbor_[0][1] << endl; - //cout << "MPI Nei\t" << "\t" << MPI_neighbor_[1][0] << endl; - EMfields = ElectroMagnFactory::create( params, domain_decomposition, vecPatch( 0 )->vecSpecies, this ); vecSpecies.resize( 0 ); @@ -1338,7 +1330,6 @@ void Patch::finalizeExchange( Field *field, int iDim ) { MPI_Status sstat [nDim_fields_][2]; MPI_Status rstat [nDim_fields_][2]; - for( int iNeighbor=0 ; iNeighborMPIbuff.srequest[iDim][iNeighbor] ), &( sstat[iDim][iNeighbor] ) ); diff --git a/src/Patch/SyncVectorPatch.cpp b/src/Patch/SyncVectorPatch.cpp index 70000d34b..914aff338 100755 --- a/src/Patch/SyncVectorPatch.cpp +++ b/src/Patch/SyncVectorPatch.cpp @@ -1592,3 +1592,606 @@ void SyncVectorPatch::finalizeExchangeAllComponentsAlongZ( std::vector } + +// --------------------------------------------------------------------------------------------------------------------- +// --------------------------------------------------------------------------------------------------------------------- +// ------------------------------------------- DEPRECATED FIELD FUNCTIONS -------------------------------------------- +// --------------------------------------------------------------------------------------------------------------------- +// --------------------------------------------------------------------------------------------------------------------- + +template +void SyncVectorPatch::exchangeAlongX( std::vector fields, VectorPatch &vecPatches, SmileiMPI *smpi ) +{ + unsigned oversize = vecPatches( 0 )->EMfields->oversize[0]; + +#ifndef _NO_MPI_TM + #pragma omp for schedule(static) +#else + #pragma omp single +#endif + for( unsigned int ipatch=0 ; ipatchis_a_MPI_neighbor( 0, iNeighbor ) ) { + fields[ipatch]->create_sub_fields ( 0, iNeighbor, oversize ); + fields[ipatch]->extract_fields_exch( 0, iNeighbor, oversize ); + } + } + if ( !dynamic_cast( fields[ipatch] ) ) + vecPatches( ipatch )->initExchange( fields[ipatch], 0, smpi ); + else + vecPatches( ipatch )->initExchangeComplex( fields[ipatch], 0, smpi ); + } + } + + unsigned int ny_( 1 ), nz_( 1 ), h0, neighbour_n_space, gsp; + T *pt1, *pt2; + h0 = vecPatches( 0 )->hindex; + + // fields[0] can be NULL, look for the 1st non null field + int IPATCH = 0; + for( unsigned int ipatch=0 ; ipatchisDual_[0] ); //Ghost size primal + + #pragma omp for schedule(static) private(pt1,pt2) + for( unsigned int ipatch=0 ; ipatchMPI_me_ == vecPatches( ipatch )->MPI_neighbor_[0][0] ) { + if (fields[ipatch]!=NULL) { + Field* f = fields[vecPatches( ipatch )->neighbor_[0][0]-h0]; + neighbour_n_space = f->dims_[0] - 1 - 2*oversize - f->isDual_[0]; // define here to be applied on usual or PML fields + if( fields[IPATCH]->dims_.size()>1 ) { + ny_ = f->dims_[1]; // Some can have ny, other ny+npml, but neighboors along X have the same number + if( fields[IPATCH]->dims_.size()>2 ) { + nz_ = f->dims_[2]; // Some can have nz, other nz+npml, but neighboors along X have the same number + } + } + pt1 = &( *static_cast( f ) )( neighbour_n_space*ny_*nz_ ); //my Xmin neighbour + pt2 = &( *static_cast( fields[ipatch] ) )( 0 ); // me + + memcpy( pt2, pt1, oversize*ny_*nz_*sizeof( T ) ); //me <= my neighbour + memcpy( pt1+gsp*ny_*nz_, pt2+gsp*ny_*nz_, oversize*ny_*nz_*sizeof( T ) );// my neighbour <= me + + + } // End if ( MPI_me_ == MPI_neighbor_[0][0] ) + } + } // End for( ipatch ) + +} + +void SyncVectorPatch::finalizeExchangeAlongX( std::vector fields, VectorPatch &vecPatches ) +{ + unsigned oversize = vecPatches( 0 )->EMfields->oversize[0]; + +#ifndef _NO_MPI_TM + #pragma omp for schedule(static) +#else + #pragma omp single +#endif + for( unsigned int ipatch=0 ; ipatchfinalizeExchange( fields[ipatch], 0 ); + for (int iNeighbor=0 ; iNeighbor<2 ; iNeighbor++) { + if ( vecPatches( ipatch )->is_a_MPI_neighbor( 0, ( iNeighbor+1 )%2 ) ) { + fields[ipatch]->inject_fields_exch( 0, iNeighbor, oversize ); + } + } + } + } + +} + +template +void SyncVectorPatch::exchangeAlongY( std::vector fields, VectorPatch &vecPatches, SmileiMPI *smpi ) +{ + unsigned oversize = vecPatches( 0 )->EMfields->oversize[1]; + +#ifndef _NO_MPI_TM + #pragma omp for schedule(static) +#else + #pragma omp single +#endif + for( unsigned int ipatch=0 ; ipatchis_a_MPI_neighbor( 1, iNeighbor ) ) { + fields[ipatch]->create_sub_fields ( 1, iNeighbor, oversize ); + fields[ipatch]->extract_fields_exch( 1, iNeighbor, oversize ); + } + } + if ( !dynamic_cast( fields[ipatch] ) ) + vecPatches( ipatch )->initExchange( fields[ipatch], 1, smpi ); + else + vecPatches( ipatch )->initExchangeComplex( fields[ipatch], 1, smpi ); + } + } + + unsigned int nx_, ny_, nz_( 1 ), h0, gsp, neighbour_n_space, my_ny; + T *pt1, *pt2; + h0 = vecPatches( 0 )->hindex; + + // fields[0] can be NULL, look for the 1st non null field + int IPATCH = 0; + for( unsigned int ipatch=0 ; ipatchisDual_[1] ); //Ghost size primal + + #pragma omp for schedule(static) private(pt1,pt2) + for( unsigned int ipatch=0 ; ipatchMPI_me_ == vecPatches( ipatch )->MPI_neighbor_[1][0] ) { + if (fields[ipatch]!=NULL) { + Field* f = fields[vecPatches( ipatch )->neighbor_[1][0]-h0]; + // In case of PML, ny can be different from e and my neighbour + nx_ = f->dims_[0]; // Some can have nx, other nx+npml, but neighbours along Y have the same number + ny_ = f->dims_[1]; // Neighbour ny_ + my_ny = fields[ipatch]->dims_[1]; + neighbour_n_space = ny_ - 1 - 2*oversize - f->isDual_[1]; // define here to be applied on usual or PML fields + if( f->dims_.size()>2 ) { + nz_ = f->dims_[2]; // Some can have nz, other nz+npml, but neighboors along Y have the same number + } + pt1 = &( *static_cast( f ) )( neighbour_n_space * nz_ );//my Ymin neighbour + pt2 = &( *static_cast( fields[ipatch] ) )( 0 );//me + for( unsigned int i = 0 ; i < nx_ ; i ++ ) { + for( unsigned int j = 0 ; j < oversize*nz_ ; j++ ) { + pt2[i*my_ny*nz_ + j] = pt1[i*ny_*nz_ + j] ; // me <= my_neighbour + pt1[i*ny_*nz_ + j + gsp*nz_] = pt2[i*my_ny*nz_ + j + gsp*nz_] ; //my_neighbour <= me + } + } + } + } + } + +} + +void SyncVectorPatch::finalizeExchangeAlongY( std::vector fields, VectorPatch &vecPatches ) +{ + unsigned oversize = vecPatches( 0 )->EMfields->oversize[1]; + +#ifndef _NO_MPI_TM + #pragma omp for schedule(static) +#else + #pragma omp single +#endif + for( unsigned int ipatch=0 ; ipatchfinalizeExchange( fields[ipatch], 1 ); + for (int iNeighbor=0 ; iNeighbor<2 ; iNeighbor++) { + if ( vecPatches( ipatch )->is_a_MPI_neighbor( 1, ( iNeighbor+1 )%2 ) ) { + fields[ipatch]->inject_fields_exch( 1, iNeighbor, oversize ); + } + } + } + } + +} + +template +void SyncVectorPatch::exchangeAlongZ( std::vector fields, VectorPatch &vecPatches, SmileiMPI *smpi ) +{ + unsigned oversize = vecPatches( 0 )->EMfields->oversize[2]; + +#ifndef _NO_MPI_TM + #pragma omp for schedule(static) +#else + #pragma omp single +#endif + for( unsigned int ipatch=0 ; ipatchis_a_MPI_neighbor( 2, iNeighbor ) ) { + fields[ipatch]->create_sub_fields ( 2, iNeighbor, oversize ); + fields[ipatch]->extract_fields_exch( 2, iNeighbor, oversize ); + } + } + if ( !dynamic_cast( fields[ipatch] ) ) + vecPatches( ipatch )->initExchange( fields[ipatch], 2, smpi ); + else + vecPatches( ipatch )->initExchangeComplex( fields[ipatch], 2, smpi ); + } + } + + unsigned int nx_, ny_, nz_, h0, gsp; + T *pt1, *pt2; + h0 = vecPatches( 0 )->hindex; + + // fields[0] can be NULL, look for the 1st non null field + int IPATCH = 0; + for( unsigned int ipatch=0 ; ipatchisDual_[2] ); //Ghost size primal + + #pragma omp for schedule(static) private(pt1,pt2) + for( unsigned int ipatch=0 ; ipatchMPI_me_ == vecPatches( ipatch )->MPI_neighbor_[2][0] ) { + if (fields[ipatch]!=NULL) { + Field* f = fields[vecPatches( ipatch )->neighbor_[2][0]-h0]; + nx_ = f->dims_[0]; // Some can have nx, other nx+npml, but neighboors along Z have the same number + ny_ = f->dims_[1]; // Some can have ny, other ny+npml, but neighboors along Z have the same number + nz_ = f->dims_[2]; // nz is same for all PML + pt1 = &( *static_cast( f ) )( nz_ - 1 - 2*oversize - f->isDual_[2] ); //my neighbour + pt2 = &( *static_cast( fields[ipatch] ) )( 0 ); //me + for( unsigned int i = 0 ; i < nx_*ny_*nz_ ; i += ny_*nz_ ) { + for( unsigned int j = 0 ; j < ny_*nz_ ; j += nz_ ) { + for( unsigned int k = 0 ; k < oversize ; k++ ) { + pt2[i+j+k] = pt1[i+j+k] ; //me <= my neighbour + pt1[i+j+k+gsp] = pt2[i+j+k+gsp] ; // my_neighbour <= me + } + } + } + } + } + } + +} + +void SyncVectorPatch::finalizeExchangeAlongZ( std::vector fields, VectorPatch &vecPatches ) +{ + unsigned oversize = vecPatches( 0 )->EMfields->oversize[2]; + +#ifndef _NO_MPI_TM + #pragma omp for schedule(static) +#else + #pragma omp single +#endif + for( unsigned int ipatch=0 ; ipatchfinalizeExchange( fields[ipatch], 2 ); + for (int iNeighbor=0 ; iNeighbor<2 ; iNeighbor++) { + if ( vecPatches( ipatch )->is_a_MPI_neighbor( 2, ( iNeighbor+1 )%2 ) ) { + fields[ipatch]->inject_fields_exch( 2, iNeighbor, oversize ); + } + } + } + } + +} + +void SyncVectorPatch::exchangeForPML( Params ¶ms, VectorPatch &vecPatches, SmileiMPI *smpi ) +{ + if ( ( params.geometry != "AMcylindrical" ) ) { + if (params.nDim_field==1) { + return; + } + else { + // Testing implementation of distributed PML on Xmin and Xmax + int iDim = 0; + for ( int min_max=0 ; min_max<2 ; min_max++ ) { + #pragma omp single + vecPatches.buildPMLList( "Bx", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongY( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "By", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongY( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Bz", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongY( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hx", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongY( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hy", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongY( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hz", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongY( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + } + if (params.nDim_field>2) { + // In 3D, distributed PML on Xmin and Xmax require synchronization along Z + for ( int min_max=0 ; min_max<2 ; min_max++ ) { + #pragma omp single + vecPatches.buildPMLList( "Bx", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongZ( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongZ( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "By", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongZ( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongZ( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Bz", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongZ( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongZ( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hx", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongZ( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongZ( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hy", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongZ( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongZ( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hz", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongZ( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongZ( vecPatches.listForPML_, vecPatches ); + } + } + // Testing implementation of distributed PML on Ymin and Ymax + iDim = 1; + for ( int min_max=0 ; min_max<2 ; min_max++ ) { + #pragma omp single + vecPatches.buildPMLList( "Bx", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongX( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "By", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongX( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Bz", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongX( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hx", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongX( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hy", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongX( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hz", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongX( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + } + if (params.nDim_field>2) { + // In 3D, distributed PML on Ymin and Ymax require synchronization along Z + for ( int min_max=0 ; min_max<2 ; min_max++ ) { + #pragma omp single + vecPatches.buildPMLList( "Bx", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongZ( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongZ( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "By", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongZ( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongZ( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Bz", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongZ( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongZ( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hx", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongZ( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongZ( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hy", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongZ( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongZ( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hz", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongZ( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongZ( vecPatches.listForPML_, vecPatches ); + } + } + + if (params.nDim_field>2) { + int iDim = 2; + for ( int min_max=0 ; min_max<2 ; min_max++ ) { + #pragma omp single + vecPatches.buildPMLList( "Bx", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongX( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "By", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongX( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Bz", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongX( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hx", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongX( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hy", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongX( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hz", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongX( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Bx", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongY( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "By", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongY( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Bz", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongY( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hx", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongY( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hy", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongY( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hz", iDim, min_max, smpi ); + + SyncVectorPatch::exchangeAlongY( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + } + } // End if (ndim_field==1) + } // End if (ndim_field>1) + } // End if (cartesian) + else { // AM + for (unsigned int imode = 0 ; imode < params.nmodes ; imode++) { + // Testing implementation of distributed PML on Xmin and Xmax + int iDim = 0; + for ( int min_max=0 ; min_max<2 ; min_max++ ) { + #pragma omp single + vecPatches.buildPMLList( "Bl", iDim, min_max, smpi, imode ); + + SyncVectorPatch::exchangeAlongY,cField>( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Br", iDim, min_max, smpi, imode ); + + SyncVectorPatch::exchangeAlongY,cField>( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Bt", iDim, min_max, smpi, imode ); + + SyncVectorPatch::exchangeAlongY,cField>( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hl", iDim, min_max, smpi, imode ); + + SyncVectorPatch::exchangeAlongY,cField>( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hr", iDim, min_max, smpi, imode ); + + SyncVectorPatch::exchangeAlongY,cField>( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Ht", iDim, min_max, smpi, imode ); + + SyncVectorPatch::exchangeAlongY,cField>( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongY( vecPatches.listForPML_, vecPatches ); + } + + // Testing implementation of distributed PML on Ymin and Ymax + iDim = 1; + for ( int min_max=0 ; min_max<2 ; min_max++ ) { + #pragma omp single + vecPatches.buildPMLList( "Bl", iDim, min_max, smpi, imode ); + + SyncVectorPatch::exchangeAlongX,cField>( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Br", iDim, min_max, smpi, imode ); + + SyncVectorPatch::exchangeAlongX,cField>( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Bt", iDim, min_max, smpi, imode ); + + SyncVectorPatch::exchangeAlongX,cField>( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hl", iDim, min_max, smpi, imode ); + + SyncVectorPatch::exchangeAlongX,cField>( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Hr", iDim, min_max, smpi, imode ); + + SyncVectorPatch::exchangeAlongX,cField>( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + + #pragma omp single + vecPatches.buildPMLList( "Ht", iDim, min_max, smpi, imode ); + + SyncVectorPatch::exchangeAlongX,cField>( vecPatches.listForPML_, vecPatches, smpi ); + SyncVectorPatch::finalizeExchangeAlongX( vecPatches.listForPML_, vecPatches ); + } + } // End for( imode ) + } // End else( AM ) +} diff --git a/src/Patch/SyncVectorPatch.h b/src/Patch/SyncVectorPatch.h index 20174cfd7..8d9ce6e0d 100755 --- a/src/Patch/SyncVectorPatch.h +++ b/src/Patch/SyncVectorPatch.h @@ -331,6 +331,16 @@ public : static void exchangeAllComponentsAlongZ( std::vector fields, VectorPatch &vecPatches, SmileiMPI *smpi ); static void finalizeExchangeAllComponentsAlongZ( std::vector fields, VectorPatch &vecPatches ); + //! Deprecated field functions + template static void exchangeAlongX( std::vector fields, VectorPatch &vecPatches, SmileiMPI *smpi ); + static void finalizeExchangeAlongX( std::vector fields, VectorPatch &vecPatches ); + template static void exchangeAlongY( std::vector fields, VectorPatch &vecPatches, SmileiMPI *smpi ); + static void finalizeExchangeAlongY( std::vector fields, VectorPatch &vecPatches ); + template static void exchangeAlongZ( std::vector fields, VectorPatch &vecPatches, SmileiMPI *smpi ); + static void finalizeExchangeAlongZ( std::vector fields, VectorPatch &vecPatches ); + + static void exchangeForPML( Params ¶ms, VectorPatch &vecPatches, SmileiMPI *smpi ); + }; #endif diff --git a/src/Patch/VectorPatch.cpp b/src/Patch/VectorPatch.cpp index 344313774..2ba09f030 100755 --- a/src/Patch/VectorPatch.cpp +++ b/src/Patch/VectorPatch.cpp @@ -21,6 +21,10 @@ #include "ElectroMagnBC.h" #include "Laser.h" +#include "ElectroMagnBC2D_PML.h" +#include "ElectroMagnBC3D_PML.h" +#include "ElectroMagnBCAM_PML.h" + #include "SyncVectorPatch.h" #include "interface.h" #include "Timers.h" @@ -1271,6 +1275,12 @@ void VectorPatch::solveMaxwell( Params ¶ms, SimWindow *simWindow, int itime, for( unsigned int ipatch=0 ; ipatchsize() ; ipatch++ ) { // Applies boundary conditions on B ( *this )( ipatch )->EMfields->boundaryConditions( itime, time_dual, ( *this )( ipatch ), params, simWindow ); + } + if ( params.EM_BCs[0][0] == "PML" ) { // If a PML on 1 border, then on all + SyncVectorPatch::exchangeForPML( params, (*this), smpi ); + } + #pragma omp for schedule(static) + for( unsigned int ipatch=0 ; ipatchsize() ; ipatch++ ) { // Computes B at time n using B and B_m. if( !params.is_spectral ) { ( *this )( ipatch )->EMfields->centerMagneticFields(); @@ -1353,6 +1363,14 @@ void VectorPatch::finalizeSyncAndBCFields( Params ¶ms, SmileiMPI *smpi, SimW // Applies boundary conditions on B if ( (!params.is_spectral) || (params.geometry!= "AMcylindrical") ) ( *this )( ipatch )->EMfields->boundaryConditions( itime, time_dual, ( *this )( ipatch ), params, simWindow ); + + } + if ( params.EM_BCs[0][0] == "PML" ) { // If a PML on 1 border, then on all + SyncVectorPatch::exchangeForPML( params, (*this), smpi ); + } + + #pragma omp for schedule(static) + for( unsigned int ipatch=0 ; ipatchsize() ; ipatch++ ) { // Computes B at time n using B and B_m. if( !params.is_spectral ) { ( *this )( ipatch )->EMfields->centerMagneticFields(); @@ -3745,6 +3763,123 @@ void VectorPatch::updateFieldList( int ispec, SmileiMPI *smpi ) } +void VectorPatch::buildPMLList( string fieldname, int idim, int min_or_max, SmileiMPI *smpi ) +{ + // min_or_max = 0 : min + // min_or_max = 1 : max + int id_bc = 2*idim + min_or_max; + + listForPML_.clear(); + if ( fieldname == "Bx" ) { + for ( unsigned int ipatch=0 ; ipatch < size() ; ipatch++ ) { + listForPML_.push_back( ( emfields(ipatch)->emBoundCond[id_bc] )->getBxPML() ); + if(listForPML_.back()){ + listForPML_.back()->MPIbuff.defineTags(patches_[ipatch], smpi, 0); + } + } + } + else if ( fieldname == "By" ) { + for ( unsigned int ipatch=0 ; ipatch < size() ; ipatch++ ) { + listForPML_.push_back( ( emfields(ipatch)->emBoundCond[id_bc] )->getByPML() ); + if(listForPML_.back()){ + listForPML_.back()->MPIbuff.defineTags(patches_[ipatch], smpi, 0); + } + } + } + else if ( fieldname == "Bz" ) { + for ( unsigned int ipatch=0 ; ipatch < size() ; ipatch++ ) { + listForPML_.push_back( ( emfields(ipatch)->emBoundCond[id_bc] )->getBzPML() ); + if(listForPML_.back()){ + listForPML_.back()->MPIbuff.defineTags(patches_[ipatch], smpi, 0); + } + } + } + else if ( fieldname == "Hx" ) { + for ( unsigned int ipatch=0 ; ipatch < size() ; ipatch++ ) { + listForPML_.push_back( ( emfields(ipatch)->emBoundCond[id_bc] )->getHxPML() ); + if(listForPML_.back()){ + listForPML_.back()->MPIbuff.defineTags(patches_[ipatch], smpi, 0); + } + } + } + else if ( fieldname == "Hy" ) { + for ( unsigned int ipatch=0 ; ipatch < size() ; ipatch++ ) { + listForPML_.push_back( ( emfields(ipatch)->emBoundCond[id_bc] )->getHyPML() ); + if(listForPML_.back()){ + listForPML_.back()->MPIbuff.defineTags(patches_[ipatch], smpi, 0); + } + } + } + else if ( fieldname == "Hz" ) { + for ( unsigned int ipatch=0 ; ipatch < size() ; ipatch++ ) { + listForPML_.push_back( ( emfields(ipatch)->emBoundCond[id_bc] )->getHzPML() ); + if(listForPML_.back()){ + listForPML_.back()->MPIbuff.defineTags(patches_[ipatch], smpi, 0); + } + } + } +} + + +void VectorPatch::buildPMLList( string fieldname, int idim, int min_or_max, SmileiMPI *smpi, int imode ) +{ + // min_or_max = 0 : min + // min_or_max = 1 : max + int id_bc = 2*idim + min_or_max; + + listForPML_.clear(); + if ( fieldname == "Bl" ) { + for ( unsigned int ipatch=0 ; ipatch < size() ; ipatch++ ) { + listForPML_.push_back( static_cast( emfields(ipatch)->emBoundCond[id_bc] )->Bl_[imode] ); + //After pushing back a pointer to PML fields, if this field is not NULL recompute the tags. + if(listForPML_.back()){ + listForPML_.back()->MPIbuff.defineTags(patches_[ipatch], smpi, 0); + } + } + } + else if ( fieldname == "Br" ) { + for ( unsigned int ipatch=0 ; ipatch < size() ; ipatch++ ) { + listForPML_.push_back( static_cast( emfields(ipatch)->emBoundCond[id_bc] )->Br_[imode] ); + if(listForPML_.back()){ + listForPML_.back()->MPIbuff.defineTags(patches_[ipatch], smpi, 0); + } + } + } + else if ( fieldname == "Bt" ) { + for ( unsigned int ipatch=0 ; ipatch < size() ; ipatch++ ) { + listForPML_.push_back( static_cast( emfields(ipatch)->emBoundCond[id_bc] )->Bt_[imode] ); + if(listForPML_.back()){ + listForPML_.back()->MPIbuff.defineTags(patches_[ipatch], smpi, 0); + } + } + } + else if ( fieldname == "Hl" ) { + for ( unsigned int ipatch=0 ; ipatch < size() ; ipatch++ ) { + listForPML_.push_back( static_cast( emfields(ipatch)->emBoundCond[id_bc] )->Hl_[imode] ); + if(listForPML_.back()){ + listForPML_.back()->MPIbuff.defineTags(patches_[ipatch], smpi, 0); + } + } + } + else if ( fieldname == "Hr" ) { + for ( unsigned int ipatch=0 ; ipatch < size() ; ipatch++ ) { + listForPML_.push_back( static_cast( emfields(ipatch)->emBoundCond[id_bc] )->Hr_[imode] ); + if(listForPML_.back()){ + listForPML_.back()->MPIbuff.defineTags(patches_[ipatch], smpi, 0); + } + } + } + else if ( fieldname == "Ht" ) { + for ( unsigned int ipatch=0 ; ipatch < size() ; ipatch++ ) { + listForPML_.push_back( static_cast( emfields(ipatch)->emBoundCond[id_bc] )->Ht_[imode] ); + if(listForPML_.back()){ + listForPML_.back()->MPIbuff.defineTags(patches_[ipatch], smpi, 0); + } + } + } +} + + void VectorPatch::applyAntennas( double time ) { #ifdef __DEBUG diff --git a/src/Patch/VectorPatch.h b/src/Patch/VectorPatch.h index 2a6bac04c..808a3b877 100755 --- a/src/Patch/VectorPatch.h +++ b/src/Patch/VectorPatch.h @@ -74,7 +74,10 @@ public : //! Resize vector of field* void updateFieldList( SmileiMPI *smpi ); void updateFieldList( int ispec, SmileiMPI *smpi ); - + + void buildPMLList( std::string fieldname, int idim, int min_or_max, SmileiMPI *smpi ); + void buildPMLList( std::string fieldname, int idim, int min_or_max, SmileiMPI *smpi, int imode ); + //! Create the diagnostic list void createDiags( Params ¶ms, SmileiMPI *smpi, OpenPMDparams &, RadiationTables * radiation_tables_ ); @@ -282,6 +285,7 @@ public : std::vector listBx_; std::vector listBy_; std::vector listBz_; + std::vector listForPML_; std::vector listA_; std::vector listA0_; diff --git a/src/Python/pyinit.py b/src/Python/pyinit.py index 74606cdd8..83acbea85 100755 --- a/src/Python/pyinit.py +++ b/src/Python/pyinit.py @@ -197,6 +197,7 @@ class Main(SmileiSingleton): EM_boundary_conditions = [["periodic"]] EM_boundary_conditions_k = [] save_magnectic_fields_for_SM = True + number_of_pml_cells = [[10,10],[10,10],[10,10]] time_fields_frozen = 0. Laser_Envelope_model = False diff --git a/src/Radiation/RadiationCorrLandauLifshitz.cpp b/src/Radiation/RadiationCorrLandauLifshitz.cpp index b65c8a6d7..6044c79ad 100755 --- a/src/Radiation/RadiationCorrLandauLifshitz.cpp +++ b/src/Radiation/RadiationCorrLandauLifshitz.cpp @@ -190,6 +190,6 @@ void RadiationCorrLandauLifshitz::operator()( // _______________________________________________________________ // Cleaning - delete [] rad_norm_energy; + free(rad_norm_energy); } diff --git a/src/Radiation/RadiationLandauLifshitz.cpp b/src/Radiation/RadiationLandauLifshitz.cpp index 443db9dab..ae15bb2a7 100755 --- a/src/Radiation/RadiationLandauLifshitz.cpp +++ b/src/Radiation/RadiationLandauLifshitz.cpp @@ -187,6 +187,6 @@ void RadiationLandauLifshitz::operator()( // _______________________________________________________________ // Cleaning - delete [] rad_norm_energy; + free(rad_norm_energy); } diff --git a/src/SmileiMPI/SmileiMPI.cpp b/src/SmileiMPI/SmileiMPI.cpp index 53cb2312d..c533383fc 100755 --- a/src/SmileiMPI/SmileiMPI.cpp +++ b/src/SmileiMPI/SmileiMPI.cpp @@ -15,6 +15,9 @@ #include "ElectroMagnBC1D_SM.h" #include "ElectroMagnBC2D_SM.h" #include "ElectroMagnBC3D_SM.h" +#include "ElectroMagnBC2D_PML.h" +#include "ElectroMagnBC3D_PML.h" +#include "ElectroMagnBCAM_PML.h" #include "Field.h" #include "Species.h" @@ -664,7 +667,7 @@ MPI_Datatype SmileiMPI::createMPIparticles( Particles *particles ) // ----------------------------------------- PATCH SEND / RECV METHODS ------------------------------------ // --------------------------------------------------------------------------------------------------------------------- // --------------------------------------------------------------------------------------------------------------------- -void SmileiMPI::isend( Patch *patch, int to, int tag, Params ¶ms ) +void SmileiMPI::isend( Patch *patch, int to, int tag, Params ¶ms, bool send_xmax_bc ) { //MPI_Request request; @@ -672,7 +675,7 @@ void SmileiMPI::isend( Patch *patch, int to, int tag, Params ¶ms ) int irequest = 0; isend_species( patch, to, irequest, tag, params ); - isend_fields ( patch, to, irequest, tag, params ); + isend_fields ( patch, to, irequest, tag, params, send_xmax_bc ); } // END isend( Patch ) @@ -727,13 +730,13 @@ void SmileiMPI::isend_species( Patch *patch, int to, int &irequest, int tag, Par irequest ++; } -void SmileiMPI::isend_fields( Patch *patch, int to, int &irequest, int tag, Params ¶ms ) +void SmileiMPI::isend_fields( Patch *patch, int to, int &irequest, int tag, Params ¶ms, bool send_xmax_bc ) { // Send fields if( params.geometry != "AMcylindrical" ) { - isend( patch->EMfields, to, irequest, patch->requests_, tag ); + isend( patch->EMfields, to, irequest, patch->requests_, tag, send_xmax_bc ); } else { - isend( patch->EMfields, to, irequest, patch->requests_, tag, static_cast( patch->EMfields )->El_.size() ); + isend( patch->EMfields, to, irequest, patch->requests_, tag, static_cast( patch->EMfields )->El_.size(), send_xmax_bc ); } // Send some scalars @@ -774,13 +777,13 @@ void SmileiMPI::waitall( Patch *patch ) } -void SmileiMPI::recv( Patch *patch, int from, int tag, Params ¶ms ) +void SmileiMPI::recv( Patch *patch, int from, int tag, Params ¶ms, bool recv_xmin_bc ) { // Receive species recv_species( patch, from, tag, params ); // Receive EM fields - recv_fields( patch, from, tag, params ); + recv_fields( patch, from, tag, params, recv_xmin_bc ); } // END recv ( Patch ) @@ -854,14 +857,14 @@ void SmileiMPI::recv_species( Patch *patch, int from, int &tag, Params ¶ms ) } } -void SmileiMPI::recv_fields( Patch *patch, int from, int &tag, Params ¶ms ) +void SmileiMPI::recv_fields( Patch *patch, int from, int &tag, Params ¶ms, bool recv_xmin_bc ) { // Receive EM fields patch->EMfields->initAntennas( patch, params ); if( params.geometry != "AMcylindrical" ) { - recv( patch->EMfields, from, tag ); + recv( patch->EMfields, from, tag, recv_xmin_bc ); } else { - recv( patch->EMfields, from, tag, static_cast( patch->EMfields )->El_.size() ); + recv( patch->EMfields, from, tag, static_cast( patch->EMfields )->El_.size(), recv_xmin_bc ); } // Receive some scalars @@ -924,7 +927,7 @@ void SmileiMPI::recv( std::vector *vec, int from, int tag ) } // End recv -void SmileiMPI::isend( ElectroMagn *EM, int to, int &irequest, vector &requests, int tag ) +void SmileiMPI::isend( ElectroMagn *EM, int to, int &irequest, vector &requests, int tag, bool send_xmax_bc ) { isend( EM->Ex_, to, tag+irequest, requests[irequest] ); @@ -1054,11 +1057,69 @@ void SmileiMPI::isend( ElectroMagn *EM, int to, int &irequest, vector( EM->emBoundCond[bcId] ) || dynamic_cast( EM->emBoundCond[bcId] )) && (bcId != 1 || send_xmax_bc) ){ + if( dynamic_cast( EM->emBoundCond[bcId] )){ + ElectroMagnBC2D_PML *embc = static_cast( EM->emBoundCond[bcId] ); + if (embc->Hx_) { + isend( embc->Hx_, to, tag + irequest, requests[tag] ); + irequest++; + isend( embc->Hy_, to, tag + irequest, requests[tag] ); + irequest++; + isend( embc->Hz_, to, tag + irequest, requests[tag] ); + irequest++; + isend( embc->Bx_, to, tag + irequest, requests[tag] ); + irequest++; + isend( embc->By_, to, tag + irequest, requests[tag] ); + irequest++; + isend( embc->Bz_, to, tag + irequest, requests[tag] ); + irequest++; + isend( embc->Ex_, to, tag + irequest, requests[tag] ); + irequest++; + isend( embc->Ey_, to, tag + irequest, requests[tag] ); + irequest++; + isend( embc->Ez_, to, tag + irequest, requests[tag] ); + irequest++; + isend( embc->Dx_, to, tag + irequest, requests[tag] ); + irequest++; + isend( embc->Dy_, to, tag + irequest, requests[tag] ); + irequest++; + isend( embc->Dz_, to, tag + irequest, requests[tag] ); + irequest++; + } + } else { + ElectroMagnBC3D_PML *embc = static_cast( EM->emBoundCond[bcId] ); + if (embc->Hx_) { + isend( embc->Hx_, to, tag+irequest, requests[tag] ); + irequest++; + isend( embc->Hy_, to, tag+irequest, requests[tag] ); + irequest++; + isend( embc->Hz_, to, tag+irequest, requests[tag] ); + irequest++; + isend( embc->Bx_, to, tag+irequest, requests[tag] ); + irequest++; + isend( embc->By_, to, tag+irequest, requests[tag] ); + irequest++; + isend( embc->Bz_, to, tag+irequest, requests[tag] ); + irequest++; + isend( embc->Ex_, to, tag+irequest, requests[tag] ); + irequest++; + isend( embc->Ey_, to, tag+irequest, requests[tag] ); + irequest++; + isend( embc->Ez_, to, tag+irequest, requests[tag] ); + irequest++; + isend( embc->Dx_, to, tag+irequest, requests[tag] ); + irequest++; + isend( embc->Dy_, to, tag+irequest, requests[tag] ); + irequest++; + isend( embc->Dz_, to, tag+irequest, requests[tag] ); + irequest++; + } + } + } } } // End isend ( ElectroMagn ) -void SmileiMPI::isend( ElectroMagn *EM, int to, int &irequest, vector &requests, int tag, unsigned int nmodes ) +void SmileiMPI::isend( ElectroMagn *EM, int to, int &irequest, vector &requests, int tag, unsigned int nmodes, bool send_xmax_bc ) { ElectroMagnAM *EMAM = static_cast( EM ); @@ -1143,56 +1204,44 @@ void SmileiMPI::isend( ElectroMagn *EM, int to, int &irequest, vectorextFields.size()>0 ) { - if( dynamic_cast( EM->emBoundCond[bcId] ) ) { - ElectroMagnBC1D_SM *embc = static_cast( EM->emBoundCond[bcId] ); - MPI_Isend( &( embc->By_val ), 1, MPI_DOUBLE, to, tag+irequest, MPI_COMM_WORLD, &requests[irequest] ); - irequest++; - MPI_Isend( &( embc->Bz_val ), 1, MPI_DOUBLE, to, tag+irequest, MPI_COMM_WORLD, &requests[irequest] ); - irequest++; - } else if( dynamic_cast( EM->emBoundCond[bcId] ) ) { - // BCs at the x-border - ElectroMagnBC2D_SM *embc = static_cast( EM->emBoundCond[bcId] ); - - if( embc->B_val[0].size() ) { - isend( &embc->B_val[0], to, tag+irequest, requests[irequest] ); + if( dynamic_cast( EM->emBoundCond[bcId] ) && (bcId != 1 || send_xmax_bc) ){ + ElectroMagnBCAM_PML *embc = static_cast( EM->emBoundCond[bcId] ); + // if I have PML && if the receiver also have PMLs <=> the hindex I send to touches the same boundary I am dealing with now + if (embc->Hl_[0] ) { + for( unsigned int imode =0; imode < nmodes; imode++ ) { + isendComplex( embc->Hl_[imode], to, tag+irequest, requests[irequest] ); irequest++; - } - if( embc->B_val[1].size() ) { - isend( &embc->B_val[1], to, tag+irequest, requests[irequest] ); + isendComplex( embc->Hr_[imode], to, tag+irequest, requests[irequest] ); irequest++; - } - if( embc->B_val[2].size() ) { - isend( &embc->B_val[2], to, tag+irequest, requests[irequest] ); + isendComplex( embc->Ht_[imode], to, tag+irequest, requests[irequest] ); irequest++; - } - - } else if( dynamic_cast( EM->emBoundCond[bcId] ) ) { - ElectroMagnBC3D_SM *embc = static_cast( EM->emBoundCond[bcId] ); - - // BCs at the border - if( embc->B_val[0] ) { - isend( embc->B_val[0], to, tag+irequest, requests[irequest] ); + isendComplex( embc->Bl_[imode], to, tag+irequest, requests[irequest] ); irequest++; - } - if( embc->B_val[1] ) { - isend( embc->B_val[1], to, tag+irequest, requests[irequest] ); + isendComplex( embc->Br_[imode], to, tag+irequest, requests[irequest] ); irequest++; - } - if( embc->B_val[2] ) { - isend( embc->B_val[2], to, tag+irequest, requests[irequest] ); + isendComplex( embc->Bt_[imode], to, tag+irequest, requests[irequest] ); + irequest++; + isendComplex( embc->El_[imode], to, tag+irequest, requests[irequest] ); + irequest++; + isendComplex( embc->Er_[imode], to, tag+irequest, requests[irequest] ); + irequest++; + isendComplex( embc->Et_[imode], to, tag+irequest, requests[irequest] ); + irequest++; + isendComplex( embc->Dl_[imode], to, tag+irequest, requests[irequest] ); + irequest++; + isendComplex( embc->Dr_[imode], to, tag+irequest, requests[irequest] ); + irequest++; + isendComplex( embc->Dt_[imode], to, tag+irequest, requests[irequest] ); irequest++; } - } } - } } // End isend ( ElectroMagn LRT ) -void SmileiMPI::recv( ElectroMagn *EM, int from, int &tag ) +void SmileiMPI::recv( ElectroMagn *EM, int from, int &tag, bool recv_xmin_bc ) { recv( EM->Ex_, from, tag ); tag++; @@ -1318,12 +1367,70 @@ void SmileiMPI::recv( ElectroMagn *EM, int from, int &tag ) } } - + if( (dynamic_cast( EM->emBoundCond[bcId] ) || dynamic_cast( EM->emBoundCond[bcId] )) && (bcId != 0 || recv_xmin_bc) ){ + if( dynamic_cast( EM->emBoundCond[bcId] )){ + ElectroMagnBC2D_PML *embc = static_cast( EM->emBoundCond[bcId] ); + if (embc->Hx_) { + recv( embc->Hx_, from, tag ); + tag++; + recv( embc->Hy_, from, tag ); + tag++; + recv( embc->Hz_, from, tag ); + tag++; + recv( embc->Bx_, from, tag ); + tag++; + recv( embc->By_, from, tag ); + tag++; + recv( embc->Bz_, from, tag ); + tag++; + recv( embc->Ex_, from, tag ); + tag++; + recv( embc->Ey_, from, tag ); + tag++; + recv( embc->Ez_, from, tag ); + tag++; + recv( embc->Dx_, from, tag ); + tag++; + recv( embc->Dy_, from, tag ); + tag++; + recv( embc->Dz_, from, tag ); + tag++; + } + } else { + ElectroMagnBC3D_PML *embc = static_cast( EM->emBoundCond[bcId] ); + if (embc->Hx_) { + recv( embc->Hx_, from, tag ); + tag++; + recv( embc->Hy_, from, tag ); + tag++; + recv( embc->Hz_, from, tag ); + tag++; + recv( embc->Bx_, from, tag ); + tag++; + recv( embc->By_, from, tag ); + tag++; + recv( embc->Bz_, from, tag ); + tag++; + recv( embc->Ex_, from, tag ); + tag++; + recv( embc->Ey_, from, tag ); + tag++; + recv( embc->Ez_, from, tag ); + tag++; + recv( embc->Dx_, from, tag ); + tag++; + recv( embc->Dy_, from, tag ); + tag++; + recv( embc->Dz_, from, tag ); + tag++; + } + } + } } } // End recv ( ElectroMagn ) -void SmileiMPI::recv( ElectroMagn *EM, int from, int &tag, unsigned int nmodes ) +void SmileiMPI::recv( ElectroMagn *EM, int from, int &tag, unsigned int nmodes, bool recv_xmin_bc ) { ElectroMagnAM *EMAM = static_cast( EM ); for( unsigned int imode =0; imode < nmodes; imode++ ) { @@ -1403,49 +1510,35 @@ void SmileiMPI::recv( ElectroMagn *EM, int from, int &tag, unsigned int nmodes ) } } - if( EM->extFields.size()>0 ) { - - if( dynamic_cast( EM->emBoundCond[bcId] ) ) { - ElectroMagnBC1D_SM *embc = static_cast( EM->emBoundCond[bcId] ); - MPI_Status status; - MPI_Recv( &( embc->By_val ), 1, MPI_DOUBLE, from, tag, MPI_COMM_WORLD, &status ); - tag++; - MPI_Recv( &( embc->Bz_val ), 1, MPI_DOUBLE, from, tag, MPI_COMM_WORLD, &status ); - tag++; - } else if( dynamic_cast( EM->emBoundCond[bcId] ) ) { - // BCs at the x-border - ElectroMagnBC2D_SM *embc = static_cast( EM->emBoundCond[bcId] ); - - if( embc->B_val[0].size() ) { - recv( &embc->B_val[0], from, tag ); + if( dynamic_cast( EM->emBoundCond[bcId] ) && (bcId != 0 || recv_xmin_bc)){ + ElectroMagnBCAM_PML *embc = static_cast( EM->emBoundCond[bcId] ); + if (embc->Hl_[0]) { + for( unsigned int imode =0; imode < nmodes; imode++ ) { + recvComplex( embc->Hl_[imode], from, tag ); tag++; - } - if( embc->B_val[1].size() ) { - recv( &embc->B_val[1], from, tag ); + recvComplex( embc->Hr_[imode], from, tag ); tag++; - } - if( embc->B_val[2].size() ) { - recv( &embc->B_val[2], from, tag ); + recvComplex( embc->Ht_[imode], from, tag ); tag++; - } - - } else if( dynamic_cast( EM->emBoundCond[bcId] ) ) { - ElectroMagnBC3D_SM *embc = static_cast( EM->emBoundCond[bcId] ); - - // BCs at the border - if( embc->B_val[0] ) { - recv( embc->B_val[0], from, tag ); + recvComplex( embc->Bl_[imode], from, tag ); tag++; - } - if( embc->B_val[1] ) { - recv( embc->B_val[1], from, tag ); + recvComplex( embc->Br_[imode], from, tag ); tag++; - } - if( embc->B_val[2] ) { - recv( embc->B_val[2], from, tag ); + recvComplex( embc->Bt_[imode], from, tag ); + tag++; + recvComplex( embc->El_[imode], from, tag ); + tag++; + recvComplex( embc->Er_[imode], from, tag ); + tag++; + recvComplex( embc->Et_[imode], from, tag ); + tag++; + recvComplex( embc->Dl_[imode], from, tag ); + tag++; + recvComplex( embc->Dr_[imode], from, tag ); + tag++; + recvComplex( embc->Dt_[imode], from, tag ); tag++; } - } } diff --git a/src/SmileiMPI/SmileiMPI.h b/src/SmileiMPI/SmileiMPI.h index e62ff275d..e1faecbee 100755 --- a/src/SmileiMPI/SmileiMPI.h +++ b/src/SmileiMPI/SmileiMPI.h @@ -68,12 +68,12 @@ class SmileiMPI // - during load balancing process // - during moving window // ----------------------------------- - void isend( Patch *patch, int to, int tag, Params ¶ms ); + void isend( Patch *patch, int to, int tag, Params ¶ms, bool send_xmax_bc = true ); void waitall( Patch *patch ); - void recv( Patch *patch, int from, int tag, Params ¶ms ); + void recv( Patch *patch, int from, int tag, Params ¶ms, bool recv_xmin_bc = true); - void isend_fields( Patch *patch, int to, int &irequest, int tag, Params ¶ms ); - void recv_fields( Patch *patch, int from, int &tag, Params ¶ms ); + void isend_fields( Patch *patch, int to, int &irequest, int tag, Params ¶ms, bool send_xmax_bc = true ); + void recv_fields( Patch *patch, int from, int &tag, Params ¶ms, bool recv_xmin_bc = true ); void isend_species( Patch *patch, int to, int &irequest, int tag, Params ¶ms ); void recv_species( Patch *patch, int from, int &tag, Params ¶ms ); @@ -85,10 +85,10 @@ class SmileiMPI void isend( std::vector *vec, int to, int tag, MPI_Request &request ); void recv( std::vector *vec, int from, int tag ); - void isend( ElectroMagn *fields, int to, int &irequest, std::vector &requests, int tag ); - void isend( ElectroMagn *fields, int to, int &irequest, std::vector &requests, int tag, unsigned int nmodes ); - void recv( ElectroMagn *fields, int from, int &tag ); - void recv( ElectroMagn *fields, int from, int &tag, unsigned int nmodes ); + void isend( ElectroMagn *fields, int to, int &irequest, std::vector &requests, int tag, bool send_xmax_bc ); + void isend( ElectroMagn *fields, int to, int &irequest, std::vector &requests, int tag, unsigned int nmodes, bool send_xmax_bc ); + void recv( ElectroMagn *fields, int from, int &tag, bool recv_xmax_bc ); + void recv( ElectroMagn *fields, int from, int &tag, unsigned int nmodes, bool recv_xmax_bc ); void isend( Field *field, int to, int tag, MPI_Request &request ); void isendComplex( Field *field, int to, int tag, MPI_Request &request ); void recv( Field *field, int from, int tag ); diff --git a/validation/analyses/validate_tst2d_18_em_pml.py b/validation/analyses/validate_tst2d_18_em_pml.py new file mode 100755 index 000000000..a2f347581 --- /dev/null +++ b/validation/analyses/validate_tst2d_18_em_pml.py @@ -0,0 +1,15 @@ +import os, re, numpy as np, math, h5py +import happi + +S = happi.Open(["./restart*"], verbose=False) + + + +# COMPARE THE Ez FIELD +Ez = S.Field.Field0.Ez(timesteps=300).getData()[0] +Validate("Ez field at iteration 300", Ez, 0.001) + +# TEST THAT Ubal_norm STAYS OK +uelm = S.Scalar.Uelm().getData(timestep=300) +Validate("Uelm is below 0.2% after absorption", uelm[0]<0.002 ) + diff --git a/validation/references/tst2d_18_em_pml.py.txt b/validation/references/tst2d_18_em_pml.py.txt new file mode 100644 index 000000000..248ef961e Binary files /dev/null and b/validation/references/tst2d_18_em_pml.py.txt differ