forked from ECP-WarpX/WarpX
-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Improving Coulomb collision method for weighted-particles (ECP-WarpX#…
…5091) * using corrected weighted-particle Coulomb collision method. * adding CI test for weighted Coulomb collisions. * using n12 in UpdateMomentumPerezElastic * updating Checksum * fixing type issue. * updating checksums * checksum
- Loading branch information
1 parent
5343aa8
commit b1e7932
Showing
13 changed files
with
313 additions
and
109 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,126 @@ | ||
| #!/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 | ||
| import yt | ||
| from scipy.constants import e | ||
|
|
||
| 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 | ||
|
|
||
| # Separate macroparticles from group A (low weight) and group B (high weight) | ||
| # by sorting based on weight | ||
| sorted_indices = data['ions','particle_weight'].argsort() | ||
| sorted_wp = data['ions', 'particle_weight'][sorted_indices].value | ||
| 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 | ||
|
|
||
| # Find the index 'Npmin' that separates macroparticles from group A and group B | ||
| 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 | ||
|
|
||
| NpA = Npmin | ||
| wpA = wpmin; | ||
| NpB = Np - Npmin | ||
| wpB = wpmax; | ||
| NpAs = 0 | ||
| NpAe = Npmin | ||
| NpBs = Npmin | ||
| NpBe = Np | ||
|
|
||
| ############# | ||
|
|
||
| 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(); | ||
|
|
||
| # 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 | ||
|
|
||
| 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('TApar at 30ps error = ', error); | ||
| print('tolerance = ', tolerance); | ||
| assert error < tolerance | ||
|
|
||
| test_name = os.path.split(os.getcwd())[1] | ||
| checksumAPI.evaluate_checksum(test_name, fn) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,90 @@ | ||
| ################################# | ||
| ########## CONSTANTS ############ | ||
| ################################# | ||
|
|
||
| my_constants.nA = 1.e25 # 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_level = 0 | ||
| amr.blocking_factor = 4 | ||
| 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.do_not_deposit = 1 | ||
|
|
||
| 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 = diag1 | ||
| diag1.intervals = 600 | ||
| diag1.diag_type = Full |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,20 @@ | ||
| { | ||
| "lev=0": { | ||
| "Bx": 0.0, | ||
| "By": 0.0, | ||
| "Bz": 0.0, | ||
| "Ex": 0.0, | ||
| "Ey": 0.0, | ||
| "Ez": 0.0, | ||
| "jx": 0.0, | ||
| "jy": 0.0, | ||
| "jz": 0.0 | ||
| }, | ||
| "ions": { | ||
| "particle_momentum_x": 3.425400072687143e-16, | ||
| "particle_momentum_y": 3.421937133999805e-16, | ||
| "particle_momentum_z": 5.522701882677923e-16, | ||
| "particle_position_x": 720.0011611411148, | ||
| "particle_weight": 1.0999999999999999e+24 | ||
| } | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.