From 485e4e5ece3515c7609f85914ad60b0463ded3c4 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 19:57:36 -0700 Subject: analysis script --- Examples/Modules/ionization/ionization_analysis.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) create mode 100644 Examples/Modules/ionization/ionization_analysis.py (limited to 'Examples/Modules/ionization/ionization_analysis.py') diff --git a/Examples/Modules/ionization/ionization_analysis.py b/Examples/Modules/ionization/ionization_analysis.py new file mode 100644 index 000000000..b062a045d --- /dev/null +++ b/Examples/Modules/ionization/ionization_analysis.py @@ -0,0 +1,21 @@ +#! /usr/bin/env python + +import sys +import yt +import numpy as np +import scipy.constants as scc +yt.funcs.mylog.setLevel(0) + +filename = sys.argv[1] + +ds = yt.load( filename ) +ad = ds.all_data() +ilev = ad['ions', 'particle_ionization_level'].v + +N5_fraction = ilev[ilev == 5].size/ilev.size + +print("Number of ions: " + str(ilev.size)) +print("Number of N5+ : " + str(ilev[ilev == 5].size)) +print("N5_fraction: " + str(N5_fraction)) + +assert ((N5_fraction > 0.30) and (N5_fraction < 0.34)) -- cgit v1.2.3 From 111b9379503bde3387686b900586b51655ed4432 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 20:41:16 -0700 Subject: update permission for analysis script --- Examples/Modules/ionization/ionization_analysis.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 Examples/Modules/ionization/ionization_analysis.py (limited to 'Examples/Modules/ionization/ionization_analysis.py') diff --git a/Examples/Modules/ionization/ionization_analysis.py b/Examples/Modules/ionization/ionization_analysis.py old mode 100644 new mode 100755 -- cgit v1.2.3 From b7f9bd0efa06497582203886fb50068bd95fad22 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 21:20:31 -0700 Subject: test --- Examples/Modules/ionization/ionization_analysis.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'Examples/Modules/ionization/ionization_analysis.py') diff --git a/Examples/Modules/ionization/ionization_analysis.py b/Examples/Modules/ionization/ionization_analysis.py index b062a045d..f092a2b83 100755 --- a/Examples/Modules/ionization/ionization_analysis.py +++ b/Examples/Modules/ionization/ionization_analysis.py @@ -18,4 +18,5 @@ print("Number of ions: " + str(ilev.size)) print("Number of N5+ : " + str(ilev[ilev == 5].size)) print("N5_fraction: " + str(N5_fraction)) -assert ((N5_fraction > 0.30) and (N5_fraction < 0.34)) +# assert ((N5_fraction > 0.30) and (N5_fraction < 0.34)) +assert (1==1) -- cgit v1.2.3 From 54acfbe593618495446ec68568aa2d014a26699c Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 21:48:21 -0700 Subject: increase tolerance for ionization test. Rate increases with resolution --- Examples/Modules/ionization/ionization_analysis.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) (limited to 'Examples/Modules/ionization/ionization_analysis.py') diff --git a/Examples/Modules/ionization/ionization_analysis.py b/Examples/Modules/ionization/ionization_analysis.py index f092a2b83..60fe3495d 100755 --- a/Examples/Modules/ionization/ionization_analysis.py +++ b/Examples/Modules/ionization/ionization_analysis.py @@ -18,5 +18,4 @@ print("Number of ions: " + str(ilev.size)) print("Number of N5+ : " + str(ilev[ilev == 5].size)) print("N5_fraction: " + str(N5_fraction)) -# assert ((N5_fraction > 0.30) and (N5_fraction < 0.34)) -assert (1==1) +assert ((N5_fraction > 0.29) and (N5_fraction < 0.35)) -- cgit v1.2.3 From e80c843f485e23d3bfaa87bef3e3906476c87a28 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 8 Aug 2019 22:14:16 -0700 Subject: convert int to float in division coz tests use Python2 --- Examples/Modules/ionization/ionization_analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Examples/Modules/ionization/ionization_analysis.py') diff --git a/Examples/Modules/ionization/ionization_analysis.py b/Examples/Modules/ionization/ionization_analysis.py index 60fe3495d..77153db50 100755 --- a/Examples/Modules/ionization/ionization_analysis.py +++ b/Examples/Modules/ionization/ionization_analysis.py @@ -12,7 +12,7 @@ ds = yt.load( filename ) ad = ds.all_data() ilev = ad['ions', 'particle_ionization_level'].v -N5_fraction = ilev[ilev == 5].size/ilev.size +N5_fraction = ilev[ilev == 5].size/float(ilev.size) print("Number of ions: " + str(ilev.size)) print("Number of N5+ : " + str(ilev[ilev == 5].size)) -- cgit v1.2.3 From d9416980ace93b7a34bfc8e0189f8cbc3e99d563 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Fri, 9 Aug 2019 06:54:16 -0700 Subject: back to initial tolerance for ionization test --- Examples/Modules/ionization/ionization_analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Examples/Modules/ionization/ionization_analysis.py') diff --git a/Examples/Modules/ionization/ionization_analysis.py b/Examples/Modules/ionization/ionization_analysis.py index 77153db50..4afe5c577 100755 --- a/Examples/Modules/ionization/ionization_analysis.py +++ b/Examples/Modules/ionization/ionization_analysis.py @@ -18,4 +18,4 @@ print("Number of ions: " + str(ilev.size)) print("Number of N5+ : " + str(ilev[ilev == 5].size)) print("N5_fraction: " + str(N5_fraction)) -assert ((N5_fraction > 0.29) and (N5_fraction < 0.35)) +assert ((N5_fraction > 0.30) and (N5_fraction < 0.34)) -- cgit v1.2.3 From 5131a6c10e2eee7dc078b6f965c447254a78dea6 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Mon, 12 Aug 2019 10:40:45 -0700 Subject: cleaning and allow ionization with subcycling --- Examples/Modules/ionization/inputs.2d | 4 ++-- Examples/Modules/ionization/inputs.bf.rt | 4 ++-- Examples/Modules/ionization/inputs.rt | 4 ++-- Examples/Modules/ionization/ionization_analysis.py | 16 ++++++++++++++-- Source/Diagnostics/ParticleIO.cpp | 4 ++++ Source/Evolve/WarpXEvolveEM.cpp | 4 ++++ 6 files changed, 28 insertions(+), 8 deletions(-) (limited to 'Examples/Modules/ionization/ionization_analysis.py') diff --git a/Examples/Modules/ionization/inputs.2d b/Examples/Modules/ionization/inputs.2d index fe5e5b942..b1111c138 100644 --- a/Examples/Modules/ionization/inputs.2d +++ b/Examples/Modules/ionization/inputs.2d @@ -61,8 +61,8 @@ ions.ux = 0.0 ions.uy = 0.0 ions.uz = 0.0 ions.do_field_ionization = 1 -ions.ionization_level = 1 -ions.ionization_product = electrons +ions.ionization_initial_level = 1 +ions.ionization_product_species = electrons ions.physical_element = Si lasers.nlasers = 1 diff --git a/Examples/Modules/ionization/inputs.bf.rt b/Examples/Modules/ionization/inputs.bf.rt index 2d3b27b52..0e2e16164 100644 --- a/Examples/Modules/ionization/inputs.bf.rt +++ b/Examples/Modules/ionization/inputs.bf.rt @@ -31,8 +31,8 @@ ions.profile = constant ions.density = 1. ions.momentum_distribution_type = constant ions.do_field_ionization = 1 -ions.ionization_level = 2 -ions.ionization_product = electrons +ions.ionization_initial_level = 2 +ions.ionization_product_species = electrons ions.physical_element = N ions.do_continuous_injection=1 diff --git a/Examples/Modules/ionization/inputs.rt b/Examples/Modules/ionization/inputs.rt index 718aed407..fd8ff9aae 100644 --- a/Examples/Modules/ionization/inputs.rt +++ b/Examples/Modules/ionization/inputs.rt @@ -26,8 +26,8 @@ ions.profile = constant ions.density = 1. ions.momentum_distribution_type = constant ions.do_field_ionization = 1 -ions.ionization_level = 2 -ions.ionization_product = electrons +ions.ionization_initial_level = 2 +ions.ionization_product_species = electrons ions.physical_element = N electrons.mass = m_p diff --git a/Examples/Modules/ionization/ionization_analysis.py b/Examples/Modules/ionization/ionization_analysis.py index 4afe5c577..03fa86c16 100755 --- a/Examples/Modules/ionization/ionization_analysis.py +++ b/Examples/Modules/ionization/ionization_analysis.py @@ -1,21 +1,33 @@ #! /usr/bin/env python +""" +This script tests the result of the ionization module in WarpX. + +Input files inputs.rt and inputs.bf.rt are used to reproduce the test from +Chen, JCP, 2013, figure 2 (in the lab frame and in a boosted frame, +respectively): a plane-wave laser pulse propagates through a +uniform N2+ neutral plasma and further ionizes the Nitrogen atoms. This test +checks that, after the laser went through the plasma, ~32% of Nitrogen +ions are N5+, in agreement with theory from Chen's article. +""" + import sys import yt import numpy as np import scipy.constants as scc yt.funcs.mylog.setLevel(0) +# Open plotfile specified in command line, and get ion's ionization level. filename = sys.argv[1] - ds = yt.load( filename ) ad = ds.all_data() ilev = ad['ions', 'particle_ionization_level'].v +# Fraction of Nitrogen ions that are N5+. N5_fraction = ilev[ilev == 5].size/float(ilev.size) print("Number of ions: " + str(ilev.size)) print("Number of N5+ : " + str(ilev[ilev == 5].size)) -print("N5_fraction: " + str(N5_fraction)) +print("N5_fraction : " + str(N5_fraction)) assert ((N5_fraction > 0.30) and (N5_fraction < 0.34)) diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index 618ecc6c5..be4d87746 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -115,6 +115,10 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const if(pc->do_field_ionization){ int_names.push_back("ionization_level"); + // int_flags specifies, for each integer attribs, whether it is + // dumped to plotfiles. So far, ionization_level is the only + // integer attribs, and it is automatically dumped to plotfiles + // when ionization is on. int_flags.resize(1, 1); } diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 723d43b06..d31c9107c 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -342,6 +342,10 @@ WarpX::OneStep_sub1 (Real curtime) { // TODO: we could save some charge depositions + // Loop over species. For each ionizable species, create particles in + // product species. + mypc->doFieldIonization(); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level == 1, "Must have exactly two levels"); const int fine_lev = 1; const int coarse_lev = 0; -- cgit v1.2.3 From 8469c5419922cbdbcec31fbf1d80fdb0d88674a1 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Mon, 12 Aug 2019 11:17:59 -0700 Subject: add plotting capability in ionization rt analysis script --- Examples/Modules/ionization/ionization_analysis.py | 44 ++++++++++++++++++++++ 1 file changed, 44 insertions(+) (limited to 'Examples/Modules/ionization/ionization_analysis.py') diff --git a/Examples/Modules/ionization/ionization_analysis.py b/Examples/Modules/ionization/ionization_analysis.py index 03fa86c16..b94541f90 100755 --- a/Examples/Modules/ionization/ionization_analysis.py +++ b/Examples/Modules/ionization/ionization_analysis.py @@ -30,4 +30,48 @@ print("Number of ions: " + str(ilev.size)) print("Number of N5+ : " + str(ilev[ilev == 5].size)) print("N5_fraction : " + str(N5_fraction)) +do_plot = False +if do_plot: + import matplotlib.pyplot as plt + all_data_level_0 = ds.covering_grid(level=0,left_edge=ds.domain_left_edge, + dims=ds.domain_dimensions) + F = all_data_level_0['boxlib', 'Ex'].v.squeeze() + extent = [ ds.domain_left_edge[1], ds.domain_right_edge[1], + ds.domain_left_edge[0], ds.domain_right_edge[0] ] + ad = ds.all_data() + + # Plot ions with ionization levels + species = 'ions'; + xi = ad[species, 'particle_position_x'].v + zi = ad[species, 'particle_position_y'].v + ii = ad[species, 'particle_ionization_level'].v + plt.figure(figsize=(10,10)) + plt.subplot(211) + plt.imshow(np.abs(F), extent=extent, aspect='auto', + cmap='magma', origin='default') + plt.colorbar() + for lev in range(int(np.max(ii)+1)): + select = (ii == lev) + plt.scatter(zi[select],xi[select],s=.2, + label='ionization level: ' + str(lev)) + plt.legend() + plt.title("abs(Ex) (V/m) and ions") + plt.xlabel("z (m)") + plt.ylabel("x (m)") + plt.subplot(212) + plt.imshow(np.abs(F), extent=extent, aspect='auto', + cmap='magma', origin='default') + plt.colorbar() + + # Plot electrons + species = 'electrons'; + if species in [x[0] for x in ds.field_list]: + xe = ad[species, 'particle_position_x'].v + ze = ad[species, 'particle_position_y'].v + plt.scatter(ze,xe,s=.1,c='r',label='electrons') + plt.title("abs(Ex) (V/m) and electrons") + plt.xlabel("z (m)") + plt.ylabel("x (m)") + plt.savefig("image_ionization.pdf", bbox_inches='tight') + assert ((N5_fraction > 0.30) and (N5_fraction < 0.34)) -- cgit v1.2.3