diff --git a/Examples/Tests/energy_conserving_thermal_plasma/analysis.py b/Examples/Tests/energy_conserving_thermal_plasma/analysis.py new file mode 100755 index 00000000000..43e9b6d9822 --- /dev/null +++ b/Examples/Tests/energy_conserving_thermal_plasma/analysis.py @@ -0,0 +1,38 @@ +#!/usr/bin/env python3 + +# Copyright 2024, Remi Lehe +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + +# This script tests the conservation of energy for a thermal plasma with periodic boundary. +# Here, we use energy-converving gather and an electrostatic solver. The energy +# is not expected to be exactly conserved, but it is expected to be better conserved +# than other gathering scheme. This tests checks that the energy does not increase by +# more than 0.3% over the duration of the simulatoin. + +import os +import sys + +import numpy as np + +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +# this will be the name of the plot file +fn = sys.argv[1] + +# Get energy as a function of time, from reduced diagnostics +EFdata = np.genfromtxt('./diags/reducedfiles/EF.txt') # Field energy +EPdata = np.genfromtxt('./diags/reducedfiles/EP.txt') # Particle energy +field_energy = EFdata[:,2] +particle_energy = EPdata[:,2] +E = field_energy + particle_energy +print(abs(E-E[0])/E[0]) +# Check that the energy is conserved to 0.3% +assert np.all( abs(E-E[0])/E[0] < 0.003 ) + +# Checksum test +test_name = os.path.split(os.getcwd())[1] +checksumAPI.evaluate_checksum(test_name, fn) diff --git a/Examples/Tests/energy_conserving_thermal_plasma/inputs_2d_electrostatic b/Examples/Tests/energy_conserving_thermal_plasma/inputs_2d_electrostatic new file mode 100644 index 00000000000..dab1b77db50 --- /dev/null +++ b/Examples/Tests/energy_conserving_thermal_plasma/inputs_2d_electrostatic @@ -0,0 +1,83 @@ +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 500 +amr.n_cell = 8 8 +amr.max_level = 0 +geometry.dims = 2 + +warpx.do_electrostatic = labframe +algo.field_gathering = energy-conserving +algo.particle_shape = 2 +warpx.use_filter = 0 + +################################# +############ CONSTANTS ############# +################################# +my_constants.n0 = 1.e30 # plasma density, m^-3 +my_constants.Ti = 100. # ion temperature, eV +my_constants.Te = 100. # electron temperature, eV +my_constants.wpe = q_e*sqrt(n0/(m_e*epsilon0)) # electron plasma frequency, radians/s +my_constants.de0 = clight/wpe # skin depth, m +my_constants.dt = ( 0.2 )/wpe # time step size, s +warpx.const_dt = dt + +################################# +####### GENERAL PARAMETERS ###### +################################# +geometry.dims = 2 +geometry.prob_lo = 0.0 0.0 +geometry.prob_hi = 10.*de0 10.*de0 +warpx.serialize_initial_conditions = 1 +warpx.verbose = 1 + +################################# +###### BOUNDARY CONDITIONS ###### +################################# +boundary.field_lo = periodic periodic +boundary.field_hi = periodic periodic +boundary.particle_lo = periodic periodic +boundary.particle_hi = periodic periodic + +################################# +############ PLASMA ############# +################################# +particles.species_names = electrons protons + +electrons.species_type = electron +electrons.injection_style = "NUniformPerCell" +electrons.num_particles_per_cell_each_dim = 2 2 +electrons.profile = constant +electrons.density = n0 +electrons.momentum_distribution_type = gaussian +electrons.ux_th = sqrt(Te*q_e/m_e)/clight +electrons.uy_th = sqrt(Te*q_e/m_e)/clight +electrons.uz_th = sqrt(Te*q_e/m_e)/clight +electrons.xmin = 0 +electrons.xmax = 10.*de0 +electrons.zmin = 0 +electrons.zmax =10.*de0 + +protons.species_type = proton +protons.injection_style = "NUniformPerCell" +protons.num_particles_per_cell_each_dim = 2 2 +protons.profile = constant +protons.density = n0 +protons.momentum_distribution_type = gaussian +protons.ux_th = sqrt(Ti*q_e/m_p)/clight +protons.uy_th = sqrt(Ti*q_e/m_p)/clight +protons.uz_th = sqrt(Ti*q_e/m_p)/clight +protons.xmin = 0 +protons.xmax = 10.*de0 +protons.zmin = 0 +protons.zmax =10.*de0 + +diagnostics.diags_names = diag1 +diag1.intervals = 500 +diag1.diag_type = Full + +warpx.reduced_diags_names = EP EF +EP.type = ParticleEnergy +EP.intervals = 100 +EF.type = FieldEnergy +EF.intervals =100 diff --git a/Regression/Checksum/benchmarks_json/BeamBeamCollision.json b/Regression/Checksum/benchmarks_json/BeamBeamCollision.json index 16c27055ab5..693faf43882 100644 --- a/Regression/Checksum/benchmarks_json/BeamBeamCollision.json +++ b/Regression/Checksum/benchmarks_json/BeamBeamCollision.json @@ -1,96 +1,96 @@ { "lev=0": { - "Bx": 971135657171.612, - "By": 971078812454.5405, - "Bz": 20140193.235893946, - "Ex": 2.9111756162943966e+20, - "Ey": 2.9113725115697712e+20, - "Ez": 1.0213536367191107e+17, - "rho_beam1": 7.970337929028706e+16, - "rho_beam2": 7.969213804851568e+16, - "rho_ele1": 343600677331163.7, - "rho_ele2": 302746939366837.25, - "rho_pos1": 333855946230626.06, - "rho_pos2": 310879461124837.44 + "Bx": 970573570349.2142, + "By": 970493693614.8823, + "Bz": 18743533.08373785, + "Ex": 2.9100804648817213e+20, + "Ey": 2.910188719414779e+20, + "Ez": 1.309938057448976e+17, + "rho_beam1": 7.969171294602906e+16, + "rho_beam2": 7.969079911431987e+16, + "rho_ele1": 234837475722889.97, + "rho_ele2": 277285508564124.12, + "rho_pos1": 229824415763789.62, + "rho_pos2": 286513388076434.4 }, "beam1": { - "particle_opticalDepthQSR": 104909.59461909423, - "particle_position_x": 0.001500222221634118, - "particle_position_y": 0.0015002445303035634, - "particle_position_z": 0.0049656251227976015, - "particle_momentum_x": 6.205341799808723e-15, - "particle_momentum_y": 6.1592257603817594e-15, - "particle_momentum_z": 6.806886719670214e-12, - "particle_weight": 635949610.5135971 + "particle_opticalDepthQSR": 104947.58821047074, + "particle_position_x": 0.0015001846430437182, + "particle_position_y": 0.0015002200838688795, + "particle_position_z": 0.004965556992831605, + "particle_momentum_x": 6.207861637949496e-15, + "particle_momentum_y": 6.16591876835476e-15, + "particle_momentum_z": 6.806755306915448e-12, + "particle_weight": 635859587.7915188 }, "beam2": { - "particle_opticalDepthQSR": 104164.848014815, - "particle_position_x": 0.0015001011957527532, - "particle_position_y": 0.001500139975740741, - "particle_position_z": 0.004965479176845744, - "particle_momentum_x": 6.200690794877584e-15, - "particle_momentum_y": 6.186048913459861e-15, - "particle_momentum_z": 6.7990490255176515e-12, - "particle_weight": 635863144.251134 + "particle_opticalDepthQSR": 104180.42997257302, + "particle_position_x": 0.0015001132145418016, + "particle_position_y": 0.001500162338878405, + "particle_position_z": 0.004965528365042041, + "particle_momentum_x": 6.201791193788119e-15, + "particle_momentum_y": 6.188946480983165e-15, + "particle_momentum_z": 6.796967564632488e-12, + "particle_weight": 635853166.1373858 }, "ele1": { - "particle_opticalDepthQSR": 435.73003117907257, - "particle_position_x": 4.882183530045367e-06, - "particle_position_y": 4.841391483672882e-06, - "particle_position_z": 1.8449175055560687e-05, - "particle_momentum_x": 5.6656608489971696e-18, - "particle_momentum_y": 5.724425295258085e-18, - "particle_momentum_z": 2.6277553331470036e-15, - "particle_weight": 2696555.9200674472 + "particle_opticalDepthQSR": 385.3959944370225, + "particle_position_x": 5.063867198487091e-06, + "particle_position_y": 5.323267522706671e-06, + "particle_position_z": 1.81974337459613e-05, + "particle_momentum_x": 5.220808698440275e-18, + "particle_momentum_y": 5.2833645967924376e-18, + "particle_momentum_z": 2.4960852396311436e-15, + "particle_weight": 1883777.3071500752 }, "ele2": { - "particle_opticalDepthQSR": 340.82229684726735, - "particle_position_x": 5.233059654483856e-06, - "particle_position_y": 4.781569085220371e-06, - "particle_position_z": 1.6293559425324337e-05, - "particle_momentum_x": 4.802611971470525e-18, - "particle_momentum_y": 4.3556825243407754e-18, - "particle_momentum_z": 2.41587659230925e-15, - "particle_weight": 2481137.860727036 + "particle_opticalDepthQSR": 391.0032239817431, + "particle_position_x": 5.442400321293304e-06, + "particle_position_y": 5.3621200149906786e-06, + "particle_position_z": 1.716649458876156e-05, + "particle_momentum_x": 4.7056270741663116e-18, + "particle_momentum_y": 4.414438514376292e-18, + "particle_momentum_z": 2.3816287305827113e-15, + "particle_weight": 2255238.1204111334 }, "pho1": { - "particle_opticalDepthBW": 9894.64959724129, - "particle_position_x": 0.00014191369823397644, - "particle_position_y": 0.00014347545392717968, - "particle_position_z": 0.00047442826029322116, + "particle_opticalDepthBW": 9912.332305353064, + "particle_position_x": 0.0001424857999095925, + "particle_position_y": 0.00014366766674682975, + "particle_position_z": 0.0004764478273708734, "particle_momentum_x": 0.0, "particle_momentum_y": 0.0, "particle_momentum_z": 0.0, - "particle_weight": 61063948.68610941 + "particle_weight": 61671336.44071695 }, "pho2": { - "particle_opticalDepthBW": 10292.1955840901, - "particle_position_x": 0.00014731892710321073, - "particle_position_y": 0.00014515617182809124, - "particle_position_z": 0.00048756513452074315, + "particle_opticalDepthBW": 10361.142481614075, + "particle_position_x": 0.00014718841339786984, + "particle_position_y": 0.00014538267635727008, + "particle_position_z": 0.00048768591004702896, "particle_momentum_x": 0.0, "particle_momentum_y": 0.0, "particle_momentum_z": 0.0, - "particle_weight": 62299636.622087024 + "particle_weight": 61889405.34553001 }, "pos1": { - "particle_opticalDepthQSR": 387.7441392212553, - "particle_position_x": 5.1462880118803425e-06, - "particle_position_y": 5.2613832293016684e-06, - "particle_position_z": 1.7054223425917483e-05, - "particle_momentum_x": 4.6437665862693495e-18, - "particle_momentum_y": 4.761862836969051e-18, - "particle_momentum_z": 2.3776996599289627e-15, - "particle_weight": 2625121.7841375084 + "particle_opticalDepthQSR": 343.0370605565725, + "particle_position_x": 5.41949101836001e-06, + "particle_position_y": 5.5292865311090835e-06, + "particle_position_z": 1.7063207973991194e-05, + "particle_momentum_x": 4.76704175252458e-18, + "particle_momentum_y": 5.049007427826035e-18, + "particle_momentum_z": 2.533578387785534e-15, + "particle_weight": 1821809.4309160584 }, "pos2": { - "particle_opticalDepthQSR": 361.943365907597, - "particle_position_x": 4.969019565031149e-06, - "particle_position_y": 4.361394970806125e-06, - "particle_position_z": 1.7413304358612675e-05, - "particle_momentum_x": 5.6348322786528905e-18, - "particle_momentum_y": 4.8171439953214205e-18, - "particle_momentum_z": 2.1937254860708963e-15, - "particle_weight": 2529794.7740602638 + "particle_opticalDepthQSR": 390.0211924195479, + "particle_position_x": 5.053169658257247e-06, + "particle_position_y": 5.039302796539821e-06, + "particle_position_z": 1.8526669235892217e-05, + "particle_momentum_x": 5.755216173987019e-18, + "particle_momentum_y": 5.056618594918483e-18, + "particle_momentum_z": 2.52324147594341e-15, + "particle_weight": 2328558.5523011778 } -} +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/ElectrostaticSphere.json b/Regression/Checksum/benchmarks_json/ElectrostaticSphere.json index 2d26136cc5e..4e9bde52a78 100644 --- a/Regression/Checksum/benchmarks_json/ElectrostaticSphere.json +++ b/Regression/Checksum/benchmarks_json/ElectrostaticSphere.json @@ -1,17 +1,17 @@ { "lev=0": { - "Ex": 6.504757806138167, - "Ey": 6.504757806138166, - "Ez": 6.501489620648247, - "rho": 2.60925680083338e-10 + "Ex": 6.504803536004636, + "Ey": 6.504803536004637, + "Ez": 6.504803536004637, + "rho": 2.609256800833379e-10 }, "electron": { - "particle_momentum_x": 1.0083594348819436e-23, - "particle_momentum_y": 1.0083594348819436e-23, - "particle_momentum_z": 1.0159789930728084e-23, - "particle_position_x": 518.4112403470117, - "particle_position_y": 518.4112403470117, - "particle_position_z": 520.2110302538686, + "particle_momentum_x": 1.0084159184928337e-23, + "particle_momentum_y": 1.008415918492834e-23, + "particle_momentum_z": 1.0084159184928338e-23, + "particle_position_x": 518.418176976446, + "particle_position_y": 518.4181769764461, + "particle_position_z": 518.418176976446, "particle_weight": 6212.501525878906 } } diff --git a/Regression/Checksum/benchmarks_json/ElectrostaticSphereLabFrame.json b/Regression/Checksum/benchmarks_json/ElectrostaticSphereLabFrame.json index d46a8cfcda4..648e805c938 100644 --- a/Regression/Checksum/benchmarks_json/ElectrostaticSphereLabFrame.json +++ b/Regression/Checksum/benchmarks_json/ElectrostaticSphereLabFrame.json @@ -1,17 +1,17 @@ { "lev=0": { - "Ex": 6.504586436288052, - "Ey": 6.504586436288051, - "Ez": 6.501896077968277, - "rho": 2.6092568008333797e-10 + "Ex": 6.504631859055906, + "Ey": 6.504631859055904, + "Ez": 6.504631859055904, + "rho": 2.609256800833379e-10 }, "electron": { - "particle_momentum_x": 9.881704085619681e-24, - "particle_momentum_y": 9.881704085619682e-24, - "particle_momentum_z": 9.969862126470096e-24, - "particle_position_x": 513.5147715983729, - "particle_position_y": 513.5147715983729, - "particle_position_z": 515.4187595183639, + "particle_momentum_x": 9.882423407883042e-24, + "particle_momentum_y": 9.882423407883041e-24, + "particle_momentum_z": 9.882423407883041e-24, + "particle_position_x": 513.522905933945, + "particle_position_y": 513.522905933945, + "particle_position_z": 513.522905933945, "particle_weight": 6212.501525878906 } } diff --git a/Regression/Checksum/benchmarks_json/ElectrostaticSphereRZ.json b/Regression/Checksum/benchmarks_json/ElectrostaticSphereRZ.json index 038ea5a6449..a94e8a4cd0b 100644 --- a/Regression/Checksum/benchmarks_json/ElectrostaticSphereRZ.json +++ b/Regression/Checksum/benchmarks_json/ElectrostaticSphereRZ.json @@ -1,17 +1,17 @@ { "lev=0": { - "Er": 0.11588932103656224, + "Er": 0.11588944816311227, "Et": 0.0, - "Ez": 0.11483033274834441, - "rho": 9.716216490680392e-12 + "Ez": 0.11492556391155323, + "rho": 9.71614990653583e-12 }, "electron": { - "particle_momentum_x": 4.54997005723341e-25, - "particle_momentum_y": 4.434893182689167e-25, - "particle_momentum_z": 6.989989085028144e-25, - "particle_position_x": 35.73642348640023, - "particle_position_y": 35.69851886414351, + "particle_momentum_x": 4.550172089304043e-25, + "particle_momentum_y": 4.435063230214811e-25, + "particle_momentum_z": 6.953391709067729e-25, + "particle_position_x": 35.73677703655427, + "particle_position_y": 35.61000307679542, "particle_theta": 823.0413054971611, "particle_weight": 6406.017620347038 } -} +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/EnergyConservingThermalPlasma.json b/Regression/Checksum/benchmarks_json/EnergyConservingThermalPlasma.json new file mode 100644 index 00000000000..3451ac42a9f --- /dev/null +++ b/Regression/Checksum/benchmarks_json/EnergyConservingThermalPlasma.json @@ -0,0 +1,29 @@ +{ + "lev=0": { + "Bx": 0.0, + "By": 0.0, + "Bz": 0.0, + "Ex": 14828575131693.746, + "Ey": 0.0, + "Ez": 15137612614715.773, + "jx": 5.572707185530917e+18, + "jy": 7.091118076558307e+18, + "jz": 6.668532755465409e+18 + }, + "protons": { + "particle_momentum_x": 3.187482462935171e-20, + "particle_momentum_y": 3.397549094666928e-20, + "particle_momentum_z": 3.300591612645518e-20, + "particle_position_x": 6.8052980975094785e-06, + "particle_position_y": 6.8011432612301085e-06, + "particle_weight": 2823958719279159.5 + }, + "electrons": { + "particle_momentum_x": 7.610483420803999e-22, + "particle_momentum_y": 8.024587481438078e-22, + "particle_momentum_z": 8.138108658606551e-22, + "particle_position_x": 6.8822971884226895e-06, + "particle_position_y": 6.8608357509110685e-06, + "particle_weight": 2823958719279159.5 + } +} \ No newline at end of file diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 4d1435e61de..9da04c2c6ca 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -4648,6 +4648,22 @@ doVis = 0 compareParticles = 1 analysisRoutine = Examples/Tests/Implicit/analysis_1d.py +[EnergyConservingThermalPlasma] +buildDir = . +inputFile = Examples/Tests/energy_conserving_thermal_plasma/inputs_2d_electrostatic +dim = 2 +addToCompileString = +cmakeSetupOpts = -DWarpX_DIMS=2 +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 0 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +analysisRoutine = Examples/Tests/energy_conserving_thermal_plasma/analysis.py + [focusing_gaussian_beam] buildDir = . inputFile = Examples/Tests/gaussian_beam/inputs_focusing_beam diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index b062501ed1b..dd6acfd2694 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -122,12 +122,12 @@ WarpX::Evolve (int numsteps) // At the beginning, we have B^{n} and E^{n}. // Particles have p^{n} and x^{n}. // is_synchronized is true. + if (is_synchronized) { - if (electrostatic_solver_id == ElectrostaticSolverAlgo::None) { - // Not called at each iteration, so exchange all guard cells - FillBoundaryE(guard_cells.ng_alloc_EB); - FillBoundaryB(guard_cells.ng_alloc_EB); - } + // Not called at each iteration, so exchange all guard cells + FillBoundaryE(guard_cells.ng_alloc_EB); + FillBoundaryB(guard_cells.ng_alloc_EB); + UpdateAuxilaryData(); FillBoundaryAux(guard_cells.ng_UpdateAux); // on first step, push p by -0.5*dt @@ -140,15 +140,15 @@ WarpX::Evolve (int numsteps) is_synchronized = false; } else { + // Beyond one step, we have E^{n} and B^{n}. + // Particles have p^{n-1/2} and x^{n}. + // E and B: enough guard cells to update Aux or call Field Gather in fp and cp + // Need to update Aux on lower levels, to interpolate to higher levels. + + // E and B are up-to-date inside the domain only + FillBoundaryE(guard_cells.ng_FieldGather); + FillBoundaryB(guard_cells.ng_FieldGather); if (electrostatic_solver_id == ElectrostaticSolverAlgo::None) { - // Beyond one step, we have E^{n} and B^{n}. - // Particles have p^{n-1/2} and x^{n}. - - // E and B are up-to-date inside the domain only - FillBoundaryE(guard_cells.ng_FieldGather); - FillBoundaryB(guard_cells.ng_FieldGather); - // E and B: enough guard cells to update Aux or call Field Gather in fp and cp - // Need to update Aux on lower levels, to interpolate to higher levels. if (fft_do_time_averaging) { FillBoundaryE_avg(guard_cells.ng_FieldGather);