diff --git a/tests/test_005_pml_boundaries.py b/tests/test_005_pml_boundaries.py index 6a1bd13..4561651 100644 --- a/tests/test_005_pml_boundaries.py +++ b/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,7 +53,82 @@ 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) @@ -56,14 +136,14 @@ def test_reflection_planewave(self): 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(): diff --git a/wakis/solverFIT3D.py b/wakis/solverFIT3D.py index 7ad4cda..04830e2 100644 --- a/wakis/solverFIT3D.py +++ b/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):