creating CI test for weighted coulomb scattering.

Examples/Tests/collision/analysis_collision_1d.py

@@ -0,0 +1,154 @@
+#!/usr/bin/env python3
+
+# Copyright 2024 Justin Angus
+#
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+#
+# This is a script that analyses the simulation results from the script `inputs_1d`. 
+# run locally: python analysis_vandb_1d.py diags/diag1000600/
+#
+# This is a 1D intra-species Coulomb scattering relaxation test consisting
+# of a low-density population streaming into a higher density population at rest.
+# Both populations belong to the same carbon12 ion species.
+# See test T1b from JCP 413 (2020) by D. Higginson, et al.
+#
+import os
+import sys
+
+import numpy as np
+from scipy.constants import epsilon_0, e
+import yt
+
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+# this will be the name of the plot file
+fn = sys.argv[1]
+ds = yt.load(fn)
+data = ds.covering_grid(level = 0, left_edge = ds.domain_left_edge, dims = ds.domain_dimensions)
+
+# carbon 12 ion (mass = 12*amu - 6*me)
+mass = 1.992100316897910e-26
+
+# compute sum of ion values
+sorted_wp = data['ions', 'particle_weight'].value
+sorted_wp.sort();
+sorted_indices = data['ions','particle_weight'].argsort()
+sorted_px = data['ions', 'particle_momentum_x'][sorted_indices].value
+sorted_py = data['ions', 'particle_momentum_y'][sorted_indices].value
+sorted_pz = data['ions', 'particle_momentum_z'][sorted_indices].value
+
+Np = len(sorted_wp)
+wpmin = sorted_wp.min();
+wpmax = sorted_wp.max();
+for i in range(len(sorted_wp)):
+    if sorted_wp[i] > wpmin:
+        Npmin = i
+        break
+
+test_tag = 'T1b'
+print(f"test_tag = {test_tag}")
+
+### use below for tests T1b and T1c, where NpA = Npmin
+if test_tag == 'T1b' or test_tag == 'T1c':
+    NpA = Npmin
+    wpA = wpmin;
+    NpB = Np - Npmin
+    wpB = wpmax;
+    NpAs = 0
+    NpAe = Npmin
+    NpBs = Npmin
+    NpBe = Np
+
+### use below for tests T1a and T1d, where NpB = Npmin
+if test_tag == 'T1a' or test_tag == 'T1d':
+    NpB = Npmin
+    wpb = wpmin;
+    NpA = Np - Npmin
+    wpA = wpmax;
+    NpBs = 0
+    NpBe = Npmin
+    NpAs = Npmin
+    NpAe = Np
+
+#############
+
+print(f"Np = {Np}, NpA = {NpA}, NpB = {NpB}")
+print(f"wpA = {wpA}, wpB = {wpB}")
+
+sorted_px_sum = np.abs(sorted_px).sum();
+sorted_py_sum = np.abs(sorted_py).sum();
+sorted_pz_sum = np.abs(sorted_pz).sum();
+sorted_wp_sum = np.abs(sorted_wp).sum();
+
+#print(f"sum(sorted_px)  =  {sorted_px_sum}")
+#print(f"sum(sorted_py)  =  {sorted_py_sum}")
+#print(f"sum(sorted_pz)  =  {sorted_pz_sum}")
+#print(f"sum(sorted_wp)  =  {sorted_wp_sum}")
+
+# compute mean velocities
+wAtot = wpA*NpA
+wBtot = wpB*NpB
+
+uBx = uBy = uBz = 0.
+for i in range(NpBs,NpBe):
+    uBx += wpB*sorted_px[i]
+    uBy += wpB*sorted_py[i]
+    uBz += wpB*sorted_pz[i]
+uBx /= (mass*wBtot) # [m/s]
+uBy /= (mass*wBtot) # [m/s]
+uBz /= (mass*wBtot) # [m/s]
+
+uAx = uAy = uAz = 0.
+for i in range(NpAs,NpAe):
+    uAx += wpA*sorted_px[i]
+    uAy += wpA*sorted_py[i]
+    uAz += wpA*sorted_pz[i]
+uAx /= (mass*wAtot) # [m/s]
+uAy /= (mass*wAtot) # [m/s]
+uAz /= (mass*wAtot) # [m/s]
+
+# compute temperatures
+TBx = TBy = TBz = 0.
+for i in range(NpBs,NpBe):
+    TBx += wpB*(sorted_px[i]/mass - uBx)**2
+    TBy += wpB*(sorted_py[i]/mass - uBy)**2
+    TBz += wpB*(sorted_pz[i]/mass - uBz)**2
+TBx *= mass/(e*wBtot)
+TBy *= mass/(e*wBtot)
+TBz *= mass/(e*wBtot)
+
+TAx = TAy = TAz = 0.
+for i in range(NpAs,NpAe):
+    TAx += wpA*(sorted_px[i]/mass - uAx)**2
+    TAy += wpA*(sorted_py[i]/mass - uAy)**2
+    TAz += wpA*(sorted_pz[i]/mass - uAz)**2
+TAx *= mass/(e*wAtot)
+TAy *= mass/(e*wAtot)
+TAz *= mass/(e*wAtot)
+
+TApar = TAz
+TAperp = (TAx + TAy)/2.0
+TA = (TAx + TAy + TAz)/3.0
+
+TBpar = TBz
+TBperp = (TBx + TBy)/2.0
+TB = (TBx + TBy + TBz)/3.0
+
+print(f"uAx = {uAx/1e3:.{5}f}; uAy = {uAy/1e3:.{5}f}; uAz = {uAz/1e3:.{5}f} [km/s]")
+print(f"uBx = {uBx/1e3:.{5}f}; uBy = {uBy/1e3:.{5}f}; uBz = {uBz/1e3:.{5}f} [km/s]")
+print(f"TA = {TA:.{5}f}; TAperp = {TAperp:.{5}f}; TApar = {TApar:.{5}f} [eV]")
+print(f"TB = {TB:.{5}f}; TBperp = {TBperp:.{5}f}; TBpar = {TBpar:.{5}f} [eV]")
+
+TApar_30ps_soln = 6.15e3 # TA parallel solution at t = 30 ps
+error = np.abs(TApar-TApar_30ps_soln)/TApar_30ps_soln
+tolerance = 0.02
+print('error = ', error);
+print('tolerance = ', tolerance);
+assert error < tolerance
+
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, fn)

