diff options
Diffstat (limited to 'Examples/Modules')
| -rwxr-xr-x | Examples/Modules/dive_cleaning/analysis.py | 107 | ||||
| -rw-r--r-- | Examples/Modules/dive_cleaning/inputs_3d | 34 |
2 files changed, 141 insertions, 0 deletions
diff --git a/Examples/Modules/dive_cleaning/analysis.py b/Examples/Modules/dive_cleaning/analysis.py new file mode 100755 index 000000000..5353a0f6f --- /dev/null +++ b/Examples/Modules/dive_cleaning/analysis.py @@ -0,0 +1,107 @@ +#!/usr/bin/env python + +# Copyright 2019-2020 Axel Huebl, Remi Lehe +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + +""" +This script checks the dive cleaning and PML routines, by initializing a +Gaussian beam in the center of the simulation box, and let the error in +space-charge field propagate away and be absorbed in the PML. + +This script verifies that the field at the end of the simulation corresponds +to the theoretical field of a Gaussian beam. +""" +import sys +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +import yt +import numpy as np +import scipy.constants as scc +from scipy.special import gammainc +yt.funcs.mylog.setLevel(0) + +# Parameters from the Simulation +Qtot = -1.e-20 +r0 = 2.e-6 + +# Open data file +filename = sys.argv[1] +ds = yt.load( filename ) +# Extract data +ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) +Ex_array = ad0['Ex'].to_ndarray().squeeze() +if ds.dimensionality == 2: + # Rename the z dimension as y, so as to make this script work for 2d and 3d + Ey_array = ad0['Ez'].to_ndarray().squeeze() + E_array = ( Ex_array**2 + Ey_array**2 )**.5 + relative_tolerance = 0.1 +elif ds.dimensionality == 3: + Ey_array = ad0['Ey'].to_ndarray() + Ez_array = ad0['Ez'].to_ndarray() + E_array = ( Ex_array**2 + Ey_array**2 + Ez_array**2 )**.5 + relative_tolerance = 0.15 + +# Extract grid coordinates +Nx, Ny, Nz = ds.domain_dimensions +xmin, ymin, zmin = ds.domain_left_edge.v +Lx, Ly, Lz = ds.domain_width.v +x = xmin + Lx/Nx*(0.5+np.arange(Nx)) +y = ymin + Ly/Ny*(0.5+np.arange(Ny)) +z = zmin + Lz/Nz*(0.5+np.arange(Nz)) + +# Compute theoretical field +if ds.dimensionality == 2: + x_2d, y_2d = np.meshgrid(x, y, indexing='ij') + r2 = x_2d**2 + y_2d**2 + factor = (Qtot/r0)/(2*np.pi*scc.epsilon_0*r2) * (1-np.exp(-r2/(2*r0**2))) + Ex_th = x_2d * factor + Ey_th = y_2d * factor + E_th = ( Ex_th**2 + Ey_th**2 )**.5 +elif ds.dimensionality == 3: + x_2d, y_2d, z_2d = np.meshgrid(x, y, z, indexing='ij') + r2 = x_2d**2 + y_2d**2 + z_2d**2 + factor = Qtot/(4*np.pi*scc.epsilon_0*r2**1.5) * gammainc(3./2, r2/(2.*r0**2)) + Ex_th = factor*x_2d + Ey_th = factor*y_2d + Ez_th = factor*z_2d + E_th = ( Ex_th**2 + Ey_th**2 + Ez_th**2 )**.5 + +# Plot theory and data +def make_2d(arr): + if arr.ndim == 3: + return arr[:,:,Nz//2] + else: + return arr +plt.figure(figsize=(10,10)) +plt.subplot(221) +plt.title('E: Theory') +plt.imshow(make_2d(E_th)) +plt.colorbar() +plt.subplot(222) +plt.title('E: Simulation') +plt.imshow(make_2d(E_array)) +plt.colorbar() +plt.subplot(223) +plt.title('E: Diff') +plt.imshow(make_2d(E_th-E_array)) +plt.colorbar() +plt.subplot(224) +plt.title('E: Relative diff') +plt.imshow(make_2d((E_th-E_array)/E_th)) +plt.colorbar() +plt.savefig('Comparison.png') + +# Automatically check the results +def check(E, E_th, label): + print( 'Relative error in %s: %.3f'%( + label, abs(E-E_th).max()/E_th.max())) + assert np.allclose( E, E_th, atol=relative_tolerance*E_th.max() ) + +check( Ex_array, Ex_th, 'Ex' ) +check( Ey_array, Ey_th, 'Ey' ) +if ds.dimensionality == 3: + check( Ez_array, Ez_th, 'Ez' ) diff --git a/Examples/Modules/dive_cleaning/inputs_3d b/Examples/Modules/dive_cleaning/inputs_3d new file mode 100644 index 000000000..97e7b4331 --- /dev/null +++ b/Examples/Modules/dive_cleaning/inputs_3d @@ -0,0 +1,34 @@ +max_step = 128 +amr.n_cell = 64 64 64 +amr.max_grid_size = 32 +amr.max_level = 0 + +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 0 0 0 # Is periodic? +geometry.prob_lo = -50.e-6 -50.e-6 -50.e-6 +geometry.prob_hi = 50.e-6 50.e-6 50.e-6 + +warpx.do_dive_cleaning = 1 +warpx.do_pml = 1 + +particles.nspecies = 1 +particles.species_names = beam +beam.charge = -q_e +beam.mass = 1.e30 +beam.injection_style = "gaussian_beam" +beam.x_rms = 2.e-6 +beam.y_rms = 2.e-6 +beam.z_rms = 2.e-6 +beam.x_m = 0. +beam.y_m = 0. +beam.z_m = 0.e-6 +beam.npart = 20000 +beam.q_tot = -1.e-20 +beam.momentum_distribution_type = "gaussian" +beam.ux_m = 0.0 +beam.uy_m = 0.0 +beam.uz_m = 0.0 + +diagnostics.diags_names = diag1 +diag1.period = 8 +diag1.diag_type = Full |