diff options
Diffstat (limited to 'Examples')
3 files changed, 389 insertions, 0 deletions
diff --git a/Examples/Tests/Langmuir/inputs.multi.rzmodes.rt b/Examples/Tests/Langmuir/inputs.multi.rzmodes.rt new file mode 100644 index 000000000..9808163f4 --- /dev/null +++ b/Examples/Tests/Langmuir/inputs.multi.rzmodes.rt @@ -0,0 +1,95 @@ +# Runs test with RZ multimode solver +# Maximum number of time steps +max_step = 40 + +# number of grid points +amr.n_cell = 64 200 + +# Maximum allowable size of each subdomain in the problem domain; +# this is used to decompose the domain for parallel calculations. +amr.max_grid_size = 32 + +# Maximum level in hierarchy (for now must be 0, i.e., one level in total) +amr.max_level = 0 + +amr.plot_int = 40 # How often to write plotfiles. "<= 0" means no plotfiles. + +# Geometry +geometry.coord_sys = 1 # 0: Cartesian +geometry.is_periodic = 0 1 # Is periodic? +geometry.prob_lo = 0.e-6 -20.e-6 # physical domain +geometry.prob_hi = 20.e-6 20.e-6 + +warpx.nmodes = 3 + +warpx.serialize_ics = 1 +warpx.plot_raw_fields = 1 + +# Verbosity +warpx.verbose = 1 + +# Interpolation +interpolation.nox = 1 +interpolation.noy = 1 +interpolation.noz = 1 + +# CFL +warpx.cfl = 1.0 + +# dive is not supported with RZ multimode +warpx.do_dive_cleaning = 0 + +# Parameters for the plasma wave +my_constants.epsilon0 = 0.001 +my_constants.epsilon1 = 0.001 +my_constants.epsilon2 = 0.001 +my_constants.kp = 266125.0928256017 +my_constants.k0 = 314159.2653589793 +my_constants.w0 = 5.e-6 +# Note: kp is calculated in SI for a density of 2e24 +# k0 is calculated so as to have 2 periods within the 40e-6 wide box. + +# Particles +particles.nspecies = 2 +particles.species_names = electrons ions + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NUniformPerCell" +electrons.num_particles_per_cell_each_dim = 2 8 2 +electrons.xmin = 0.e-6 +electrons.xmax = 18.e-6 +electrons.zmin = -20.e-6 +electrons.zmax = +20.e-6 + +electrons.profile = constant +electrons.density = 2.e24 # number of electrons per m^3 +electrons.momentum_distribution_type = parse_momentum_function +electrons.momentum_function_ux(x,y,z) = "+ epsilon0/kp*2*x/w0**2*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + - epsilon1/kp*2/w0*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + + epsilon1/kp*4*x**2/w0**3*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + - epsilon2/kp*8*x/w0**2*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + + epsilon2/kp*8*x*(x**2-y**2)/w0**4*exp(-(x**2+y**2)/w0**2)*sin(k0*z)" +electrons.momentum_function_uy(x,y,z) = "+ epsilon0/kp*2*y/w0**2*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + + epsilon1/kp*4*x*y/w0**3*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + + epsilon2/kp*8*y/w0**2*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + + epsilon2/kp*8*y*(x**2-y**2)/w0**4*exp(-(x**2+y**2)/w0**2)*sin(k0*z)" +electrons.momentum_function_uz(x,y,z) = "- epsilon0/kp*k0*exp(-(x**2+y**2)/w0**2)*cos(k0*z) + - epsilon1/kp*k0*2*x/w0*exp(-(x**2+y**2)/w0**2)*cos(k0*z) + - epsilon2/kp*k0*4*(x**2-y**2)/w0**2*exp(-(x**2+y**2)/w0**2)*cos(k0*z)" + +ions.charge = q_e +ions.mass = m_p +ions.injection_style = "NUniformPerCell" +ions.num_particles_per_cell_each_dim = 2 8 2 +ions.xmin = 0.e-6 +ions.xmax = 18.e-6 +ions.zmin = -20.e-6 +ions.zmax = +20.e-6 + +ions.profile = constant +ions.density = 2.e24 # number of ions per m^3 +ions.momentum_distribution_type = constant +ions.ux = 0. +ions.uy = 0. +ions.uz = 0. diff --git a/Examples/Tests/Langmuir/langmuir_PICMI_rz_multimode_analyze.py b/Examples/Tests/Langmuir/langmuir_PICMI_rz_multimode_analyze.py new file mode 100644 index 000000000..9d9fa40f6 --- /dev/null +++ b/Examples/Tests/Langmuir/langmuir_PICMI_rz_multimode_analyze.py @@ -0,0 +1,170 @@ +# This is a script that analyses the multimode simulation results. +# This simulates a RZ multimode periodic plasma wave. +# The electric field from the simulation is compared to the analytic value + +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from pywarpx import picmi + +nr = 64 +nz = 200 + +rmin = 0.e0 +zmin = 0.e0 +rmax = +20.e-6 +zmax = +40.e-6 + +# Parameters describing particle distribution +density = 2.e24 +epsilon0 = 0.001 +epsilon1 = 0.001 +epsilon2 = 0.001 +w0 = 5.e-6 +n_osc_z = 3 + +# Wave vector of the wave +k0 = 2.*np.pi*n_osc_z/(zmax - zmin) + +# Plasma frequency +wp = np.sqrt((density*picmi.q_e**2)/(picmi.m_e*picmi.ep0)) +kp = wp/picmi.c + +uniform_plasma = picmi.UniformDistribution(density = density, + upper_bound = [+18e-6, +18e-6, None], + directed_velocity = [0., 0., 0.]) + +momentum_expressions = ["""+ epsilon0/kp*2*x/w0**2*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + - epsilon1/kp*2/w0*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + + epsilon1/kp*4*x**2/w0**3*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + - epsilon2/kp*8*x/w0**2*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + + epsilon2/kp*8*x*(x**2-y**2)/w0**4*exp(-(x**2+y**2)/w0**2)*sin(k0*z)""", + """+ epsilon0/kp*2*y/w0**2*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + + epsilon1/kp*4*x*y/w0**3*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + + epsilon2/kp*8*y/w0**2*exp(-(x**2+y**2)/w0**2)*sin(k0*z) + + epsilon2/kp*8*y*(x**2-y**2)/w0**4*exp(-(x**2+y**2)/w0**2)*sin(k0*z)""", + """- epsilon0/kp*k0*exp(-(x**2+y**2)/w0**2)*cos(k0*z) + - epsilon1/kp*k0*2*x/w0*exp(-(x**2+y**2)/w0**2)*cos(k0*z) + - epsilon2/kp*k0*4*(x**2-y**2)/w0**2*exp(-(x**2+y**2)/w0**2)*cos(k0*z)"""] + +analytic_plasma = picmi.AnalyticDistribution(density_expression = density, + upper_bound = [+18e-6, +18e-6, None], + epsilon0 = epsilon0, + epsilon1 = epsilon1, + epsilon2 = epsilon2, + kp = kp, + k0 = k0, + w0 = w0, + momentum_expressions = momentum_expressions) + +electrons = picmi.Species(particle_type='electron', name='electrons', initial_distribution=analytic_plasma) +ions = picmi.Species(particle_type='proton', name='protons', initial_distribution=uniform_plasma) + +grid = picmi.CylindricalGrid(number_of_cells = [nr, nz], + n_azimuthal_modes = 3, + lower_bound = [rmin, zmin], + upper_bound = [rmax, zmax], + lower_boundary_conditions = ['dirichlet', 'periodic'], + upper_boundary_conditions = ['dirichlet', 'periodic'], + moving_window_velocity = [0., 0.], + warpx_max_grid_size=64) + +solver = picmi.ElectromagneticSolver(grid=grid, cfl=1.) + +sim = picmi.Simulation(solver = solver, + max_steps = 40, + verbose = 1, + warpx_plot_int = 40, + warpx_current_deposition_algo = 'direct', + warpx_field_gathering_algo = 'standard', + warpx_particle_pusher_algo = 'boris') + +sim.add_species(electrons, layout=picmi.GriddedLayout(n_macroparticle_per_cell=[2,8,2], grid=grid)) +sim.add_species(ions, layout=picmi.GriddedLayout(n_macroparticle_per_cell=[2,8,2], grid=grid)) + +# write_inputs will create an inputs file that can be used to run +# with the compiled version. +#sim.write_input_file(file_name='inputsrz_from_PICMI') + +# Alternatively, sim.step will run WarpX, controlling it from Python +sim.step() + + +# Below is WarpX specific code to check the results. +import pywarpx +from pywarpx.fields import * + +def calcEr( z, r, k0, w0, wp, t, epsilons) : + """ + Return the radial electric field as an array + of the same length as z and r, in the half-plane theta=0 + """ + Er_array = ( + epsilons[0] * picmi.m_e*picmi.c**2/picmi.q_e * 2*r/w0**2 * + np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.sin( wp*t ) + - epsilons[1] * picmi.m_e*picmi.c**2/picmi.q_e * 2/w0 * + np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.sin( wp*t ) + + epsilons[1] * picmi.m_e*picmi.c**2/picmi.q_e * 4*r**2/w0**3 * + np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.sin( wp*t ) + - epsilons[2] * picmi.m_e*picmi.c**2/picmi.q_e * 8*r/w0**2 * + np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.sin( wp*t ) + + epsilons[2] * picmi.m_e*picmi.c**2/picmi.q_e * 8*r**3/w0**4 * + np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.sin( wp*t )) + return( Er_array ) + +def calcEz( z, r, k0, w0, wp, t, epsilons) : + """ + Return the longitudinal electric field as an array + of the same length as z and r, in the half-plane theta=0 + """ + Ez_array = ( + - epsilons[0] * picmi.m_e*picmi.c**2/picmi.q_e * k0 * + np.exp( -r**2/w0**2 ) * np.cos( k0*z ) * np.sin( wp*t ) + - epsilons[1] * picmi.m_e*picmi.c**2/picmi.q_e * k0 * 2*r/w0 * + np.exp( -r**2/w0**2 ) * np.cos( k0*z ) * np.sin( wp*t ) + - epsilons[2] * picmi.m_e*picmi.c**2/picmi.q_e * k0 * 4*r**2/w0**2 * + np.exp( -r**2/w0**2 ) * np.cos( k0*z ) * np.sin( wp*t )) + return( Ez_array ) + +# Get node centered coordinates +dr = (rmax - rmin)/nr +dz = (zmax - zmin)/nz +coords = np.indices([nr+1, nz+1], 'd') +rr = rmin + coords[0]*dr +zz = zmin + coords[1]*dz + +# Current time of the simulation +t0 = pywarpx._libwarpx.libwarpx.warpx_gett_new(0) + +# Get the raw field data. Note that these are the real and imaginary +# parts of the fields for each azimuthal mode. +Ex_sim_modes = ExWrapper()[...] +Ez_sim_modes = EzWrapper()[...] + +# Sum the real components to get the field along x-axis (theta = 0) +Er_sim = np.sum(Ex_sim_modes[:,:,::2], axis=2) +Ez_sim = np.sum(Ez_sim_modes[:,:,::2], axis=2) + +# The analytical solutions +Er_th = calcEr(zz[:-1,:], rr[:-1,:] + dr/2., k0, w0, wp, t0, [epsilon0, epsilon1, epsilon2]) +Ez_th = calcEz(zz[:,:-1] + dz/2., rr[:,:-1], k0, w0, wp, t0, [epsilon0, epsilon1, epsilon2]) + +max_error_Er = abs(Er_sim - Er_th).max()/abs(Er_th).max() +max_error_Ez = abs(Ez_sim - Ez_th).max()/abs(Ez_th).max() +print("Max error Er %e"%max_error_Er) +print("Max error Ez %e"%max_error_Er) + +# Plot the last field from the loop (Ez at iteration 40) +plt.subplot2grid( (1,2), (0,0) ) +plt.imshow( Ez_sim ) +plt.colorbar() +plt.title('Ez, last iteration\n(simulation)') +plt.subplot2grid( (1,2), (0,1) ) +plt.imshow( Ez_th ) +plt.colorbar() +plt.title('Ez, last iteration\n(theory)') +plt.tight_layout() +plt.savefig('langmuir_multi_rz_multimode_analysis.png') + +assert max(max_error_Er, max_error_Ez) < 0.02 diff --git a/Examples/Tests/Langmuir/langmuir_multi_rz_multimode_analysis.py b/Examples/Tests/Langmuir/langmuir_multi_rz_multimode_analysis.py new file mode 100644 index 000000000..b87dd0fae --- /dev/null +++ b/Examples/Tests/Langmuir/langmuir_multi_rz_multimode_analysis.py @@ -0,0 +1,124 @@ +#! /usr/bin/env python + +# This is a script that analyses the simulation results from +# the script `inputs.multi.rzmodes.rt`. This simulates a RZ multimode periodic plasma wave. +# The electric field from the simulation is compared to the analytic value +import sys +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +import yt +yt.funcs.mylog.setLevel(50) +import numpy as np +from scipy.constants import e, m_e, epsilon_0, c + +# this will be the name of the plot file +fn = sys.argv[1] + +# Parameters (these parameters must match the parameters in `inputs.multi.rzmm.rt`) +epsilon0 = 0.001 +epsilon1 = 0.001 +epsilon2 = 0.001 +n = 2.e24 +w0 = 5.e-6 +n_osc_z = 2 +gscale = 1 +rmin = 0e-6; rmax = 20.e-6; Nr = 64//gscale +zmin = -20e-6; zmax = 20.e-6; Nz = 200//gscale + +# Wave vector of the wave +k0 = 2.*np.pi*n_osc_z/(zmax-zmin) +# Plasma frequency +wp = np.sqrt((n*e**2)/(m_e*epsilon_0)) +kp = wp/c + +def Jx( z, r, k0, w0, wp, t) : + x = r + y = 0 + return -n*e*c*(+ epsilon0 /kp * 2*x/w0**2 * np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.cos( wp*t ) + - epsilon1 /kp * 2/w0 * np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.cos( wp*t ) + + epsilon1 /kp * 4*x**2/w0**3 * np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.cos( wp*t ) + - epsilon2 /kp * 8*x/w0**2 * np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.cos( wp*t ) + + epsilon2 /kp * 8*x*(x**2-y**2)/w0**4 * np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.cos( wp*t )) + +def Jz( z, r, k0, w0, wp, t) : + x = r + y = 0 + return +n*e*c*(- epsilon0 /kp * k0 * np.exp( -r**2/w0**2 ) * np.cos( k0*z ) * np.cos( wp*t ) + - epsilon1 /kp * k0 * 2*x/w0 * np.exp( -r**2/w0**2 ) * np.cos( k0*z ) * np.cos( wp*t ) + - epsilon2 /kp * k0 * 4*(x**2-y**2)/w0**2 * np.exp( -r**2/w0**2 ) * np.cos( k0*z ) * np.cos( wp*t )) + +def Er( z, r, k0, w0, wp, t) : + """ + Return the radial electric field as an array + of the same length as z and r, in the half-plane theta=0 + """ + Er_array = ( + epsilon0 * m_e*c**2/e * 2*r/w0**2 * + np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.sin( wp*t ) + - epsilon1 * m_e*c**2/e * 2/w0 * + np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.sin( wp*t ) + + epsilon1 * m_e*c**2/e * 4*r**2/w0**3 * + np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.sin( wp*t ) + - epsilon2 * m_e*c**2/e * 8*r/w0**2 * + np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.sin( wp*t ) + + epsilon2 * m_e*c**2/e * 8*r**3/w0**4 * + np.exp( -r**2/w0**2 ) * np.sin( k0*z ) * np.sin( wp*t )) + return( Er_array ) + +def Ez( z, r, k0, w0, wp, t) : + """ + Return the longitudinal electric field as an array + of the same length as z and r, in the half-plane theta=0 + """ + Ez_array = ( + - epsilon0 * m_e*c**2/e * k0 * + np.exp( -r**2/w0**2 ) * np.cos( k0*z ) * np.sin( wp*t ) + - epsilon1 * m_e*c**2/e * k0 * 2*r/w0 * + np.exp( -r**2/w0**2 ) * np.cos( k0*z ) * np.sin( wp*t ) + - epsilon2 * m_e*c**2/e * k0 * 4*r**2/w0**2 * + np.exp( -r**2/w0**2 ) * np.cos( k0*z ) * np.sin( wp*t )) + return( Ez_array ) + +# Read the file +ds = yt.load(fn) +t0 = ds.current_time.to_ndarray().mean() +data = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, + dims=ds.domain_dimensions) + +# Get cell centered coordinates +dr = (rmax - rmin)/Nr +dz = (zmax - zmin)/Nz +coords = np.indices([Nr, Nz], 'd') +rr = rmin + (coords[0] + 0.5)*dr +zz = zmin + (coords[1] + 0.5)*dz + +# Check the validity of the fields +overall_max_error = 0 +Er_sim = data['Ex'].to_ndarray()[:,:,0] +Er_th = Er(zz, rr, k0, w0, wp, t0) +max_error = abs(Er_sim-Er_th).max()/abs(Er_th).max() +print('Er: Max error: %.2e' %(max_error)) +overall_max_error = max( overall_max_error, max_error ) + +Ez_sim = data['Ez'].to_ndarray()[:,:,0] +Ez_th = Ez(zz, rr, k0, w0, wp, t0) +max_error = abs(Ez_sim-Ez_th).max()/abs(Ez_th).max() +print('Ez: Max error: %.2e' %(max_error)) +overall_max_error = max( overall_max_error, max_error ) + +# Plot the last field from the loop (Ez at iteration 40) +plt.subplot2grid( (1,2), (0,0) ) +plt.imshow( Ez_sim ) +plt.colorbar() +plt.title('Ez, last iteration\n(simulation)') +plt.subplot2grid( (1,2), (0,1) ) +plt.imshow( Ez_th ) +plt.colorbar() +plt.title('Ez, last iteration\n(theory)') +plt.tight_layout() +plt.savefig('langmuir_multi_rz_analysis.png') + +# Automatically check the validity +assert overall_max_error < 0.02 + |