Examples/Tests/collision/inputs_1d

@@ -0,0 +1,97 @@
+#################################
+########## CONSTANTS ############
+#################################
+
+my_constants.epsilon = 1.0e-12
+#
+my_constants.nA = (1.+epsilon)*1e25       # m^-3
+my_constants.NpA = 400        # m^-3
+my_constants.UA = 6.55e5      # m/s
+my_constants.TA = 500.        # eV
+#
+my_constants.nB = 1.e26       # m^-3
+my_constants.NpB = 400        # m^-3
+my_constants.UB = 0.          # m/s
+my_constants.TB = 500.        # eV
+#
+my_constants.q_c12 = 6.*q_e
+my_constants.m_c12 = 12.*m_u - 6.*m_e
+
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 600
+amr.n_cell = 180
+amr.max_grid_size = 4
+amr.blocking_factor = 4
+amr.max_level = 0
+geometry.dims = 1
+geometry.prob_lo = 0.
+geometry.prob_hi = 0.01
+
+#################################
+###### Boundary Condition #######
+#################################
+boundary.field_lo = periodic
+boundary.field_hi = periodic
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.serialize_initial_conditions = 1
+warpx.verbose = 1
+warpx.const_dt = 0.05e-12
+warpx.use_filter = 0
+
+# Do not evolve the E and B fields
+algo.maxwell_solver = none
+
+# Order of particle shape factors
+algo.particle_shape = 1
+
+#################################
+############ PLASMA #############
+#################################
+particles.species_names = ions
+
+ions.charge = q_c12
+ions.mass = m_c12
+
+ions.injection_sources = groupA groupB
+
+ions.groupA.injection_style = "NUniformPerCell"
+ions.groupA.num_particles_per_cell_each_dim = NpA
+ions.groupA.profile = constant
+ions.groupA.density = nA  # number per m^3
+ions.groupA.momentum_distribution_type = "gaussian"
+ions.groupA.uz_m = UA/clight
+ions.groupA.ux_th = sqrt(TA*q_e/m_c12)/clight
+ions.groupA.uy_th = sqrt(TA*q_e/m_c12)/clight
+ions.groupA.uz_th = sqrt(TA*q_e/m_c12)/clight
+
+ions.groupB.injection_style = "NUniformPerCell"
+ions.groupB.num_particles_per_cell_each_dim = NpB
+ions.groupB.profile = constant
+ions.groupB.density = nB  # number per m^3
+ions.groupB.momentum_distribution_type = "gaussian"
+ions.groupB.uz_m = UB/clight
+ions.groupB.ux_th = sqrt(TB*q_e/m_c12)/clight
+ions.groupB.uy_th = sqrt(TB*q_e/m_c12)/clight
+ions.groupB.uz_th = sqrt(TB*q_e/m_c12)/clight
+
+#################################
+############ COLLISION ##########
+#################################
+collisions.collision_names = collision1
+collision1.species = ions ions
+collision1.CoulombLog = 10.0
+
+# Diagnostics
+diagnostics.diags_names = diags1
+
+diags1.fields_to_plot = "none"
+diags1.species = "ions"
+diags1.intervals = 600
+diags1.diag_type = Full
+diags1.ions.variables = w z ux uy uz
+

Regression/Checksum/benchmarks_json/collisionZ.json

@@ -0,0 +1,8 @@
+{
+  "ions": {
+    "particle_momentum_x": 3.4266327177706956e-16,
+    "particle_momentum_y": 3.435391408305914e-16,
+    "particle_momentum_z": 5.531650381094717e-16,
+    "particle_weight": 1.1000000000000997e+24
+  }
+}

Regression/WarpX-tests.ini

@@ -195,6 +195,20 @@ useOMP = 1
 numthreads = 1
 analysisRoutine = Examples/Tests/collider_relevant_diags/analysis_multiple_particles.py
 
+[collisionZ]
+buildDir = .
+inputFile = Examples/Tests/collision/inputs_1d
+runtime_params =
+dim = 1
+addToCompileString =
+cmakeSetupOpts = -DWarpX_DIMS=1
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 1
+analysisRoutine = Examples/Tests/collision/analysis_collision_1d.py
+
 [collisionISO]
 buildDir = .
 inputFile = Examples/Tests/collision/inputs_3d_isotropization