Customize the PML to account for beam path in non vacuum material

tests/test_005_pml_boundaries.py

@@ -21,7 +21,7 @@ def test_reflection_planewave(self):
         # Domain bounds and grid
         xmin, xmax = -1., 1.
         ymin, ymax = -1., 1.
-        zmin, zmax = 0., 10.
+        zmin, zmax = 0., 1.
 
         Nx, Ny = 20, 20
         Nz = 200
@@ -33,10 +33,15 @@ def test_reflection_planewave(self):
         bc_low = ['periodic', 'periodic', 'pec']
         bc_high = ['periodic', 'periodic', 'pml']
 
+        # Test different eps_r and sigma case
+        eps_r = 1.0
+        sigma = 0.0
 
+        # Solver
         solver = wakis.SolverFIT3D(grid, use_stl=False,
-                                bc_low=bc_low, bc_high=bc_high,
-                                n_pml=10)
+                                    bg=[eps_r, 1.0, sigma],
+                                    bc_low=bc_low, bc_high=bc_high,
+                                    n_pml=30)
 
         # Source
         amplitude = 100.
@@ -48,22 +53,97 @@ def test_reflection_planewave(self):
         solver.dt = 1/f/200 #ensure right amplitude
 
         # Simulation
-        Nt = int(2.0*(zmax-zmin-solver.n_pml*solver.dz)/c/solver.dt)
+
+        # Simulation time is extended by a factor eps_r 
+        # to ensure the wave is fully absorbed in the PML
+        Nt = int(eps_r*2.0*(zmax-zmin-solver.n_pml*solver.dz)/c/solver.dt)
+
+        for n in tqdm(range(Nt)):
+            planeWave.update(solver, n*solver.dt)
+            solver.one_step()
+
+            if flag_interactive and n%int(Nt/100) == 0:
+                solver.plot1D('Hy', ylim=(-amplitude, amplitude), pos=[0.5, 0.35, 0.2, 0.1],
+                            off_screen=True, title=f'005_Hy', n=n)
+                solver.plot1D('Ex', ylim=(-amplitude*c*mu_0, amplitude*c*mu_0), pos=[0.5, 0.35, 0.2, 0.1],
+                            off_screen=True, title=f'005_Ex', n=n)                
+
+        reflectionH = solver.H.get_abs()[Nx//2,Ny//2,:].max()
+        reflectionE = solver.E.get_abs()[Nx//2,Ny//2,:].max()/(mu_0*c)
+        assert reflectionH <= 10, f"PML H reflection >10% with eps_r={eps_r}, sigma={sigma}"
+        assert reflectionE <= 10, f"PML E reflection >10% with eps_r={eps_r}, sigma={sigma}"
+
+        if flag_interactive:
+            os.system(f'convert -delay 10 -loop 0 005_Hy*.png 005_Hy_planewave.gif')
+            os.system(f'convert -delay 10 -loop 0 005_Ex*.png 005_Ex_planewave.gif')
+            os.system(f'rm 005_Hy*.png')
+            os.system(f'rm 005_Ex*.png')
+
+            solver.plot2D('Ex', plane='ZX', pos=0.5, cmap='bwr', 
+                interpolation='spline36', n=n, vmin=-amplitude*c*mu_0, vmax=amplitude*c*mu_0,
+                off_screen=True, title=f'005_Ex2d') 
+
+            solver.plot2D('Hy', plane='ZX', pos=0.5, cmap='bwr', 
+                interpolation='spline36', n=n, vmin=-amplitude, vmax=amplitude,
+                off_screen=True, title=f'005_Hy2d') 
+
+
+    def test_reflection_planewave_resistive_material(self):
+        print("\n---------- Initializing simulation ------------------")
+        # Domain bounds and grid
+        xmin, xmax = -1., 1.
+        ymin, ymax = -1., 1.
+        zmin, zmax = 0., 1.
+
+        Nx, Ny = 20, 20
+        Nz = 200
+
+        grid = wakis.GridFIT3D(xmin, xmax, ymin, ymax, zmin, zmax, 
+                            Nx, Ny, Nz)
+
+        # Boundary conditions and solver
+        bc_low = ['periodic', 'periodic', 'pec']
+        bc_high = ['periodic', 'periodic', 'pml']
+
+        # Test different eps_r and sigma case
+        eps_r = 3.0
+        sigma = 1.0e-3
+
+        # Solver
+        solver = wakis.SolverFIT3D(grid, use_stl=False,
+                                    bg=[eps_r, 1.0, sigma],
+                                    bc_low=bc_low, bc_high=bc_high,
+                                    n_pml=30)
+
+        # Source
+        amplitude = 100.
+        nodes = 7
+        f = 15/((solver.z.max()-solver.z.min())/c)
+        planeWave = wakis.sources.PlaneWave(xs=slice(0, Nx), ys=slice(0,Ny), zs=0, 
+                                            beta=1.0, amplitude=amplitude,
+                                            f=f, nodes=nodes, phase=np.pi/2)
+        solver.dt = 1/f/200 #ensure right amplitude
+
+        # Simulation
+
+        # Simulation time is extended by a factor eps_r 
+        # to ensure the wave is fully absorbed in the PML
+        Nt = int(eps_r*2.0*(zmax-zmin-solver.n_pml*solver.dz)/c/solver.dt)
 
         for n in tqdm(range(Nt)):
             planeWave.update(solver, n*solver.dt)
             solver.one_step()
 
             if flag_interactive and n%int(Nt/100) == 0:
                 solver.plot1D('Hy', ylim=(-amplitude, amplitude), pos=[0.5, 0.35, 0.2, 0.1],
-                               off_screen=True, title=f'005_Hy', n=n)
+                            off_screen=True, title=f'005_Hy', n=n)
                 solver.plot1D('Ex', ylim=(-amplitude*c*mu_0, amplitude*c*mu_0), pos=[0.5, 0.35, 0.2, 0.1],
-                               off_screen=True, title=f'005_Ex', n=n)                
+                            off_screen=True, title=f'005_Ex', n=n)                
 
         reflectionH = solver.H.get_abs()[Nx//2,Ny//2,:].max()
         reflectionE = solver.E.get_abs()[Nx//2,Ny//2,:].max()/(mu_0*c)
-        assert reflectionH == pytest.approx(10, 0.1), "PML H reflection >10%"
-        assert reflectionE == pytest.approx(12, 0.2), "PML E reflection >10%"
+        assert reflectionH <= 10, f"PML H reflection >10% with eps_r={eps_r}, sigma={sigma}"
+        assert reflectionE <= 10, f"PML E reflection >10% with eps_r={eps_r}, sigma={sigma}"
 
         if flag_interactive:
             os.system(f'convert -delay 10 -loop 0 005_Hy*.png 005_Hy_planewave.gif')
@@ -78,7 +158,78 @@ def test_reflection_planewave(self):
             solver.plot2D('Hy', plane='ZX', pos=0.5, cmap='bwr', 
                 interpolation='spline36', n=n, vmin=-amplitude, vmax=amplitude,
                 off_screen=True, title=f'005_Hy2d') 
+
+
+    def test_reflection_planewave_high_resistivity_material(self):
+        print("\n---------- Initializing simulation ------------------")
+        # Domain bounds and grid
+        xmin, xmax = -1., 1.
+        ymin, ymax = -1., 1.
+        zmin, zmax = 0., 1.
+
+        Nx, Ny = 20, 20
+        Nz = 200
+
+        grid = wakis.GridFIT3D(xmin, xmax, ymin, ymax, zmin, zmax, 
+                            Nx, Ny, Nz)
+
+        # Boundary conditions and solver
+        bc_low = ['periodic', 'periodic', 'pec']
+        bc_high = ['periodic', 'periodic', 'pml']
+
+        # Test different eps_r and sigma case
+        eps_r = 1.0
+        sigma = 1
+
+        # Solver
+        solver = wakis.SolverFIT3D(grid, use_stl=False,
+                                    bg=[eps_r, 1.0, sigma],
+                                    bc_low=bc_low, bc_high=bc_high,
+                                    n_pml=30)
+
+        # Source
+        amplitude = 100.
+        nodes = 7
+        f = 15/((solver.z.max()-solver.z.min())/c)
+        planeWave = wakis.sources.PlaneWave(xs=slice(0, Nx), ys=slice(0,Ny), zs=0, 
+                                            beta=1.0, amplitude=amplitude,
+                                            f=f, nodes=nodes, phase=np.pi/2)
+        solver.dt = 1/f/200 #ensure right amplitude
+
+        # Simulation
+
+        # Simulation time is extended by a factor eps_r 
+        # to ensure the wave is fully absorbed in the PML
+        Nt = int(eps_r*2.0*(zmax-zmin-solver.n_pml*solver.dz)/c/solver.dt)
+
+        for n in tqdm(range(Nt)):
+            planeWave.update(solver, n*solver.dt)
+            solver.one_step()
+
+            if flag_interactive and n%int(Nt/100) == 0:
+                solver.plot1D('Hy', ylim=(-amplitude, amplitude), pos=[0.5, 0.35, 0.2, 0.1],
+                            off_screen=True, title=f'005_Hy', n=n)
+                solver.plot1D('Ex', ylim=(-amplitude*c*mu_0, amplitude*c*mu_0), pos=[0.5, 0.35, 0.2, 0.1],
+                            off_screen=True, title=f'005_Ex', n=n)                
 
+        reflectionH = solver.H.get_abs()[Nx//2,Ny//2,:].max()
+        reflectionE = solver.E.get_abs()[Nx//2,Ny//2,:].max()/(mu_0*c)
+        assert reflectionH <= 10, f"PML H reflection >10% with eps_r={eps_r}, sigma={sigma}"
+        assert reflectionE <= 10, f"PML E reflection >10% with eps_r={eps_r}, sigma={sigma}"
+
+        if flag_interactive:
+            os.system(f'convert -delay 10 -loop 0 005_Hy*.png 005_Hy_planewave.gif')
+            os.system(f'convert -delay 10 -loop 0 005_Ex*.png 005_Ex_planewave.gif')
+            os.system(f'rm 005_Hy*.png')
+            os.system(f'rm 005_Ex*.png')
+
+            solver.plot2D('Ex', plane='ZX', pos=0.5, cmap='bwr', 
+                interpolation='spline36', n=n, vmin=-amplitude*c*mu_0, vmax=amplitude*c*mu_0,
+                off_screen=True, title=f'005_Ex2d') 
+
+            solver.plot2D('Hy', plane='ZX', pos=0.5, cmap='bwr', 
+                interpolation='spline36', n=n, vmin=-amplitude, vmax=amplitude,
+                off_screen=True, title=f'005_Hy2d') 
 
     def _pml_func():
 

wakis/solverFIT3D.py

@@ -66,6 +66,8 @@ def __init__(self, grid, wake=None, cfln=0.5, dt=None,
             If true, activates all the solids and materials passed to the `grid` object
         use_gpu: bool, default False,
             Using cupyx, enables GPU accelerated computation of every timestep
+        n_pml: int, default 10,
+            Number of PML cells at the boundaries of the simulation box
         bg: list, default [1.0, 1.0]
             Background material for the simulation box [eps_r, mu_r, sigma]. Default is vacuum.
             It supports any material from the material library in `materials.py`, of a
@@ -207,9 +209,10 @@ def __init__(self, grid, wake=None, cfln=0.5, dt=None,
             if verbose:
                 print('Filling PML sigmas...')
             self.n_pml = n_pml
-            self.pml_lo = 5e-3
-            self.pml_hi = 1.e-1
+            self.pml_lo = 5.0e-3
+            self.pml_hi = 10.0
             self.pml_func = np.geomspace
+            self.pml_eps_r = 1.0
             self.fill_pml_sigmas()
             self.update_logger(['n_pml'])
 
@@ -758,59 +761,89 @@ def fill_pml_sigmas(self):
         # Fill
         if self.bc_low[0].lower() == 'pml':
             #sx[0:self.n_pml] = eps_0/(2*self.dt)*((self.x[self.n_pml] - self.x[:self.n_pml])/(self.n_pml*self.dx))**pml_exp
-            sx[0:self.n_pml] = np.linspace( self.pml_hi, self.pml_lo, self.n_pml)
+            sx[0:self.n_pml] = self.pml_func(self.pml_hi, self.pml_lo, self.n_pml)
             for d in ['x', 'y', 'z']:
+                # Get the properties from the layer before the PML
+                # Take the values at the center of the yz plane
+                ieps_0_pml = self.ieps[self.n_pml+1, self.Ny//2, self.Nz//2, d]
+                sigma_0_pml = self.sigma[self.n_pml+1, self.Ny//2, self.Nz//2, d]
+                sigma_mult_pml = 1 if sigma_0_pml < 1 else sigma_0_pml # avoid null sigma in PML for relaxation time computation
                 for i in range(self.n_pml):
-                    self.ieps[i, :, :, d] = 1./eps_0
-                    self.sigma[i, :, :, d] = sx[i]
-                    #if sx[i] > 0 : self.ieps[i, :, :, d] = 1/(eps_0+sx[i]*(2*self.dt))
+                    self.ieps[i, :, :, d] = ieps_0_pml
+                    self.sigma[i, :, :, d] = sigma_0_pml + sigma_mult_pml * sx[i]
+                    #if sx[i] > 0 : self.ieps[i, :, :, d] = 1/(eps_0+sx[i]*(2*self.dt)) 
 
         if self.bc_low[1].lower() == 'pml':
             #sy[0:self.n_pml] = 1/(2*self.dt)*((self.y[self.n_pml] - self.y[:self.n_pml])/(self.n_pml*self.dy))**pml_exp
-            sy[0:self.n_pml] = self.pml_func( self.pml_hi, self.pml_lo, self.n_pml)
+            sy[0:self.n_pml] = self.pml_func(self.pml_hi, self.pml_lo, self.n_pml)
             for d in ['x', 'y', 'z']:
+                # Get the properties from the layer before the PML
+                # Take the values at the center of the xz plane
+                ieps_0_pml = self.ieps[self.Nx//2, self.n_pml+1, self.Nz//2, d]
+                sigma_0_pml = self.sigma[self.Nx//2, self.n_pml+1, self.Nz//2, d]
+                sigma_mult_pml = 1 if sigma_0_pml < 1 else sigma_0_pml # avoid null sigma in PML for relaxation time computation
                 for j in range(self.n_pml):
-                    self.ieps[:, j, :, d] = 1./eps_0
-                    self.sigma[:, j, :, d] = sy[j]
-                    #if sy[j] > 0 : self.ieps[:, j, :, d] = 1/(eps_0+sy[j]*(2*self.dt))
+                    self.ieps[:, j, :, d] = ieps_0_pml
+                    self.sigma[:, j, :, d] = sigma_0_pml + sigma_mult_pml * sy[j]
+                    #if sy[j] > 0 : self.ieps[:, j, :, d] = 1/(eps_0+sy[j]*(2*self.dt)) 
 
         if self.bc_low[2].lower() == 'pml':
             #sz[0:self.n_pml] = eps_0/(2*self.dt)*((self.z[self.n_pml] - self.z[:self.n_pml])/(self.n_pml*self.dz))**pml_exp
-            sz[0:self.n_pml] = self.pml_func( self.pml_hi, self.pml_lo, self.n_pml)
+            sz[0:self.n_pml] = self.pml_func(self.pml_hi, self.pml_lo, self.n_pml)
             for d in ['x', 'y', 'z']:
+                # Get the properties from the layer before the PML
+                # Take the values at the center of the xy plane
+                ieps_0_pml = self.ieps[self.Nx//2, self.Ny//2, self.n_pml+1, d]
+                sigma_0_pml = self.sigma[self.Nx//2, self.Ny//2, self.n_pml+1, d]
+                sigma_mult_pml = 1 if sigma_0_pml < 1 else sigma_0_pml # avoid null sigma in PML for relaxation time computation
                 for k in range(self.n_pml):
-                    self.ieps[:, :, k, d] = 1./eps_0
-                    self.sigma[:, :, k, d] = sz[k]
-                    #if sz[k] > 0. : self.ieps[:, :, k, d] = 1/(np.mean(sz[:self.n_pml])*eps_0)
+                    self.ieps[:, :, k, d] = ieps_0_pml
+                    self.sigma[:, :, k, d] = sigma_0_pml + sigma_mult_pml * sz[k]
+                    #if sz[k] > 0. : self.ieps[:, :, k, d] = 1/(np.mean(sz[:self.n_pml])*eps_0) 
 
         if self.bc_high[0].lower() == 'pml':
             #sx[-self.n_pml:] = 1/(2*self.dt)*((self.x[-self.n_pml:] - self.x[-self.n_pml])/(self.n_pml*self.dx))**pml_exp
-            sx[-self.n_pml:] = self.pml_func( self.pml_lo, self.pml_hi, self.n_pml)
+            sx[-self.n_pml:] = self.pml_func(self.pml_lo, self.pml_hi, self.n_pml)
             for d in ['x', 'y', 'z']:
+                # Get the properties from the layer before the PML
+                # Take the values at the center of the yz plane
+                ieps_0_pml = self.ieps[-(self.n_pml+1), self.Ny//2, self.Nz//2, d]
+                sigma_0_pml = self.sigma[-(self.n_pml+1), self.Ny//2, self.Nz//2, d]
+                sigma_mult_pml = 1 if sigma_0_pml < 1 else sigma_0_pml # avoid null sigma in PML for relaxation time computation
                 for i in range(self.n_pml):
                     i +=1
-                    self.ieps[-i, :, :, d] = 1./eps_0
-                    self.sigma[-i, :, :, d] = sx[-i]
-                    #if sx[-i] > 0 : self.ieps[-i, :, :, d] = 1/(eps_0+sx[-i]*(2*self.dt))
+                    self.ieps[-i, :, :, d] = ieps_0_pml
+                    self.sigma[-i, :, :, d] = sigma_0_pml + sigma_mult_pml * sx[-i]
+                    #if sx[-i] > 0 : self.ieps[-i, :, :, d] = 1/(eps_0+sx[-i]*(2*self.dt)) 
 
         if self.bc_high[1].lower() == 'pml':
             #sy[-self.n_pml:] = 1/(2*self.dt)*((self.y[-self.n_pml:] - self.y[-self.n_pml])/(self.n_pml*self.dy))**pml_exp
-            sy[-self.n_pml:] = self.pml_func( self.pml_lo, self.pml_hi, self.n_pml)
+            sy[-self.n_pml:] = self.pml_func(self.pml_lo, self.pml_hi, self.n_pml)
             for d in ['x', 'y', 'z']:
+                # Get the properties from the layer before the PML
+                # Take the values at the center of the xz plane
+                ieps_0_pml = self.ieps[self.Nx//2, -(self.n_pml+1), self.Nz//2, d]
+                sigma_0_pml = self.sigma[self.Nx//2, -(self.n_pml+1), self.Nz//2, d]
+                sigma_mult_pml = 1 if sigma_0_pml < 1 else sigma_0_pml # avoid null sigma in PML for relaxation time computation
                 for j in range(self.n_pml):
                     j +=1
-                    self.ieps[:, -j, :, d] = 1./eps_0
-                    self.sigma[:, -j, :, d] = sy[-j]
-                    #if sy[-j] > 0 : self.ieps[:, -j, :, d] = 1/(eps_0+sy[-j]*(2*self.dt))
+                    self.ieps[:, -j, :, d] = ieps_0_pml
+                    self.sigma[:, -j, :, d] = sigma_0_pml + sigma_mult_pml * sy[-j]
+                    #if sy[-j] > 0 : self.ieps[:, -j, :, d] = 1/(eps_0+sy[-j]*(2*self.dt)) 
 
         if self.bc_high[2].lower() == 'pml':
             #sz[-self.n_pml:] = eps_0/(2*self.dt)*((self.z[-self.n_pml:] - self.z[-self.n_pml])/(self.n_pml*self.dz))**pml_exp
-            sz[-self.n_pml:] = self.pml_func( self.pml_lo, self.pml_hi, self.n_pml)
+            sz[-self.n_pml:] = self.pml_func(self.pml_lo, self.pml_hi, self.n_pml)
             for d in ['x', 'y', 'z']:
+                # Get the properties from the layer before the PML
+                # Take the values at the center of the xy plane
+                ieps_0_pml = self.ieps[self.Nx//2, self.Ny//2, -(self.n_pml+1), d]
+                sigma_0_pml = self.sigma[self.Nx//2, self.Ny//2, -(self.n_pml+1), d]
+                sigma_mult_pml = 1 if sigma_0_pml < 1 else sigma_0_pml # avoid null sigma in PML for relaxation time computation
                 for k in range(self.n_pml):
                     k +=1
-                    self.ieps[:, :, -k, d] = 1./eps_0
-                    self.sigma[:, :, -k, d] = sz[-k]
+                    self.ieps[:, :, -k, d] = ieps_0_pml
+                    self.sigma[:, :, -k, d] = sigma_0_pml + sigma_mult_pml * sz[-k]
                     #self.ieps[:, :, -k, d] = 1/(np.mean(sz[-self.n_pml:])*eps_0)
 
     def get_abc(self):