diff options
| author |  Edoardo Zoni <59625522+EZoni@users.noreply.github.com>
| 2022-12-07 15:40:02 -0800 |
|---|---|---|
| committer |  GitHub <noreply@github.com>
| 2022-12-07 15:40:02 -0800 |
| commit | 4073384c7b66b1848bcc94e6c986f7d532c7da11 (patch) | |
| tree | a3a7d152098eff3f8c049638ac40b93a40551108 | |
| parent | 02447ce0f59e729865a8cbe9502bf6ca0c91e2cd (diff) | |
| download | WarpX-4073384c7b66b1848bcc94e6c986f7d532c7da11.tar.gz WarpX-4073384c7b66b1848bcc94e6c986f7d532c7da11.tar.zst WarpX-4073384c7b66b1848bcc94e6c986f7d532c7da11.zip | |
PSATD: Implement First-Order Equations (#3466)
* Implement First-Order PSATD Equations
* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci
* Fix Unused Parameter Warning
* Fix RZ Build
* Fix Normalization of G to Match PML
* Add CI Test: 3D Uniform Plasma
* Cleaning
* Update 2D CI Checksums
* Update 3D CI Checksums
* Add F,G to CI Checksums of `uniform_plasma_multiJ`
* Allow User to Choose First-Order v. Second-Order
* Add WARPX_ALWAYS_ASSERT_WITH_MESSAGE
* Rename New Class `PsatdAlgorithmFirstOrder`
* Remove Inline Comment
* Update RZ CI Checksums
* Fix inline comment
* Use auxiliary variables to avoid divisions
* Use auxiliary variables to avoid divisions
* Make `nci_psatd_stability` dir and merge inputs
* Move all Galilean tests to `nci_psatd_stability`
* Fix CI
* Fix index for backward FFT of J
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
37 files changed, 957 insertions, 292 deletions
diff --git a/Examples/Tests/galilean/analysis.py b/Examples/Tests/nci_psatd_stability/analysis_galilean.py index 9fe6bab72..9fe6bab72 100755 --- a/Examples/Tests/galilean/analysis.py +++ b/Examples/Tests/nci_psatd_stability/analysis_galilean.py diff --git a/Examples/Tests/nci_psatd_stability/analysis_multiJ.py b/Examples/Tests/nci_psatd_stability/analysis_multiJ.py new file mode 100755 index 000000000..1c68b114c --- /dev/null +++ b/Examples/Tests/nci_psatd_stability/analysis_multiJ.py @@ -0,0 +1,49 @@ +#!/usr/bin/env python3 +""" +This script is used to test the results of the multi-J PSATD +first-order equations, with one J deposition. It compares the +energy of the electric field with a given reference energy. + +The reference energy is computed by running the same test with J constant +in time, rho linear in time, and without divergence cleaning. The reference +energy corresponds to unstable results due to NCI (suppressed by the use of +both J and rho constant in time, and with divergence cleaning). +""" +import os +import sys + +import numpy as np +import scipy.constants as scc + +import yt ; yt.funcs.mylog.setLevel(0) +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +filename = sys.argv[1] + +ds = yt.load(filename) + +# yt 4.0+ has rounding issues with our domain data: +# RuntimeError: yt attempted to read outside the boundaries +# of a non-periodic domain along dimension 0. +if 'force_periodicity' in dir(ds): ds.force_periodicity() + +all_data = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) +Ex = all_data['boxlib', 'Ex'].squeeze().v +Ey = all_data['boxlib', 'Ey'].squeeze().v +Ez = all_data['boxlib', 'Ez'].squeeze().v + +# Set reference energy values, and tolerances for numerical stability and charge conservation +tol_energy = 1e-8 +energy_ref = 66e6 + +# Check numerical stability by comparing electric field energy to reference energy +energy = np.sum(scc.epsilon_0/2*(Ex**2+Ey**2+Ez**2)) +err_energy = energy / energy_ref +print('\nCheck numerical stability:') +print(f'err_energy = {err_energy}') +print(f'tol_energy = {tol_energy}') +assert(err_energy < tol_energy) + +test_name = os.path.split(os.getcwd())[1] +checksumAPI.evaluate_checksum(test_name, filename) diff --git a/Examples/Tests/galilean/inputs_2d b/Examples/Tests/nci_psatd_stability/inputs_2d index f19bfdde5..f19bfdde5 100644 --- a/Examples/Tests/galilean/inputs_2d +++ b/Examples/Tests/nci_psatd_stability/inputs_2d diff --git a/Examples/Tests/galilean/inputs_2d_hybrid b/Examples/Tests/nci_psatd_stability/inputs_2d_hybrid index a9d1ecd14..a9d1ecd14 100644 --- a/Examples/Tests/galilean/inputs_2d_hybrid +++ b/Examples/Tests/nci_psatd_stability/inputs_2d_hybrid diff --git a/Examples/Tests/galilean/inputs_3d b/Examples/Tests/nci_psatd_stability/inputs_3d index 7a0d80748..53305af12 100644 --- a/Examples/Tests/galilean/inputs_3d +++ b/Examples/Tests/nci_psatd_stability/inputs_3d @@ -1,84 +1,81 @@ -################################# -####### GENERAL PARAMETERS ###### -################################# -max_step = 300 - -amr.n_cell = 32 32 32 -warpx.numprocs = 1 1 2 -amr.max_level = 0 -psatd.v_galilean = 0. 0. 0.99498743710662 - -geometry.dims = 3 -geometry.prob_lo = -9.67 -9.67 -9.67 -geometry.prob_hi = 9.67 9.67 9.67 - -################################# -###### Boundary condition ####### -################################# -boundary.field_lo = periodic periodic periodic -boundary.field_hi = periodic periodic periodic - -################################# -############ NUMERICS ########### -################################# -warpx.verbose = 1 - +# algo algo.current_deposition = direct algo.maxwell_solver = psatd algo.particle_pusher = vay - -warpx.cfl = 1. algo.particle_shape = 3 -################################# -############ PLASMA ############# -################################# -particles.species_names = electrons ions +# amr +amr.max_level = 0 +amr.n_cell = 32 32 32 -warpx.do_nodal = 1 -warpx.use_filter = 1 +# boundary +boundary.field_hi = periodic periodic periodic +boundary.field_lo = periodic periodic periodic -psatd.nox = 8 -psatd.noy = 8 -psatd.noz = 8 +# diag1 +diag1.diag_type = Full +diag1.intervals = 300 +# diagnostics +diagnostics.diags_names = diag1 + +# electrons electrons.charge = -q_e -electrons.mass = m_e +electrons.density = 282197938148984.7 electrons.injection_style = "NUniformPerCell" +electrons.mass = m_e +electrons.momentum_distribution_type = "gaussian" electrons.num_particles_per_cell_each_dim = 1 1 1 electrons.profile = constant -electrons.density = 282197938148984.7 -electrons.momentum_distribution_type = "gaussian" -electrons.uz_m = 9.9498743710661994 -electrons.xmin = -9.67 -electrons.xmax = 9.67 -electrons.ymin = -9.67 -electrons.ymax = 9.67 -electrons.zmin = -9.67 -electrons.zmax = 9.67 electrons.ux_th = 0.0001 electrons.uy_th = 0.0001 +electrons.uz_m = 9.9498743710661994 electrons.uz_th = 0.0001 +electrons.xmax = 9.67 +electrons.xmin = -9.67 +electrons.ymax = 9.67 +electrons.ymin = -9.67 +electrons.zmax = 9.67 +electrons.zmin = -9.67 + +# geometry +geometry.dims = 3 +geometry.prob_hi = 9.67 9.67 9.67 +geometry.prob_lo = -9.67 -9.67 -9.67 +# ions ions.charge = q_e -ions.mass = m_p +ions.density = 282197938148984.7 ions.injection_style = "NUniformPerCell" +ions.mass = m_p +ions.momentum_distribution_type = "gaussian" ions.num_particles_per_cell_each_dim = 1 1 1 ions.profile = constant -ions.density = 282197938148984.7 -ions.momentum_distribution_type = "gaussian" -ions.uz_m = 9.9498743710661994 -ions.xmin = -9.67 -ions.xmax = 9.67 -ions.ymin = -9.67 -ions.ymax = 9.67 -ions.zmin = -9.67 -ions.zmax = 9.67 ions.ux_th = 0.0001 ions.uy_th = 0.0001 +ions.uz_m = 9.9498743710661994 ions.uz_th = 0.0001 +ions.xmax = 9.67 +ions.xmin = -9.67 +ions.ymax = 9.67 +ions.ymin = -9.67 +ions.zmax = 9.67 +ions.zmin = -9.67 -# Diagnostics -diagnostics.diags_names = diag1 -diag1.intervals = 100 -diag1.diag_type = Full +# max_step +max_step = 300 + +# particles +particles.species_names = electrons ions + +# psatd +psatd.nox = 8 +psatd.noy = 8 +psatd.noz = 8 + +# warpx +warpx.cfl = 1. +warpx.do_nodal = 1 +warpx.numprocs = 1 1 2 +warpx.use_filter = 1 +warpx.verbose = 1 diff --git a/Examples/Tests/averaged_galilean/inputs_avg_2d b/Examples/Tests/nci_psatd_stability/inputs_avg_2d index 2a2611906..2a2611906 100644 --- a/Examples/Tests/averaged_galilean/inputs_avg_2d +++ b/Examples/Tests/nci_psatd_stability/inputs_avg_2d diff --git a/Examples/Tests/averaged_galilean/inputs_avg_3d b/Examples/Tests/nci_psatd_stability/inputs_avg_3d index bb76e4c47..bb76e4c47 100644 --- a/Examples/Tests/averaged_galilean/inputs_avg_3d +++ b/Examples/Tests/nci_psatd_stability/inputs_avg_3d diff --git a/Examples/Tests/galilean/inputs_rz b/Examples/Tests/nci_psatd_stability/inputs_rz index 772711ddc..772711ddc 100644 --- a/Examples/Tests/galilean/inputs_rz +++ b/Examples/Tests/nci_psatd_stability/inputs_rz diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_multiJ.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_multiJ.json index b0c362f0f..c342d9dca 100644 --- a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_multiJ.json +++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_multiJ.json @@ -1,29 +1,29 @@ { "electrons": { - "particle_momentum_x": 5.664739488600762e-20, + "particle_momentum_x": 5.663705657675969e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 5.664739488600764e-20, - "particle_position_x": 0.6553599999999999, + "particle_momentum_z": 5.663705657675969e-20, + "particle_position_x": 0.65536, "particle_position_y": 0.65536, "particle_weight": 3200000000000000.5 }, "lev=0": { "Bx": 0.0, - "By": 3.4900393205053586, + "By": 3.4892704618136277, "Bz": 0.0, - "Ex": 3771422651410.755, + "Ex": 3771082786646.7104, "Ey": 0.0, - "Ez": 3771422651410.742, - "jx": 1.0095457953459832e+16, + "Ez": 3771082786646.702, + "jx": 1.0093631772735916e+16, "jy": 0.0, - "jz": 1.0095457953459836e+16 + "jz": 1.0093631772735912e+16 }, "positrons": { - "particle_momentum_x": 5.664739488600754e-20, + "particle_momentum_x": 5.663705657675971e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 5.664739488600756e-20, - "particle_position_x": 0.6553599999999999, - "particle_position_y": 0.6553599999999999, + "particle_momentum_z": 5.663705657675969e-20, + "particle_position_x": 0.65536, + "particle_position_y": 0.65536, "particle_weight": 3200000000000000.5 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_multiJ_nodal.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_multiJ_nodal.json index 66c8e3e80..2aaa9e0f2 100644 --- a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_multiJ_nodal.json +++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_multiJ_nodal.json @@ -1,29 +1,29 @@ { "electrons": { - "particle_momentum_x": 5.668522616618711e-20, + "particle_momentum_x": 5.668409581183492e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 5.668522616618715e-20, - "particle_position_x": 0.6553600000002356, - "particle_position_y": 0.6553600000002355, + "particle_momentum_z": 5.668409581183495e-20, + "particle_position_x": 0.6553600000004717, + "particle_position_y": 0.6553600000004718, "particle_weight": 3200000000000000.5 }, "lev=0": { "Bx": 0.0, - "By": 5.6351165293218966, + "By": 5.634624433795568, "Bz": 0.0, - "Ex": 3747153697353.926, + "Ex": 3747068321640.6587, "Ey": 0.0, - "Ez": 3747153697353.9287, - "jx": 1.0088631639558242e+16, + "Ez": 3747068321640.659, + "jx": 1.0088430505673088e+16, "jy": 0.0, - "jz": 1.0088631639558248e+16 + "jz": 1.0088430505673096e+16 }, "positrons": { - "particle_momentum_x": 5.66852261661871e-20, + "particle_momentum_x": 5.668409581183492e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 5.668522616618714e-20, - "particle_position_x": 0.6553600000002356, - "particle_position_y": 0.6553600000002356, + "particle_momentum_z": 5.668409581183495e-20, + "particle_position_x": 0.6553600000004718, + "particle_position_y": 0.6553600000004717, "particle_weight": 3200000000000000.5 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_multiJ.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_multiJ.json index c76d7cfc5..487d4a898 100644 --- a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_multiJ.json +++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_multiJ.json @@ -1,28 +1,28 @@ { "electrons": { - "particle_momentum_x": 9.629015522300135e-20, - "particle_position_x": 2.621440000009873, - "particle_position_y": 2.621440000009873, - "particle_position_z": 2.6214399999999998, + "particle_momentum_x": 9.633869745818886e-20, + "particle_position_x": 2.621440000001177, + "particle_position_y": 2.6214400000011784, + "particle_position_z": 2.6214400000000007, "particle_weight": 128000000000.00002 }, "lev=0": { - "Bx": 79.96476923345703, - "By": 79.96476923350225, - "Bz": 79.96690317049361, - "Ex": 84753887916472.72, - "Ey": 84753887916472.66, - "Ez": 84753877853695.67, - "jx": 6.081254778189634e+16, - "jy": 6.081254778189637e+16, - "jz": 6.081251943036953e+16, + "Bx": 80.96035538655111, + "By": 80.96035538657691, + "Bz": 80.96271445263956, + "Ex": 84777489275096.88, + "Ey": 84777489275096.88, + "Ez": 84777485856239.4, + "jx": 6.08447015360442e+16, + "jy": 6.084470153604425e+16, + "jz": 6.084470085554113e+16, "part_per_cell": 524288.0, - "rho": 703417424.2676101 + "rho": 703546536.8089281 }, "positrons": { - "particle_momentum_z": 9.629011306229332e-20, - "particle_position_x": 2.621440000009873, - "particle_position_y": 2.621440000009873, - "particle_position_z": 2.6214399999999998 + "particle_momentum_z": 9.63386961193585e-20, + "particle_position_x": 2.621440000001177, + "particle_position_y": 2.621440000001179, + "particle_position_z": 2.6214400000000007 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_multiJ_nodal.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_multiJ_nodal.json index 1f89c4dcb..e97556e55 100644 --- a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_multiJ_nodal.json +++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_multiJ_nodal.json @@ -1,28 +1,28 @@ { "electrons": { - "particle_momentum_x": 9.3282651765877e-20, + "particle_momentum_x": 9.320515980761094e-20, "particle_position_x": 2.6214400000000015, "particle_position_y": 2.621440000000001, - "particle_position_z": 2.621440000000001, + "particle_position_z": 2.6214400000000007, "particle_weight": 128000000000.00002 }, "lev=0": { - "Bx": 17.338468210649435, - "By": 17.338468210679473, - "Bz": 17.338468210708463, - "Ex": 86130544037694.12, - "Ey": 86130544037694.16, - "Ez": 86130544037694.16, - "jx": 5.808322546347548e+16, - "jy": 5.80832254634755e+16, - "jz": 5.8083225463475464e+16, + "Bx": 17.319790986722108, + "By": 17.319790986747606, + "Bz": 17.319790986759728, + "Ex": 86079985956930.31, + "Ey": 86079985956930.3, + "Ez": 86079985956930.31, + "jx": 5.8033630197761816e+16, + "jy": 5.803363019776185e+16, + "jz": 5.803363019776182e+16, "part_per_cell": 524288.0, - "rho": 721143170.1131016 + "rho": 720717041.7116292 }, "positrons": { - "particle_momentum_z": 9.328265176587699e-20, + "particle_momentum_z": 9.320515980761093e-20, "particle_position_x": 2.6214400000000015, "particle_position_y": 2.621440000000001, - "particle_position_z": 2.621440000000001 + "particle_position_z": 2.6214400000000007 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/multi_J_rz_psatd.json b/Regression/Checksum/benchmarks_json/multi_J_rz_psatd.json index 458117917..2a9620903 100644 --- a/Regression/Checksum/benchmarks_json/multi_J_rz_psatd.json +++ b/Regression/Checksum/benchmarks_json/multi_J_rz_psatd.json @@ -9,31 +9,31 @@ "particle_weight": 6241484108.424456 }, "lev=0": { - "By": 24912.66260587033, - "Ex": 4667306677660.305, - "Ez": 4307437890419.4253, - "jx": 362735950033724.2, - "jz": 1937267340131275.2, - "rho": 5308546.075566203, - "rho_driver": 6288266.101815152, - "rho_plasma_e": 49569864.00850832, - "rho_plasma_p": 50769174.61530346 + "By": 24045.34926330333, + "Ex": 4530500183998.0205, + "Ez": 4297045713383.818, + "jx": 361149490291532.8, + "jz": 1805428826325930.2, + "rho": 4895064.444869195, + "rho_driver": 6288266.101815153, + "rho_plasma_e": 49569825.13450765, + "rho_plasma_p": 50769176.974483095 }, "plasma_e": { - "particle_momentum_x": 6.65811141013141e-20, - "particle_momentum_y": 6.738987045495091e-20, - "particle_momentum_z": 2.846571109435123e-20, - "particle_position_x": 1.1423367495493797, - "particle_position_y": 0.6139715590509269, - "particle_theta": 20188.939948727297, + "particle_momentum_x": 6.649609416554171e-20, + "particle_momentum_y": 6.724424134488497e-20, + "particle_momentum_z": 2.81433599356851e-20, + "particle_position_x": 1.14233458508442, + "particle_position_y": 0.6140029351346352, + "particle_theta": 20188.939948727293, "particle_weight": 1002457942911.3788 }, "plasma_p": { - "particle_momentum_x": 6.640192720118765e-20, - "particle_momentum_y": 6.767588557043428e-20, - "particle_momentum_z": 5.58476124317102e-20, - "particle_position_x": 1.1365201600226575, - "particle_position_y": 0.6152066982817419, + "particle_momentum_x": 6.635739630451454e-20, + "particle_momentum_y": 6.761235868190358e-20, + "particle_momentum_z": 5.475884783890083e-20, + "particle_position_x": 1.1365201524804165, + "particle_position_y": 0.6152066555308389, "particle_theta": 20286.92798337582, "particle_weight": 1002457942911.3788 } diff --git a/Regression/Checksum/benchmarks_json/uniform_plasma_multiJ.json b/Regression/Checksum/benchmarks_json/uniform_plasma_multiJ.json new file mode 100644 index 000000000..e215428fe --- /dev/null +++ b/Regression/Checksum/benchmarks_json/uniform_plasma_multiJ.json @@ -0,0 +1,35 @@ +{ + "electrons": { + "particle_momentum_x": 2.4258140600451355e-21, + "particle_momentum_y": 2.4115944515722047e-21, + "particle_momentum_z": 8.903860561145549e-17, + "particle_position_x": 158433.52093622342, + "particle_position_y": 158432.7590925775, + "particle_position_z": 158433.27723231513, + "particle_weight": 2.041377132710917e+18 + }, + "ions": { + "particle_momentum_x": 1.3150385606388975e-18, + "particle_momentum_y": 1.3042988334363742e-18, + "particle_momentum_z": 1.6348805659682514e-13, + "particle_position_x": 158433.5805512471, + "particle_position_y": 158432.80155515939, + "particle_position_z": 158433.27779468897, + "particle_weight": 2.041377132710917e+18 + }, + "lev=0": { + "Bx": 0.054117593911342535, + "By": 0.053964421360199535, + "Bz": 0.0005552030959216721, + "divE": 18542936.02490211, + "Ex": 16260950.356941158, + "Ey": 16306834.079553638, + "Ez": 5784750.2181154955, + "F": 2.003020981926516e-02, + "G": 1.349588434641251e+03, + "jx": 984.2101557440694, + "jy": 1030.1973073214422, + "jz": 20803.708075753144, + "rho": 8.4139043087009e-05 + } +} diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index c3b1ad67a..57930ea44 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -412,7 +412,7 @@ analysisOutputImage = langmuir_multi_analysis.png [Langmuir_multi_psatd_multiJ] buildDir = . inputFile = Examples/Tests/langmuir/inputs_3d_multi_rt -runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.5773502691896258 algo.current_deposition=direct psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium +runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.5773502691896258 algo.current_deposition=direct psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.solution_type=first-order psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium dim = 3 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_PSATD=ON @@ -431,7 +431,7 @@ analysisOutputImage = Langmuir_multi_psatd_multiJ.png [Langmuir_multi_psatd_multiJ_nodal] buildDir = . inputFile = Examples/Tests/langmuir/inputs_3d_multi_rt -runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.5773502691896258 algo.current_deposition=direct psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium warpx.do_nodal=1 +runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.5773502691896258 algo.current_deposition=direct psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.solution_type=first-order psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium warpx.do_nodal=1 dim = 3 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_PSATD=ON @@ -697,7 +697,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png [Langmuir_multi_2d_psatd_multiJ] buildDir = . inputFile = Examples/Tests/langmuir/inputs_2d_multi_rt -runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.7071067811865475 psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium +runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.7071067811865475 psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.solution_type=first-order psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium dim = 2 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON @@ -716,7 +716,7 @@ analysisOutputImage = Langmuir_multi_2d_psatd_multiJ.png [Langmuir_multi_2d_psatd_multiJ_nodal] buildDir = . inputFile = Examples/Tests/langmuir/inputs_2d_multi_rt -runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.7071067811865475 psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium warpx.do_nodal=1 +runtime_params = algo.maxwell_solver=psatd warpx.cfl=0.7071067811865475 psatd.update_with_rho=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=2 psatd.solution_type=first-order psatd.J_in_time=linear warpx.abort_on_warning_threshold=medium warpx.do_nodal=1 dim = 2 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON @@ -2517,7 +2517,7 @@ analysisRoutine = Examples/Tests/particle_fields_diags/analysis_particle_diags_s [galilean_2d_psatd] buildDir = . -inputFile = Examples/Tests/galilean/inputs_2d +inputFile = Examples/Tests/nci_psatd_stability/inputs_2d runtime_params = warpx.do_nodal=1 algo.current_deposition=direct psatd.current_correction=0 warpx.abort_on_warning_threshold=medium dim = 2 addToCompileString = USE_PSATD=TRUE @@ -2531,11 +2531,11 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons ions -analysisRoutine = Examples/Tests/galilean/analysis.py +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_galilean.py [galilean_2d_psatd_current_correction_psb] buildDir = . -inputFile = Examples/Tests/galilean/inputs_2d +inputFile = Examples/Tests/nci_psatd_stability/inputs_2d runtime_params = psatd.periodic_single_box_fft=1 psatd.update_with_rho=0 psatd.current_correction=1 diag1.fields_to_plot=Ex Ey Ez Bx By Bz jx jy jz rho divE dim = 2 addToCompileString = USE_PSATD=TRUE @@ -2549,11 +2549,11 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons ions -analysisRoutine = Examples/Tests/galilean/analysis.py +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_galilean.py [galilean_2d_psatd_current_correction] buildDir = . -inputFile = Examples/Tests/galilean/inputs_2d +inputFile = Examples/Tests/nci_psatd_stability/inputs_2d runtime_params = psatd.periodic_single_box_fft=0 psatd.update_with_rho=0 psatd.current_correction=1 diag1.fields_to_plot=Ex Ey Ez Bx By Bz jx jy jz rho divE amr.max_grid_size=64 amr.blocking_factor=64 dim = 2 addToCompileString = USE_PSATD=TRUE @@ -2567,11 +2567,11 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons ions -analysisRoutine = Examples/Tests/galilean/analysis.py +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_galilean.py [galilean_2d_psatd_hybrid] buildDir = . -inputFile = Examples/Tests/galilean/inputs_2d_hybrid +inputFile = Examples/Tests/nci_psatd_stability/inputs_2d_hybrid runtime_params = psatd.current_correction=0 warpx.abort_on_warning_threshold=medium dim = 2 addToCompileString = USE_PSATD=TRUE @@ -2607,7 +2607,7 @@ analysisRoutine = Examples/analysis_default_regression.py [galilean_rz_psatd] buildDir = . -inputFile = Examples/Tests/galilean/inputs_rz +inputFile = Examples/Tests/nci_psatd_stability/inputs_rz runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_initial_conditions=1 electrons.random_theta=0 ions.random_theta=0 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium dim = 2 addToCompileString = USE_RZ=TRUE USE_PSATD=TRUE BLAS_LIB=-lblas LAPACK_LIB=-llapack @@ -2621,11 +2621,11 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons ions -analysisRoutine = Examples/Tests/galilean/analysis.py +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_galilean.py [galilean_rz_psatd_current_correction_psb] buildDir = . -inputFile = Examples/Tests/galilean/inputs_rz +inputFile = Examples/Tests/nci_psatd_stability/inputs_rz runtime_params = psatd.periodic_single_box_fft=1 psatd.current_correction=1 electrons.random_theta=0 ions.random_theta=0 dim = 2 addToCompileString = USE_RZ=TRUE USE_PSATD=TRUE BLAS_LIB=-lblas LAPACK_LIB=-llapack @@ -2639,11 +2639,11 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons ions -analysisRoutine = Examples/Tests/galilean/analysis.py +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_galilean.py [galilean_rz_psatd_current_correction] buildDir = . -inputFile = Examples/Tests/galilean/inputs_rz +inputFile = Examples/Tests/nci_psatd_stability/inputs_rz runtime_params = psatd.periodic_single_box_fft=0 psatd.current_correction=1 electrons.random_theta=0 ions.random_theta=0 amr.max_grid_size=32 amr.blocking_factor=32 dim = 2 addToCompileString = USE_RZ=TRUE USE_PSATD=TRUE BLAS_LIB=-lblas LAPACK_LIB=-llapack @@ -2657,12 +2657,12 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons ions -analysisRoutine = Examples/Tests/galilean/analysis.py +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_galilean.py [galilean_3d_psatd] buildDir = . -inputFile = Examples/Tests/galilean/inputs_3d -runtime_params = warpx.do_nodal=1 algo.current_deposition=direct psatd.current_correction=0 warpx.abort_on_warning_threshold=medium +inputFile = Examples/Tests/nci_psatd_stability/inputs_3d +runtime_params = psatd.v_galilean=0. 0. 0.99498743710662 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium dim = 3 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_PSATD=ON @@ -2675,12 +2675,12 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons ions -analysisRoutine = Examples/Tests/galilean/analysis.py +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_galilean.py [galilean_3d_psatd_current_correction_psb] buildDir = . -inputFile = Examples/Tests/galilean/inputs_3d -runtime_params = warpx.numprocs=1 1 1 psatd.periodic_single_box_fft=1 psatd.update_with_rho=0 psatd.current_correction=1 diag1.fields_to_plot=Ex Ey Ez Bx By Bz jx jy jz rho divE +inputFile = Examples/Tests/nci_psatd_stability/inputs_3d +runtime_params = psatd.v_galilean=0. 0. 0.99498743710662 warpx.numprocs=1 1 1 psatd.periodic_single_box_fft=1 psatd.update_with_rho=0 psatd.current_correction=1 diag1.fields_to_plot=Ex Ey Ez Bx By Bz jx jy jz rho divE dim = 3 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_PSATD=ON @@ -2693,12 +2693,12 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons ions -analysisRoutine = Examples/Tests/galilean/analysis.py +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_galilean.py [galilean_3d_psatd_current_correction] buildDir = . -inputFile = Examples/Tests/galilean/inputs_3d -runtime_params = warpx.numprocs=1 1 2 psatd.periodic_single_box_fft=0 psatd.update_with_rho=0 psatd.current_correction=1 diag1.fields_to_plot=Ex Ey Ez Bx By Bz jx jy jz rho divE +inputFile = Examples/Tests/nci_psatd_stability/inputs_3d +runtime_params = psatd.v_galilean=0. 0. 0.99498743710662 warpx.numprocs=1 1 2 psatd.periodic_single_box_fft=0 psatd.update_with_rho=0 psatd.current_correction=1 diag1.fields_to_plot=Ex Ey Ez Bx By Bz jx jy jz rho divE dim = 3 addToCompileString = USE_PSATD=TRUE cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_PSATD=ON @@ -2711,11 +2711,11 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons ions -analysisRoutine = Examples/Tests/galilean/analysis.py +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_galilean.py [averaged_galilean_2d_psatd] buildDir = . -inputFile = Examples/Tests/averaged_galilean/inputs_avg_2d +inputFile = Examples/Tests/nci_psatd_stability/inputs_avg_2d runtime_params = psatd.current_correction=0 warpx.abort_on_warning_threshold=medium dim = 2 addToCompileString = USE_PSATD=TRUE @@ -2729,11 +2729,11 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons ions -analysisRoutine = Examples/Tests/galilean/analysis.py +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_galilean.py [averaged_galilean_2d_psatd_hybrid] buildDir = . -inputFile = Examples/Tests/averaged_galilean/inputs_avg_2d +inputFile = Examples/Tests/nci_psatd_stability/inputs_avg_2d runtime_params = amr.max_grid_size_x = 128 amr.max_grid_size_y = 64 warpx.do_nodal = 0 algo.field_gathering = momentum-conserving interpolation.field_centering_nox = 8 interpolation.field_centering_noz = 8 warpx.do_current_centering = 1 interpolation.current_centering_nox = 8 interpolation.current_centering_noz = 8 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium dim = 2 addToCompileString = USE_PSATD=TRUE @@ -2747,11 +2747,11 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons ions -analysisRoutine = Examples/Tests/galilean/analysis.py +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_galilean.py [averaged_galilean_3d_psatd] buildDir = . -inputFile = Examples/Tests/averaged_galilean/inputs_avg_3d +inputFile = Examples/Tests/nci_psatd_stability/inputs_avg_3d runtime_params = psatd.current_correction=0 warpx.abort_on_warning_threshold=medium dim = 3 addToCompileString = USE_PSATD=TRUE @@ -2765,11 +2765,11 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons ions -analysisRoutine = Examples/Tests/galilean/analysis.py +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_galilean.py [averaged_galilean_3d_psatd_hybrid] buildDir = . -inputFile = Examples/Tests/averaged_galilean/inputs_avg_3d +inputFile = Examples/Tests/nci_psatd_stability/inputs_avg_3d runtime_params = warpx.do_nodal = 0 algo.field_gathering = momentum-conserving interpolation.field_centering_nox = 8 interpolation.field_centering_noy = 8 interpolation.field_centering_noz = 8 warpx.do_current_centering = 1 interpolation.current_centering_nox = 8 interpolation.current_centering_noy = 8 interpolation.current_centering_noz = 8 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium dim = 3 addToCompileString = USE_PSATD=TRUE @@ -2783,7 +2783,7 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = electrons ions -analysisRoutine = Examples/Tests/galilean/analysis.py +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_galilean.py [multi_J_rz_psatd] buildDir = . @@ -3748,3 +3748,21 @@ numthreads = 1 compileTest = 0 doVis = 0 analysisRoutine = Examples/Tests/btd_rz/analysis_BTD_laser_antenna.py + +[uniform_plasma_multiJ] +buildDir = . +inputFile = Examples/Tests/nci_psatd_stability/inputs_3d +runtime_params = psatd.solution_type=first-order psatd.J_in_time=constant psatd.rho_in_time=constant warpx.do_dive_cleaning=1 warpx.do_divb_cleaning=1 warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=1 diag1.fields_to_plot=Bx By Bz divE Ex Ey Ez F G jx jy jz rho warpx.abort_on_warning_threshold=medium +dim = 3 +addToCompileString = USE_PSATD=TRUE +cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_PSATD=ON +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = electrons ions +analysisRoutine = Examples/Tests/nci_psatd_stability/analysis_multiJ.py diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index a1ce21a91..e519298ca 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -133,7 +133,7 @@ public: int ncell, int delta, amrex::IntVect ref_ratio, amrex::Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, int do_moving_window, int pml_has_particles, int do_pml_in_domain, - const int J_in_time, const int rho_in_time, + const int psatd_solution_type, const int J_in_time, const int rho_in_time, const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning, const amrex::IntVect& fill_guards_fields, const amrex::IntVect& fill_guards_current, diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index c7b773123..5e03822f1 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -548,7 +548,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri int ncell, int delta, amrex::IntVect ref_ratio, Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, int do_moving_window, int /*pml_has_particles*/, int do_pml_in_domain, - const int J_in_time, const int rho_in_time, + const int psatd_solution_type, const int J_in_time, const int rho_in_time, const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning, const amrex::IntVect& fill_guards_fields, const amrex::IntVect& fill_guards_current, @@ -724,7 +724,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { #ifndef WARPX_USE_PSATD - amrex::ignore_unused(lev, dt, J_in_time, rho_in_time); + amrex::ignore_unused(lev, dt, psatd_solution_type, J_in_time, rho_in_time); # if(AMREX_SPACEDIM!=3) amrex::ignore_unused(noy_fft); # endif @@ -745,7 +745,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri spectral_solver_fp = std::make_unique<SpectralSolver>(lev, realspace_ba, dm, nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, dx, dt, in_pml, periodic_single_box, update_with_rho, - fft_do_time_averaging, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning); + fft_do_time_averaging, psatd_solution_type, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning); #endif } @@ -852,7 +852,7 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri spectral_solver_cp = std::make_unique<SpectralSolver>(lev, realspace_cba, cdm, nox_fft, noy_fft, noz_fft, do_nodal, v_galilean_zero, v_comoving_zero, cdx, dt, in_pml, periodic_single_box, update_with_rho, - fft_do_time_averaging, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning); + fft_do_time_averaging, psatd_solution_type, J_in_time, rho_in_time, m_dive_cleaning, m_divb_cleaning); #endif } } diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index e8e6a025b..8abfa0e7f 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -549,14 +549,14 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time) // 3) Deposit rho (in rho_new, since it will be moved during the loop) // (after checking that pointer to rho_fp on MR level 0 is not null) - if (rho_fp[0]) + if (rho_fp[0] && rho_in_time == RhoInTime::Linear) { // Deposit rho at relative time -dt // (dt[0] denotes the time step on mesh refinement level 0) mypc->DepositCharge(rho_fp, -dt[0]); // Filter, exchange boundary, and interpolate across levels SyncRho(); - // Forward FFT of rho_new + // Forward FFT of rho PSATDForwardTransformRho(rho_fp, rho_cp, 0, 1); } @@ -587,17 +587,14 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time) // Loop over multi-J depositions for (int i_depose = 0; i_depose < n_loop; i_depose++) { - // Move J deposited previously, from new to old - if (J_in_time == JInTime::Linear) - { - PSATDMoveJNewToJOld(); - } + // Move J from new to old if J is linear in time + if (J_in_time == JInTime::Linear) PSATDMoveJNewToJOld(); const amrex::Real t_depose_current = (J_in_time == JInTime::Linear) ? (i_depose-n_depose+1)*sub_dt : (i_depose-n_depose+0.5_rt)*sub_dt; - // TODO Update this when rho quadratic in time is implemented - const amrex::Real t_depose_charge = (i_depose-n_depose+1)*sub_dt; + const amrex::Real t_depose_charge = (rho_in_time == RhoInTime::Linear) ? + (i_depose-n_depose+1)*sub_dt : (i_depose-n_depose+0.5_rt)*sub_dt; // Deposit new J at relative time t_depose_current with time step dt // (dt[0] denotes the time step on mesh refinement level 0) @@ -616,14 +613,14 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time) // (after checking that pointer to rho_fp on MR level 0 is not null) if (rho_fp[0]) { - // Move rho deposited previously, from new to old - PSATDMoveRhoNewToRhoOld(); + // Move rho from new to old if rho is linear in time + if (rho_in_time == RhoInTime::Linear) PSATDMoveRhoNewToRhoOld(); // Deposit rho at relative time t_depose_charge mypc->DepositCharge(rho_fp, t_depose_charge); // Filter, exchange boundary, and interpolate across levels SyncRho(); - // Forward FFT of rho_new + // Forward FFT of rho PSATDForwardTransformRho(rho_fp, rho_cp, 0, 1); } diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt index a370d4b2d..912ed47c4 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt @@ -1,5 +1,6 @@ target_sources(WarpX PRIVATE + PsatdAlgorithmFirstOrder.cpp PsatdAlgorithmJConstantInTime.cpp PsatdAlgorithmJLinearInTime.cpp PsatdAlgorithmPml.cpp diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package index c798ffb01..40f9d0e9a 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -1,4 +1,5 @@ CEXE_sources += SpectralBaseAlgorithm.cpp +CEXE_sources += PsatdAlgorithmFirstOrder.cpp CEXE_sources += PsatdAlgorithmJConstantInTime.cpp CEXE_sources += PsatdAlgorithmJLinearInTime.cpp CEXE_sources += PsatdAlgorithmPml.cpp diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp index 1d6248f8d..dfd7ffe92 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp @@ -103,9 +103,9 @@ PsatdAlgorithmComoving::pushSpectralFields (SpectralFieldData& f) const const Complex Bz_old = fields(i,j,k,Idx.Bz); // Shortcuts for the values of J and rho - const Complex Jx = fields(i,j,k,Idx.Jx); - const Complex Jy = fields(i,j,k,Idx.Jy); - const Complex Jz = fields(i,j,k,Idx.Jz); + const Complex Jx = fields(i,j,k,Idx.Jx_mid); + const Complex Jy = fields(i,j,k,Idx.Jy_mid); + const Complex Jz = fields(i,j,k,Idx.Jz_mid); const Complex rho_old = fields(i,j,k,Idx.rho_old); const Complex rho_new = fields(i,j,k,Idx.rho_new); @@ -447,9 +447,9 @@ void PsatdAlgorithmComoving::CurrentCorrection (SpectralFieldData& field_data) amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Shortcuts for the values of J and rho - const Complex Jx = fields(i,j,k,Idx.Jx); - const Complex Jy = fields(i,j,k,Idx.Jy); - const Complex Jz = fields(i,j,k,Idx.Jz); + const Complex Jx = fields(i,j,k,Idx.Jx_mid); + const Complex Jy = fields(i,j,k,Idx.Jy_mid); + const Complex Jz = fields(i,j,k,Idx.Jz_mid); const Complex rho_old = fields(i,j,k,Idx.rho_old); const Complex rho_new = fields(i,j,k,Idx.rho_new); @@ -482,15 +482,15 @@ void PsatdAlgorithmComoving::CurrentCorrection (SpectralFieldData& field_data) const Complex theta = amrex::exp(- I * k_dot_v * dt * 0.5_rt); const Complex den = 1._rt - theta * theta; - fields(i,j,k,Idx.Jx) = Jx - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kx_mod / (knorm_mod * knorm_mod); - fields(i,j,k,Idx.Jy) = Jy - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * ky_mod / (knorm_mod * knorm_mod); - fields(i,j,k,Idx.Jz) = Jz - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kz_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jx_mid) = Jx - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kx_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jy_mid) = Jy - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * ky_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jz_mid) = Jz - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kz_mod / (knorm_mod * knorm_mod); } else { - fields(i,j,k,Idx.Jx) = Jx - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kx_mod / (knorm_mod * knorm_mod); - fields(i,j,k,Idx.Jy) = Jy - (kmod_dot_J - I * (rho_new - rho_old) / dt) * ky_mod / (knorm_mod * knorm_mod); - fields(i,j,k,Idx.Jz) = Jz - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kz_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jx_mid) = Jx - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kx_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jy_mid) = Jy - (kmod_dot_J - I * (rho_new - rho_old) / dt) * ky_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jz_mid) = Jz - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kz_mod / (knorm_mod * knorm_mod); } } }); diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.H new file mode 100644 index 000000000..c90b19e1a --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.H @@ -0,0 +1,100 @@ +/* Copyright 2019 + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PSATD_ALGORITHM_FIRST_ORDER_H_ +#define WARPX_PSATD_ALGORITHM_FIRST_ORDER_H_ + +#include "FieldSolver/SpectralSolver/SpectralFieldData.H" +#include "FieldSolver/SpectralSolver/SpectralKSpace.H" +#include "SpectralBaseAlgorithm.H" + +#include <AMReX_Array.H> +#include <AMReX_Config.H> +#include <AMReX_REAL.H> + +#include <AMReX_BaseFwd.H> + +#include <array> +#include <memory> + +#if WARPX_USE_PSATD +/* + * \brief Class that updates the fields in spectral space according to the first-order PSATD equations. + */ +class PsatdAlgorithmFirstOrder : public SpectralBaseAlgorithm +{ + public: + + /** + * \brief Constructor of the class PsatdAlgorithmFirstOrder + * + * \param[in] spectral_kspace spectral space + * \param[in] dm distribution mapping + * \param[in] spectral_index object containing indices to access data in spectral space + * \param[in] norder_x order of the spectral solver along x + * \param[in] norder_y order of the spectral solver along y + * \param[in] norder_z order of the spectral solver along z + * \param[in] nodal whether the E and B fields are defined on a fully nodal grid or a Yee grid + * \param[in] dt time step of the simulation + * \param[in] div_cleaning whether to use divergence correction for both E and B (thus, F and G) + * \param[in] J_in_time time dependency of J (currently supported: constant, linear) + * \param[in] rho_in_time time dependency of rho (currently supported: constant, linear) + */ + PsatdAlgorithmFirstOrder ( + const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const SpectralFieldIndex& spectral_index, + const int norder_x, + const int norder_y, + const int norder_z, + const bool nodal, + const amrex::Real dt, + const bool div_cleaning, + const int J_in_time, + const int rho_in_time); + + /** + * \brief Updates E, B, F, and G fields in spectral space, + * according to the first-order PSATD equations + * + * \param[in,out] f all the fields in spectral space + */ + virtual void pushSpectralFields (SpectralFieldData& f) const override final; + + /** + * \brief Virtual function for current correction in Fourier space + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This function overrides the virtual function \c CurrentCorrection in the + * base class \c SpectralBaseAlgorithm and cannot be overridden by further + * derived classes. + * + * \param[in,out] field_data All fields in Fourier space + */ + virtual void CurrentCorrection (SpectralFieldData& field_data) override final; + + /** + * \brief Virtual function for Vay current deposition in Fourier space + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This function overrides the virtual function \c VayDeposition in the + * base class \c SpectralBaseAlgorithm and cannot be overridden by further + * derived classes. + * + * \param[in,out] field_data All fields in Fourier space + */ + virtual void VayDeposition (SpectralFieldData& field_data) override final; + + private: + + SpectralFieldIndex m_spectral_index; + + // Other member variables + amrex::Real m_dt; + bool m_div_cleaning; + int m_J_in_time; + int m_rho_in_time; +}; +#endif // WARPX_USE_PSATD +#endif // WARPX_PSATD_ALGORITHM_FIRST_ORDER_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp new file mode 100644 index 000000000..d32604760 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp @@ -0,0 +1,375 @@ +/* Copyright 2019 + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "PsatdAlgorithmFirstOrder.H" + +#include "Utils/TextMsg.H" +#include "Utils/WarpXAlgorithmSelection.H" +#include "Utils/WarpXConst.H" +#include "Utils/WarpX_Complex.H" + +#include <AMReX_Array4.H> +#include <AMReX_BLProfiler.H> +#include <AMReX_BaseFab.H> +#include <AMReX_BoxArray.H> +#include <AMReX_GpuComplex.H> +#include <AMReX_GpuLaunch.H> +#include <AMReX_GpuQualifiers.H> +#include <AMReX_IntVect.H> +#include <AMReX_MFIter.H> +#include <AMReX_PODVector.H> + +#include <cmath> + +#if WARPX_USE_PSATD + +using namespace amrex::literals; + +PsatdAlgorithmFirstOrder::PsatdAlgorithmFirstOrder( + const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const SpectralFieldIndex& spectral_index, + const int norder_x, + const int norder_y, + const int norder_z, + const bool nodal, + const amrex::Real dt, + const bool div_cleaning, + const int J_in_time, + const int rho_in_time) + // Initializer list + : SpectralBaseAlgorithm(spectral_kspace, dm, spectral_index, norder_x, norder_y, norder_z, nodal), + m_spectral_index(spectral_index), + m_dt(dt), + m_div_cleaning(div_cleaning), + m_J_in_time(J_in_time), + m_rho_in_time(rho_in_time) +{} + +void +PsatdAlgorithmFirstOrder::pushSpectralFields (SpectralFieldData& f) const +{ + const bool div_cleaning = m_div_cleaning; + + const bool J_constant = (m_J_in_time == JInTime::Constant) ? true : false; + const bool J_linear = (m_J_in_time == JInTime::Linear ) ? true : false; + const bool rho_constant = (m_rho_in_time == RhoInTime::Constant) ? true : false; + const bool rho_linear = (m_rho_in_time == RhoInTime::Linear ) ? true : false; + + const amrex::Real dt = m_dt; + const amrex::Real dt2 = dt*dt; + + const SpectralFieldIndex& Idx = m_spectral_index; + + // Loop over boxes + for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi) + { + const amrex::Box& bx = f.fields[mfi].box(); + + // Extract arrays for the fields to be updated + amrex::Array4<Complex> fields = f.fields[mfi].array(); + + // Extract pointers for the k vectors + const amrex::Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr(); +#if defined(WARPX_DIM_3D) + const amrex::Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr(); +#endif + const amrex::Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr(); + + // Loop over indices within one box + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Record old values of the fields to be updated + const Complex Ex_old = fields(i,j,k,Idx.Ex); + const Complex Ey_old = fields(i,j,k,Idx.Ey); + const Complex Ez_old = fields(i,j,k,Idx.Ez); + const Complex Bx_old = fields(i,j,k,Idx.Bx); + const Complex By_old = fields(i,j,k,Idx.By); + const Complex Bz_old = fields(i,j,k,Idx.Bz); + + // Shortcuts for the values of J and rho + const Complex Jx_mid = (J_constant) ? fields(i,j,k,Idx.Jx_mid) : 0._rt; + const Complex Jy_mid = (J_constant) ? fields(i,j,k,Idx.Jy_mid) : 0._rt; + const Complex Jz_mid = (J_constant) ? fields(i,j,k,Idx.Jz_mid) : 0._rt; + const Complex Jx_old = (J_linear ) ? fields(i,j,k,Idx.Jx_old) : 0._rt; + const Complex Jy_old = (J_linear ) ? fields(i,j,k,Idx.Jy_old) : 0._rt; + const Complex Jz_old = (J_linear ) ? fields(i,j,k,Idx.Jz_old) : 0._rt; + const Complex Jx_new = (J_linear ) ? fields(i,j,k,Idx.Jx_new) : 0._rt; + const Complex Jy_new = (J_linear ) ? fields(i,j,k,Idx.Jy_new) : 0._rt; + const Complex Jz_new = (J_linear ) ? fields(i,j,k,Idx.Jz_new) : 0._rt; + + const Complex Jx_c0 = (J_constant) ? Jx_mid : Jx_old; + const Complex Jy_c0 = (J_constant) ? Jy_mid : Jy_old; + const Complex Jz_c0 = (J_constant) ? Jz_mid : Jz_old; + const Complex Jx_c1 = (J_linear ) ? (Jx_new-Jx_old)/dt : 0._rt; + const Complex Jy_c1 = (J_linear ) ? (Jy_new-Jy_old)/dt : 0._rt; + const Complex Jz_c1 = (J_linear ) ? (Jz_new-Jz_old)/dt : 0._rt; + + Complex rho_mid, rho_old, rho_new, F_old, G_old; + Complex rho_c0, rho_c1; + if (div_cleaning) + { + rho_mid = (rho_constant) ? fields(i,j,k,Idx.rho_mid) : 0._rt; + rho_old = (rho_linear ) ? fields(i,j,k,Idx.rho_old) : 0._rt; + rho_new = (rho_linear ) ? fields(i,j,k,Idx.rho_new) : 0._rt; + + F_old = fields(i,j,k,Idx.F); + G_old = fields(i,j,k,Idx.G); + + rho_c0 = (rho_constant) ? rho_mid : rho_old; + rho_c1 = (rho_linear ) ? (rho_new-rho_old)/dt : 0._rt; + } + + // k vector values + const amrex::Real kx = modified_kx_arr[i]; +#if defined(WARPX_DIM_3D) + const amrex::Real ky = modified_ky_arr[j]; + const amrex::Real kz = modified_kz_arr[k]; +#else + constexpr amrex::Real ky = 0._rt; + const amrex::Real kz = modified_kz_arr[j]; +#endif + // Physical constants and imaginary unit + constexpr amrex::Real c = PhysConst::c; + constexpr amrex::Real c2 = c*c; + constexpr amrex::Real inv_c = 1._rt/c; + constexpr amrex::Real mu0 = PhysConst::mu0; + constexpr Complex I = Complex{0._rt, 1._rt}; + + const amrex::Real kx2 = kx*kx; + const amrex::Real ky2 = ky*ky; + const amrex::Real kz2 = kz*kz; + + const amrex::Real knorm = std::sqrt(kx2 + ky2 + kz2); + const amrex::Real knorm2 = knorm*knorm; + const amrex::Real knorm4 = knorm2*knorm2; + + // Auxiliary variables + const amrex::Real inv_knorm = 1._rt/knorm; + const amrex::Real inv_knorm2 = 1._rt/knorm2; + const amrex::Real inv_knorm4 = 1._rt/knorm4; + + const amrex::Real C = std::cos(c*knorm*dt); + const amrex::Real S = std::sin(c*knorm*dt); + + // Update equations + + if (knorm == 0._rt) + { + fields(i,j,k,Idx.Ex) = Ex_old - mu0*c2*dt*Jx_c0 - 0.5_rt*mu0*c2*dt2*Jx_c1; + fields(i,j,k,Idx.Ey) = Ey_old - mu0*c2*dt*Jy_c0 - 0.5_rt*mu0*c2*dt2*Jy_c1; + fields(i,j,k,Idx.Ez) = Ez_old - mu0*c2*dt*Jz_c0 - 0.5_rt*mu0*c2*dt2*Jz_c1; + + if (div_cleaning) + { + fields(i,j,k,Idx.F) = F_old - mu0*c2*dt*rho_c0 - 0.5_rt*mu0*c2*dt2*rho_c1; + } + } + else // knorm != 0 + { + Complex C01, C02, C03, C04, C05, C06, C07, C08, + C09, C10, C11, C12, C13, C14, C15, C16; + + // Ex + C01 = (div_cleaning) ? C : (kx2+ky2*C+kz2*C)*inv_knorm2; + C02 = (div_cleaning) ? 0._rt : kx*ky*(1._rt-C)*inv_knorm2; + C03 = (div_cleaning) ? 0._rt : kx*kz*(1._rt-C)*inv_knorm2; + C04 = 0._rt; + C05 = -I*c*kz*S*inv_knorm; + C06 = I*c*ky*S*inv_knorm; + C07 = (div_cleaning) ? I*c*kx*S*inv_knorm : 0._rt; + C09 = (div_cleaning) ? -mu0*c*S*inv_knorm : -mu0*c*(dt*c*kx2*knorm2+ky2*knorm*S+kz2*knorm*S)*inv_knorm4; + C10 = (div_cleaning) ? 0._rt : mu0*c*kx*ky*(knorm*S-dt*c*knorm2)*inv_knorm4; + C11 = (div_cleaning) ? 0._rt : mu0*c*kx*kz*(knorm*S-dt*c*knorm2)*inv_knorm4; + C12 = 0._rt; // This is not redundant, do not remove this + if (J_linear) C12 = (div_cleaning) ? mu0*(C-1._rt)*inv_knorm2 : mu0*(2._rt*ky2*(C-1._rt)+2._rt*kz2*(C-1._rt)-dt2*c2*kx2*knorm2)*inv_knorm4*0.5_rt; + C13 = (J_linear && !div_cleaning) ? mu0*kx*ky*(2._rt*(1._rt-C)-dt2*c2*knorm2)*inv_knorm4*0.5_rt : 0._rt; + C14 = (J_linear && !div_cleaning) ? mu0*kx*kz*(2._rt*(1._rt-C)-dt2*c2*knorm2)*inv_knorm4*0.5_rt : 0._rt; + C15 = (div_cleaning) ? I*mu0*c2*kx*(C-1._rt)*inv_knorm2 : 0._rt; + C16 = (div_cleaning && rho_linear) ? I*mu0*c*kx*(knorm*S-dt*c*knorm2)*inv_knorm4 : 0._rt; + + fields(i,j,k,Idx.Ex) = C01*Ex_old + C02*Ey_old + C03*Ez_old + + C04*Bx_old + C05*By_old + C06*Bz_old + + C07*F_old // only with div cleaning + + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0 + + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1 // only with J linear in time + + C15*rho_c0 // only with div cleaning + + C16*rho_c1; // only with div cleaning and rho linear in time + + // Ey + C01 = (div_cleaning) ? 0._rt : kx*ky*(1._rt-C)*inv_knorm2; + C02 = (div_cleaning) ? C : (kx2*C+ky2+kz2*C)*inv_knorm2; + C03 = (div_cleaning) ? 0._rt : ky*kz*(1._rt-C)*inv_knorm2; + C04 = I*c*kz*S*inv_knorm; + C05 = 0._rt; + C06 = -I*c*kx*S*inv_knorm; + C07 = (div_cleaning) ? I*c*ky*S*inv_knorm : 0._rt; + C09 = (div_cleaning) ? 0._rt : mu0*c*kx*ky*(knorm*S-dt*c*knorm2)*inv_knorm4; + C10 = (div_cleaning) ? -mu0*c*S*inv_knorm : -mu0*c*(dt*c*ky2*knorm2+kx2*knorm*S+kz2*knorm*S)*inv_knorm4; + C11 = (div_cleaning) ? 0._rt : mu0*c*ky*kz*(knorm*S-dt*c*knorm2)*inv_knorm4; + C12 = (J_linear && !div_cleaning) ? mu0*kx*ky*(2._rt*(1._rt-C)-dt2*c2*knorm2)*inv_knorm4*0.5_rt : 0._rt; + C13 = 0._rt; // This is not redundant, do not remove this + if (J_linear) C13 = (div_cleaning) ? mu0*(C-1._rt)*inv_knorm2 : mu0*(2._rt*kx2*(C-1._rt)+2._rt*kz2*(C-1._rt)-dt2*c2*ky2*knorm2)*inv_knorm4*0.5_rt; + C14 = (J_linear && !div_cleaning) ? mu0*ky*kz*(2._rt*(1._rt-C)-dt2*c2*knorm2)*inv_knorm4*0.5_rt : 0._rt; + C15 = (div_cleaning) ? I*mu0*c2*ky*(C-1._rt)*inv_knorm2 : 0._rt; + C16 = (div_cleaning && rho_linear) ? I*mu0*c*ky*(knorm*S-dt*c*knorm2)*inv_knorm4 : 0._rt; + + fields(i,j,k,Idx.Ey) = C01*Ex_old + C02*Ey_old + C03*Ez_old + + C04*Bx_old + C05*By_old + C06*Bz_old + + C07*F_old // only with div cleaning + + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0 + + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1 // only with J linear in time + + C15*rho_c0 // only with div cleaning + + C16*rho_c1; // only with div cleaning and rho linear in time + + // Ez + C01 = (div_cleaning) ? 0._rt : kx*kz*(1._rt-C)*inv_knorm2; + C02 = (div_cleaning) ? 0._rt : ky*kz*(1._rt-C)*inv_knorm2; + C03 = (div_cleaning) ? C : (kx2*C+ky2*C+kz2)*inv_knorm2; + C04 = -I*c*ky*S*inv_knorm; + C05 = I*c*kx*S*inv_knorm; + C06 = 0._rt; + C07 = (div_cleaning) ? I*c*kz*S*inv_knorm : 0._rt; + C09 = (div_cleaning) ? 0._rt : mu0*c*kx*kz*(knorm*S-dt*c*knorm2)*inv_knorm4; + C10 = (div_cleaning) ? 0._rt : mu0*c*ky*kz*(knorm*S-dt*c*knorm2)*inv_knorm4; + C11 = (div_cleaning) ? -mu0*c*S*inv_knorm : -mu0*c*(dt*c*kz2*knorm2+kx2*knorm*S+ky2*knorm*S)*inv_knorm4; + C12 = (J_linear && !div_cleaning) ? mu0*kx*kz*(2._rt*(1._rt-C)-dt2*c2*knorm2)*inv_knorm4*0.5_rt : 0._rt; + C13 = (J_linear && !div_cleaning) ? mu0*ky*kz*(2._rt*(1._rt-C)-dt2*c2*knorm2)*inv_knorm4*0.5_rt : 0._rt; + C14 = 0._rt; // This is not redundant, do not remove this + if (J_linear) C14 = (div_cleaning) ? mu0*(C-1._rt)*inv_knorm2 : mu0*(2._rt*kx2*(C-1._rt)+2._rt*ky2*(C-1._rt)-dt2*c2*kz2*knorm2)*inv_knorm4*0.5_rt; + C15 = (div_cleaning) ? I*mu0*c2*kz*(C-1._rt)*inv_knorm2 : 0._rt; + C16 = (div_cleaning && rho_linear) ? I*mu0*c*kz*(knorm*S-dt*c*knorm2)*inv_knorm4 : 0._rt; + + fields(i,j,k,Idx.Ez) = C01*Ex_old + C02*Ey_old + C03*Ez_old + + C04*Bx_old + C05*By_old + C06*Bz_old + + C07*F_old // only with div cleaning + + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0 + + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1 // only with J linear in time + + C15*rho_c0 // only with div cleaning + + C16*rho_c1; // only with div cleaning and rho linear in time + + // Bx + C01 = 0._rt; + C02 = I*kz*S*inv_knorm*inv_c; + C03 = -I*ky*S*inv_knorm*inv_c; + C04 = (div_cleaning) ? C : (kx2+ky2*C+kz2*C)*inv_knorm2; + C05 = (div_cleaning) ? 0._rt : kx*ky*(1._rt-C)*inv_knorm2; + C06 = (div_cleaning) ? 0._rt : kx*kz*(1._rt-C)*inv_knorm2; + C08 = (div_cleaning) ? I*kx*S*inv_knorm*inv_c : 0._rt; + C09 = 0._rt; + C10 = I*mu0*kz*(C-1._rt)*inv_knorm2; + C11 = -I*mu0*ky*(C-1._rt)*inv_knorm2; + C12 = 0._rt; + C13 = (J_linear) ? I*mu0*kz*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt; + C14 = (J_linear) ? -I*mu0*ky*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt; + + fields(i,j,k,Idx.Bx) = C01*Ex_old + C02*Ey_old + C03*Ez_old + + C04*Bx_old + C05*By_old + C06*Bz_old + + C08*G_old // only with div cleaning + + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0 + + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1; // only with J linear in time + + // By + C01 = -I*kz*S*inv_knorm*inv_c; + C02 = 0._rt; + C03 = I*kx*S*inv_knorm*inv_c; + C04 = (div_cleaning) ? 0._rt : kx*ky*(1._rt-C)*inv_knorm2; + C05 = (div_cleaning) ? C : (kx2*C+ky2+kz2*C)*inv_knorm2; + C06 = (div_cleaning) ? 0._rt : ky*kz*(1._rt-C)*inv_knorm2; + C08 = (div_cleaning) ? I*ky*S*inv_knorm*inv_c : 0._rt; + C09 = -I*mu0*kz*(C-1._rt)*inv_knorm2; + C10 = 0._rt; + C11 = I*mu0*kx*(C-1._rt)*inv_knorm2; + C12 = (J_linear) ? -I*mu0*kz*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt; + C13 = 0._rt; + C14 = (J_linear) ? I*mu0*kx*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt; + + fields(i,j,k,Idx.By) = C01*Ex_old + C02*Ey_old + C03*Ez_old + + C04*Bx_old + C05*By_old + C06*Bz_old + + C08*G_old // only with div cleaning + + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0 + + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1; // only with J linear in time + + // Bz + C01 = I*ky*S*inv_knorm*inv_c; + C02 = -I*kx*S*inv_knorm*inv_c; + C03 = 0._rt; + C04 = (div_cleaning) ? 0._rt : kx*kz*(1._rt-C)*inv_knorm2; + C05 = (div_cleaning) ? 0._rt : ky*kz*(1._rt-C)*inv_knorm2; + C06 = (div_cleaning) ? C : (kx2*C+ky2*C+kz2)*inv_knorm2; + C08 = (div_cleaning) ? I*kz*S*inv_knorm*inv_c : 0._rt; + C09 = I*mu0*ky*(C-1._rt)*inv_knorm2; + C10 = -I*mu0*kx*(C-1._rt)*inv_knorm2; + C11 = 0._rt; + C12 = (J_linear) ? I*mu0*ky*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt; + C13 = (J_linear) ? -I*mu0*kx*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt; + C14 = 0._rt; + + fields(i,j,k,Idx.Bz) = C01*Ex_old + C02*Ey_old + C03*Ez_old + + C04*Bx_old + C05*By_old + C06*Bz_old + + C08*G_old // only with div cleaning + + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0 + + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1; // only with J linear in time + + if (div_cleaning) + { + // F + C01 = I*kx*S*inv_knorm*inv_c; + C02 = I*ky*S*inv_knorm*inv_c; + C03 = I*kz*S*inv_knorm*inv_c; + C07 = C; + C09 = I*mu0*kx*(C-1._rt)*inv_knorm2; + C10 = I*mu0*ky*(C-1._rt)*inv_knorm2; + C11 = I*mu0*kz*(C-1._rt)*inv_knorm2; + C12 = (J_linear) ? I*mu0*kx*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt; + C13 = (J_linear) ? I*mu0*ky*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt; + C14 = (J_linear) ? I*mu0*kz*(knorm*S-dt*c*knorm2)*inv_knorm4*inv_c : 0._rt; + C15 = -mu0*c*S*inv_knorm; + C16 = (rho_linear) ? mu0*(C-1._rt)*inv_knorm2 : 0._rt; + + fields(i,j,k,Idx.F) = C01*Ex_old + C02*Ey_old + C03*Ez_old + + C07*F_old + + C09*Jx_c0 + C10*Jy_c0 + C11*Jz_c0 + + C12*Jx_c1 + C13*Jy_c1 + C14*Jz_c1 // only with J linear in time + + C15*rho_c0 + + C16*rho_c1; // only with rho linear in time + + // G + C04 = I*c*kx*S*inv_knorm; + C05 = I*c*ky*S*inv_knorm; + C06 = I*c*kz*S*inv_knorm; + C08 = C; + + fields(i,j,k,Idx.G) = C04*Bx_old + C05*By_old + C06*Bz_old + + C08*G_old; + } + } + }); + } +} + +void PsatdAlgorithmFirstOrder::CurrentCorrection (SpectralFieldData& field_data) +{ + // Profiling + BL_PROFILE("PsatdAlgorithmFirstOrder::CurrentCorrection"); + + amrex::ignore_unused(field_data); + amrex::Abort(Utils::TextMsg::Err( + "Current correction not implemented for first-order PSATD equations")); +} + +void +PsatdAlgorithmFirstOrder::VayDeposition (SpectralFieldData& field_data) +{ + // Profiling + BL_PROFILE("PsatdAlgorithmFirstOrder::VayDeposition()"); + + amrex::ignore_unused(field_data); + amrex::Abort(Utils::TextMsg::Err( + "Vay deposition not implemented for first-order PSATD equations")); +} + +#endif // WARPX_USE_PSATD diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp index be850e252..af3e2e336 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmGalileanRZ.cpp @@ -103,9 +103,9 @@ PsatdAlgorithmGalileanRZ::pushSpectralFields (SpectralFieldDataRZ & f) auto const Bp_m = Idx.Bx + Idx.n_fields*mode; auto const Bm_m = Idx.By + Idx.n_fields*mode; auto const Bz_m = Idx.Bz + Idx.n_fields*mode; - auto const Jp_m = Idx.Jx + Idx.n_fields*mode; - auto const Jm_m = Idx.Jy + Idx.n_fields*mode; - auto const Jz_m = Idx.Jz + Idx.n_fields*mode; + auto const Jp_m = Idx.Jx_mid + Idx.n_fields*mode; + auto const Jm_m = Idx.Jy_mid + Idx.n_fields*mode; + auto const Jz_m = Idx.Jz_mid + Idx.n_fields*mode; auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode; auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode; @@ -323,9 +323,9 @@ PsatdAlgorithmGalileanRZ::CurrentCorrection (SpectralFieldDataRZ& field_data) [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept { // All of the fields of each mode are grouped together - auto const Jp_m = Idx.Jx + Idx.n_fields*mode; - auto const Jm_m = Idx.Jy + Idx.n_fields*mode; - auto const Jz_m = Idx.Jz + Idx.n_fields*mode; + auto const Jp_m = Idx.Jx_mid + Idx.n_fields*mode; + auto const Jm_m = Idx.Jy_mid + Idx.n_fields*mode; + auto const Jz_m = Idx.Jz_mid + Idx.n_fields*mode; auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode; auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp index 9d1fbf3e1..c52437bf8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp @@ -190,9 +190,9 @@ PsatdAlgorithmJConstantInTime::pushSpectralFields (SpectralFieldData& f) const const Complex Bz_old = fields(i,j,k,Idx.Bz); // Shortcuts for the values of J - const Complex Jx = fields(i,j,k,Idx.Jx); - const Complex Jy = fields(i,j,k,Idx.Jy); - const Complex Jz = fields(i,j,k,Idx.Jz); + const Complex Jx = fields(i,j,k,Idx.Jx_mid); + const Complex Jy = fields(i,j,k,Idx.Jy_mid); + const Complex Jz = fields(i,j,k,Idx.Jz_mid); Complex F_old; if (dive_cleaning) @@ -751,9 +751,9 @@ void PsatdAlgorithmJConstantInTime::CurrentCorrection (SpectralFieldData& field_ ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Shortcuts for the values of J and rho - const Complex Jx = fields(i,j,k,Idx.Jx); - const Complex Jy = fields(i,j,k,Idx.Jy); - const Complex Jz = fields(i,j,k,Idx.Jz); + const Complex Jx = fields(i,j,k,Idx.Jx_mid); + const Complex Jy = fields(i,j,k,Idx.Jy_mid); + const Complex Jz = fields(i,j,k,Idx.Jz_mid); const Complex rho_old = fields(i,j,k,Idx.rho_old); const Complex rho_new = fields(i,j,k,Idx.rho_new); @@ -787,25 +787,25 @@ void PsatdAlgorithmJConstantInTime::CurrentCorrection (SpectralFieldData& field_ const Complex rho_old_mod = rho_old * amrex::exp(I * k_dot_vg * dt); const Complex den = 1._rt - amrex::exp(I * k_dot_vg * dt); - fields(i,j,k,Idx.Jx) = Jx - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) + fields(i,j,k,Idx.Jx_mid) = Jx - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) * kx / (k_norm * k_norm); - fields(i,j,k,Idx.Jy) = Jy - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) + fields(i,j,k,Idx.Jy_mid) = Jy - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) * ky / (k_norm * k_norm); - fields(i,j,k,Idx.Jz) = Jz - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) + fields(i,j,k,Idx.Jz_mid) = Jz - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) * kz / (k_norm * k_norm); } else { - fields(i,j,k,Idx.Jx) = Jx - (k_dot_J - I * (rho_new - rho_old) / dt) + fields(i,j,k,Idx.Jx_mid) = Jx - (k_dot_J - I * (rho_new - rho_old) / dt) * kx / (k_norm * k_norm); - fields(i,j,k,Idx.Jy) = Jy - (k_dot_J - I * (rho_new - rho_old) / dt) + fields(i,j,k,Idx.Jy_mid) = Jy - (k_dot_J - I * (rho_new - rho_old) / dt) * ky / (k_norm * k_norm); - fields(i,j,k,Idx.Jz) = Jz - (k_dot_J - I * (rho_new - rho_old) / dt) + fields(i,j,k,Idx.Jz_mid) = Jz - (k_dot_J - I * (rho_new - rho_old) / dt) * kz / (k_norm * k_norm); } } @@ -840,11 +840,11 @@ PsatdAlgorithmJConstantInTime::VayDeposition (SpectralFieldData& field_data) ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Shortcuts for the values of D - const Complex Dx = fields(i,j,k,Idx.Jx); + const Complex Dx = fields(i,j,k,Idx.Jx_mid); #if defined(WARPX_DIM_3D) - const Complex Dy = fields(i,j,k,Idx.Jy); + const Complex Dy = fields(i,j,k,Idx.Jy_mid); #endif - const Complex Dz = fields(i,j,k,Idx.Jz); + const Complex Dz = fields(i,j,k,Idx.Jz_mid); // Imaginary unit constexpr Complex I = Complex{0._rt, 1._rt}; @@ -859,18 +859,18 @@ PsatdAlgorithmJConstantInTime::VayDeposition (SpectralFieldData& field_data) #endif // Compute Jx - if (kx_mod != 0._rt) fields(i,j,k,Idx.Jx) = I * Dx / kx_mod; - else fields(i,j,k,Idx.Jx) = 0._rt; + if (kx_mod != 0._rt) fields(i,j,k,Idx.Jx_mid) = I * Dx / kx_mod; + else fields(i,j,k,Idx.Jx_mid) = 0._rt; #if defined(WARPX_DIM_3D) // Compute Jy - if (ky_mod != 0._rt) fields(i,j,k,Idx.Jy) = I * Dy / ky_mod; - else fields(i,j,k,Idx.Jy) = 0._rt; + if (ky_mod != 0._rt) fields(i,j,k,Idx.Jy_mid) = I * Dy / ky_mod; + else fields(i,j,k,Idx.Jy_mid) = 0._rt; #endif // Compute Jz - if (kz_mod != 0._rt) fields(i,j,k,Idx.Jz) = I * Dz / kz_mod; - else fields(i,j,k,Idx.Jz) = 0._rt; + if (kz_mod != 0._rt) fields(i,j,k,Idx.Jz_mid) = I * Dz / kz_mod; + else fields(i,j,k,Idx.Jz_mid) = 0._rt; }); } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp index bd9df977e..861001c78 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp @@ -120,9 +120,9 @@ PsatdAlgorithmJLinearInTime::pushSpectralFields (SpectralFieldData& f) const const Complex Bz_old = fields(i,j,k,Idx.Bz); // Shortcuts for the values of J and rho - const Complex Jx_old = fields(i,j,k,Idx.Jx); - const Complex Jy_old = fields(i,j,k,Idx.Jy); - const Complex Jz_old = fields(i,j,k,Idx.Jz); + const Complex Jx_old = fields(i,j,k,Idx.Jx_old); + const Complex Jy_old = fields(i,j,k,Idx.Jy_old); + const Complex Jz_old = fields(i,j,k,Idx.Jz_old); const Complex Jx_new = fields(i,j,k,Idx.Jx_new); const Complex Jy_new = fields(i,j,k,Idx.Jy_new); const Complex Jz_new = fields(i,j,k,Idx.Jz_new); diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp index 55b58821c..effb1cc2b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp @@ -83,7 +83,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) const bool update_with_rho = m_update_with_rho; const bool time_averaging = m_time_averaging; - const bool J_in_time_linear = (m_J_in_time == JInTime::Linear) ? true : false; + const bool J_linear = (m_J_in_time == JInTime::Linear) ? true : false; const bool dive_cleaning = m_dive_cleaning; const bool divb_cleaning = m_divb_cleaning; @@ -112,7 +112,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) amrex::Array4<const amrex::Real> X5_arr; amrex::Array4<const amrex::Real> X6_arr; - if (time_averaging && J_in_time_linear) + if (time_averaging && J_linear) { X5_arr = X5_coef[mfi].array(); X6_arr = X6_coef[mfi].array(); @@ -131,6 +131,9 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) amrex::ParallelFor(bx, modes, [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept { + int idx_jx = (J_linear) ? static_cast<int>(Idx.Jx_old) : static_cast<int>(Idx.Jx_mid); + int idx_jy = (J_linear) ? static_cast<int>(Idx.Jy_old) : static_cast<int>(Idx.Jy_mid); + int idx_jz = (J_linear) ? static_cast<int>(Idx.Jz_old) : static_cast<int>(Idx.Jz_mid); // All of the fields of each mode are grouped together int const Ep_m = Idx.Ex + Idx.n_fields*mode; @@ -139,9 +142,9 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) int const Bp_m = Idx.Bx + Idx.n_fields*mode; int const Bm_m = Idx.By + Idx.n_fields*mode; int const Bz_m = Idx.Bz + Idx.n_fields*mode; - int const Jp_m = Idx.Jx + Idx.n_fields*mode; - int const Jm_m = Idx.Jy + Idx.n_fields*mode; - int const Jz_m = Idx.Jz + Idx.n_fields*mode; + int const Jp_m = idx_jx + Idx.n_fields*mode; + int const Jm_m = idx_jy + Idx.n_fields*mode; + int const Jz_m = idx_jz + Idx.n_fields*mode; int const rho_old_m = Idx.rho_old + Idx.n_fields*mode; int const rho_new_m = Idx.rho_new + Idx.n_fields*mode; @@ -238,7 +241,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) G_old = fields(i,j,k,G_m); } - if (J_in_time_linear) + if (J_linear) { const int Jp_m_new = Idx.Jx_new + Idx.n_fields*mode; const int Jm_m_new = Idx.Jy_new + Idx.n_fields*mode; @@ -335,7 +338,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const & f) { const bool time_averaging = m_time_averaging; - const bool J_in_time_linear = (m_J_in_time == JInTime::Linear) ? true : false; + const bool J_linear = (m_J_in_time == JInTime::Linear) ? true : false; // Fill them with the right values: // Loop over boxes and allocate the corresponding coefficients @@ -356,7 +359,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const amrex::Array4<amrex::Real> X5; amrex::Array4<amrex::Real> X6; - if (time_averaging && J_in_time_linear) + if (time_averaging && J_linear) { X5 = X5_coef[mfi].array(); X6 = X6_coef[mfi].array(); @@ -395,7 +398,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const X3(i,j,k,mode) = - c*c * dt*dt / (3._rt*ep0); } - if (time_averaging && J_in_time_linear) + if (time_averaging && J_linear) { constexpr amrex::Real c2 = PhysConst::c; const amrex::Real dt3 = dt * dt * dt; @@ -450,9 +453,9 @@ PsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data) [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept { // All of the fields of each mode are grouped together - auto const Jp_m = Idx.Jx + Idx.n_fields*mode; - auto const Jm_m = Idx.Jy + Idx.n_fields*mode; - auto const Jz_m = Idx.Jz + Idx.n_fields*mode; + auto const Jp_m = Idx.Jx_mid + Idx.n_fields*mode; + auto const Jm_m = Idx.Jy_mid + Idx.n_fields*mode; + auto const Jz_m = Idx.Jz_mid + Idx.n_fields*mode; auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode; auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 4ab88f1a3..c7848d731 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -83,18 +83,19 @@ class SpectralFieldIndex // Always int Ex = -1, Ey = -1, Ez = -1; int Bx = -1, By = -1, Bz = -1; - int Jx = -1, Jy = -1, Jz = -1; - int rho_old = -1, rho_new = -1, divE = -1; + int divE = -1; // Time averaging int Ex_avg = -1, Ey_avg = -1, Ez_avg = -1; int Bx_avg = -1, By_avg = -1, Bz_avg = -1; - // J linear in time + // J + int Jx_old = -1, Jy_old = -1, Jz_old = -1; + int Jx_mid = -1, Jy_mid = -1, Jz_mid = -1; int Jx_new = -1, Jy_new = -1, Jz_new = -1; - // rho quadratic in time - int rho_mid = -1; + // rho + int rho_old = -1, rho_mid = -1, rho_new = -1; // div(E) and div(B) cleaning int F = -1, G = -1; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 56189e0ff..e1db4e5ab 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -51,10 +51,8 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, { Ex = c++; Ey = c++; Ez = c++; Bx = c++; By = c++; Bz = c++; - Jx = c++; Jy = c++; Jz = c++; // TODO Allocate rho_old and rho_new only when needed - rho_old = c++; rho_new = c++; // Reuse data corresponding to index Bx = 3 to avoid storing extra memory divE = 3; @@ -69,17 +67,25 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, if (divb_cleaning) G = c++; - if (J_in_time == JInTime::Linear) + if (J_in_time == JInTime::Constant) { - Jx_new = c++; - Jy_new = c++; - Jz_new = c++; + Jx_mid = c++; Jy_mid = c++; Jz_mid = c++; + } + else if (J_in_time == JInTime::Linear) + { + Jx_old = c++; Jy_old = c++; Jz_old = c++; + Jx_new = c++; Jy_new = c++; Jz_new = c++; } - if (rho_in_time == RhoInTime::Quadratic) + if (rho_in_time == RhoInTime::Constant) { rho_mid = c++; } + else if (rho_in_time == RhoInTime::Linear) + { + rho_old = c++; + rho_new = c++; + } if (pml_rz) { diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 38b242010..da4b9687b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -58,6 +58,8 @@ class SpectralSolver * (no domain decomposition) * \param[in] update_with_rho whether rho is used in the field update equations * \param[in] fft_do_time_averaging whether the time averaging algorithm is used + * \param[in] psatd_solution_type whether the PSATD equations are derived + * from a first-order or second-order model * \param[in] J_in_time integer that corresponds to the time dependency of J * (constant, linear) for the PSATD algorithm * \param[in] rho_in_time integer that corresponds to the time dependency of rho @@ -80,6 +82,7 @@ class SpectralSolver const bool periodic_single_box, const bool update_with_rho, const bool fft_do_time_averaging, + const int psatd_solution_type, const int J_in_time, const int rho_in_time, const bool dive_cleaning, diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index c0d7b8941..f29a0bd03 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -8,10 +8,12 @@ #include "FieldSolver/SpectralSolver/SpectralFieldData.H" #include "SpectralAlgorithms/PsatdAlgorithmComoving.H" #include "SpectralAlgorithms/PsatdAlgorithmPml.H" +#include "SpectralAlgorithms/PsatdAlgorithmFirstOrder.H" #include "SpectralAlgorithms/PsatdAlgorithmJConstantInTime.H" #include "SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H" #include "SpectralKSpace.H" #include "SpectralSolver.H" +#include "Utils/TextMsg.H" #include "Utils/WarpXProfilerWrapper.H" #include <memory> @@ -30,6 +32,7 @@ SpectralSolver::SpectralSolver( const bool pml, const bool periodic_single_box, const bool update_with_rho, const bool fft_do_time_averaging, + const int psatd_solution_type, const int J_in_time, const int rho_in_time, const bool dive_cleaning, @@ -49,13 +52,13 @@ SpectralSolver::SpectralSolver( // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space - if (pml) // PSATD equations in the PML grids + if (pml) // PSATD equations in the PML region { algorithm = std::make_unique<PsatdAlgorithmPml>( k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, dt, dive_cleaning, divb_cleaning); } - else // PSATD equations in the regulard grids + else // PSATD equations in the regular domain { // Comoving PSATD algorithm if (v_comoving[0] != 0. || v_comoving[1] != 0. || v_comoving[2] != 0.) @@ -64,7 +67,31 @@ SpectralSolver::SpectralSolver( k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, v_comoving, dt, update_with_rho); } - else // PSATD algorithms: standard, Galilean, averaged Galilean, multi-J + // Galilean PSATD algorithm (only J constant in time) + else if (v_galilean[0] != 0. || v_galilean[1] != 0. || v_galilean[2] != 0.) + { + algorithm = std::make_unique<PsatdAlgorithmJConstantInTime>( + k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, + v_galilean, dt, update_with_rho, fft_do_time_averaging, + dive_cleaning, divb_cleaning); + } + else if (psatd_solution_type == PSATDSolutionType::FirstOrder) + { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + !fft_do_time_averaging, + "psatd.do_time_averaging=1 not supported when psatd.solution_type=first-order"); + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + (!dive_cleaning && !divb_cleaning) || (dive_cleaning && divb_cleaning), + "warpx.do_dive_cleaning and warpx.do_divb_cleaning must be equal when psatd.solution_type=first-order"); + + const bool div_cleaning = (dive_cleaning && divb_cleaning); + + algorithm = std::make_unique<PsatdAlgorithmFirstOrder>( + k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, + dt, div_cleaning, J_in_time, rho_in_time); + } + else if (psatd_solution_type == PSATDSolutionType::SecondOrder) { if (J_in_time == JInTime::Constant) { @@ -73,7 +100,7 @@ SpectralSolver::SpectralSolver( v_galilean, dt, update_with_rho, fft_do_time_averaging, dive_cleaning, divb_cleaning); } - else // J linear in time + else if (J_in_time == JInTime::Linear) { algorithm = std::make_unique<PsatdAlgorithmJLinearInTime>( k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 3d68e8e52..9df4bb21b 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -280,9 +280,9 @@ void WarpX::PSATDForwardTransformJ ( { Idx = spectral_solver_fp[lev]->m_spectral_index; - idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx); - idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy); - idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz); + idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx_mid); + idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy_mid); + idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz_mid); ForwardTransformVect(lev, *spectral_solver_fp[lev], J_fp[lev], idx_jx, idx_jy, idx_jz); @@ -290,9 +290,9 @@ void WarpX::PSATDForwardTransformJ ( { Idx = spectral_solver_cp[lev]->m_spectral_index; - idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx); - idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy); - idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz); + idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx_mid); + idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy_mid); + idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz_mid); ForwardTransformVect(lev, *spectral_solver_cp[lev], J_cp[lev], idx_jx, idx_jy, idx_jz); } @@ -304,11 +304,23 @@ void WarpX::PSATDForwardTransformJ ( { for (int lev = 0; lev <= finest_level; ++lev) { - spectral_solver_fp[lev]->ApplyFilter(lev, Idx.Jx, Idx.Jy, Idx.Jz); + Idx = spectral_solver_fp[lev]->m_spectral_index; + + idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx_mid); + idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy_mid); + idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz_mid); + + spectral_solver_fp[lev]->ApplyFilter(lev, idx_jx, idx_jy, idx_jz); if (spectral_solver_cp[lev]) { - spectral_solver_cp[lev]->ApplyFilter(lev, Idx.Jx, Idx.Jy, Idx.Jz); + Idx = spectral_solver_cp[lev]->m_spectral_index; + + idx_jx = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jx_new) : static_cast<int>(Idx.Jx_mid); + idx_jy = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jy_new) : static_cast<int>(Idx.Jy_mid); + idx_jz = (J_in_time == JInTime::Linear) ? static_cast<int>(Idx.Jz_new) : static_cast<int>(Idx.Jz_mid); + + spectral_solver_cp[lev]->ApplyFilter(lev, idx_jx, idx_jy, idx_jz); } } } @@ -328,9 +340,11 @@ void WarpX::PSATDBackwardTransformJ ( { Idx = spectral_solver_fp[lev]->m_spectral_index; - idx_jx = static_cast<int>(Idx.Jx); - idx_jy = static_cast<int>(Idx.Jy); - idx_jz = static_cast<int>(Idx.Jz); + // Note that these backward FFTs are currently called only + // with algorithms that do not support J linear in time + idx_jx = static_cast<int>(Idx.Jx_mid); + idx_jy = static_cast<int>(Idx.Jy_mid); + idx_jz = static_cast<int>(Idx.Jz_mid); BackwardTransformVect(lev, *spectral_solver_fp[lev], J_fp[lev], idx_jx, idx_jy, idx_jz, m_fill_guards_current); @@ -339,9 +353,11 @@ void WarpX::PSATDBackwardTransformJ ( { Idx = spectral_solver_cp[lev]->m_spectral_index; - idx_jx = static_cast<int>(Idx.Jx); - idx_jy = static_cast<int>(Idx.Jy); - idx_jz = static_cast<int>(Idx.Jz); + // Note that these backward FFTs are currently called only + // with algorithms that do not support J linear in time + idx_jx = static_cast<int>(Idx.Jx_mid); + idx_jy = static_cast<int>(Idx.Jy_mid); + idx_jz = static_cast<int>(Idx.Jz_mid); BackwardTransformVect(lev, *spectral_solver_cp[lev], J_cp[lev], idx_jx, idx_jy, idx_jz, m_fill_guards_current); @@ -359,7 +375,15 @@ void WarpX::PSATDForwardTransformRho ( const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; // Select index in k space - const int dst_comp = (dcomp == 0) ? Idx.rho_old : Idx.rho_new; + int dst_comp; + if (rho_in_time == RhoInTime::Constant) + { + dst_comp = Idx.rho_mid; + } + else // rho_in_time == RhoInTime::Linear + { + dst_comp = (dcomp == 0) ? Idx.rho_old : Idx.rho_new; + } for (int lev = 0; lev <= finest_level; ++lev) { @@ -568,15 +592,15 @@ WarpX::PSATDMoveJNewToJOld () for (int lev = 0; lev <= finest_level; ++lev) { - spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jx_new, Idx.Jx); - spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jy_new, Idx.Jy); - spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jz_new, Idx.Jz); + spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jx_new, Idx.Jx_old); + spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jy_new, Idx.Jy_old); + spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jz_new, Idx.Jz_old); if (spectral_solver_cp[lev]) { - spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jx_new, Idx.Jx); - spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jy_new, Idx.Jy); - spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jz_new, Idx.Jz); + spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jx_new, Idx.Jx_old); + spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jy_new, Idx.Jy_old); + spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jz_new, Idx.Jz_old); } } } diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 0d3e919a0..32cd11dfd 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -480,7 +480,7 @@ WarpX::InitPML () pml_ncell, pml_delta, amrex::IntVect::TheZeroVector(), dt[0], nox_fft, noy_fft, noz_fft, do_nodal, do_moving_window, pml_has_particles, do_pml_in_domain, - J_in_time, rho_in_time, + psatd_solution_type, J_in_time, rho_in_time, do_pml_dive_cleaning, do_pml_divb_cleaning, amrex::IntVect(0), amrex::IntVect(0), guard_cells.ng_FieldSolver.max(), @@ -517,7 +517,7 @@ WarpX::InitPML () pml_ncell, pml_delta, refRatio(lev-1), dt[lev], nox_fft, noy_fft, noz_fft, do_nodal, do_moving_window, pml_has_particles, do_pml_in_domain, - J_in_time, rho_in_time, do_pml_dive_cleaning, do_pml_divb_cleaning, + psatd_solution_type, J_in_time, rho_in_time, do_pml_dive_cleaning, do_pml_divb_cleaning, amrex::IntVect(0), amrex::IntVect(0), guard_cells.ng_FieldSolver.max(), v_particle_pml, diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index c3160cfad..936fd3d2b 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -84,6 +84,13 @@ struct GatheringAlgo { }; }; +struct PSATDSolutionType { + enum { + FirstOrder = 0, + SecondOrder = 1 + }; +}; + struct JInTime { enum { Constant = 0, @@ -93,8 +100,8 @@ struct JInTime { struct RhoInTime { enum { - Linear = 1, - Quadratic = 2 + Constant = 0, + Linear = 1 }; }; diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index b99459b46..c5ef16743 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -63,6 +63,12 @@ const std::map<std::string, int> gathering_algo_to_int = { {"default", GatheringAlgo::EnergyConserving } }; +const std::map<std::string, int> psatd_solution_type_to_int = { + {"first-order", PSATDSolutionType::FirstOrder}, + {"second-order", PSATDSolutionType::SecondOrder}, + {"default", PSATDSolutionType::SecondOrder} +}; + const std::map<std::string, int> J_in_time_to_int = { {"constant", JInTime::Constant}, {"linear", JInTime::Linear}, @@ -70,8 +76,8 @@ const std::map<std::string, int> J_in_time_to_int = { }; const std::map<std::string, int> rho_in_time_to_int = { + {"constant", RhoInTime::Constant}, {"linear", RhoInTime::Linear}, - {"quadratic", RhoInTime::Quadratic}, {"default", RhoInTime::Linear} }; @@ -145,6 +151,8 @@ GetAlgorithmInteger( amrex::ParmParse& pp, const char* pp_search_key ){ algo_to_int = charge_deposition_algo_to_int; } else if (0 == std::strcmp(pp_search_key, "field_gathering")) { algo_to_int = gathering_algo_to_int; + } else if (0 == std::strcmp(pp_search_key, "solution_type")) { + algo_to_int = psatd_solution_type_to_int; } else if (0 == std::strcmp(pp_search_key, "J_in_time")) { algo_to_int = J_in_time_to_int; } else if (0 == std::strcmp(pp_search_key, "rho_in_time")) { diff --git a/Source/WarpX.H b/Source/WarpX.H index 20d4e4210..2d1124c06 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -187,6 +187,11 @@ public: */ static amrex::Vector<ParticleBoundaryType> particle_boundary_hi; + //! Integer that corresponds to the order of the PSATD solution + //! (whether the PSATD equations are derived from first-order or + //! second-order solution) + static short psatd_solution_type; + //! Integers that correspond to the time dependency of J (constant, linear) //! and rho (linear, quadratic) for the PSATD algorithm static short J_in_time; @@ -1639,7 +1644,7 @@ private: const int icomp, const int dcomp, const bool apply_kspace_filter=true); /** - * \brief Copy rho_new to rho_old in spectral space + * \brief Copy rho_new to rho_old in spectral space (when rho is linear in time) */ void PSATDMoveRhoNewToRhoOld (); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 579387433..e51197089 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -120,6 +120,7 @@ short WarpX::charge_deposition_algo; short WarpX::field_gathering_algo; short WarpX::particle_pusher_algo; short WarpX::electromagnetic_solver_id; +short WarpX::psatd_solution_type; short WarpX::J_in_time; short WarpX::rho_in_time; short WarpX::load_balance_costs_update_algo; @@ -1148,6 +1149,11 @@ WarpX::ReadParameters () WARPX_ALWAYS_ASSERT_WITH_MESSAGE(noz_fft > 0, "PSATD order must be finite unless psatd.periodic_single_box_fft is used"); } + // Integer that corresponds to the order of the PSATD solution + // (whether the PSATD equations are derived from first-order or + // second-order solution) + psatd_solution_type = GetAlgorithmInteger(pp_psatd, "solution_type"); + // Integers that correspond to the time dependency of J (constant, linear) // and rho (linear, quadratic) for the PSATD algorithm J_in_time = GetAlgorithmInteger(pp_psatd, "J_in_time"); @@ -1311,13 +1317,6 @@ WarpX::ReadParameters () ); } - if (J_in_time == JInTime::Constant) - { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - rho_in_time == RhoInTime::Linear, - "psatd.J_in_time=constant supports only psatd.rho_in_time=linear"); - } - if (J_in_time == JInTime::Linear) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( @@ -1331,6 +1330,14 @@ WarpX::ReadParameters () WARPX_ALWAYS_ASSERT_WITH_MESSAGE( v_comoving_is_zero, "psatd.J_in_time=linear not implemented with comoving PSATD"); + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + !current_correction, + "psatd.current_correction=1 not implemented with psatd.J_in_time=linear"); + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + current_deposition_algo != CurrentDepositionAlgo::Vay, + "algo.current_deposition=vay not implemented with psatd.J_in_time=linear"); } for (int dir = 0; dir < AMREX_SPACEDIM; dir++) @@ -2334,6 +2341,7 @@ void WarpX::AllocLevelSpectralSolver (amrex::Vector<std::unique_ptr<SpectralSolv fft_periodic_single_box, update_with_rho, fft_do_time_averaging, + psatd_solution_type, J_in_time, rho_in_time, do_dive_cleaning, |