diff --git a/Examples/Tests/PythonWrappers/PICMI_inputs_2d.py b/Examples/Tests/PythonWrappers/PICMI_inputs_2d.py index c34c9cf8035..2f9b169529d 100755 --- a/Examples/Tests/PythonWrappers/PICMI_inputs_2d.py +++ b/Examples/Tests/PythonWrappers/PICMI_inputs_2d.py @@ -256,7 +256,7 @@ def check_values(benchmark, data, rtol, atol): # E check_values(1013263608.6369569, Ex[:,:], rtol, atol) -check_values(717278253.4505507 , Ey[:,:], rtol, atol) +check_values(717278256.7957529 , Ey[:,:], rtol, atol) check_values(717866566.5718911 , Ez[:,:], rtol, atol) # B check_values(3.0214509313437636, Bx[:,:], rtol, atol) @@ -267,22 +267,22 @@ def check_values(benchmark, data, rtol, atol): check_values(1013672631.8764204, G[:,:], rtol, atol) # E in PML check_values(364287936.1526477 , Expml[:,:,0], rtol, atol) -check_values(183582351.3212558 , Expml[:,:,1], rtol, atol) +check_values(183582352.20753333, Expml[:,:,1], rtol, atol) check_values(190065766.41491824, Expml[:,:,2], rtol, atol) -check_values(440581905.9336025 , Eypml[:,:,0], rtol, atol) -check_values(178117293.6629357 , Eypml[:,:,1], rtol, atol) +check_values(440581907.0828975 , Eypml[:,:,0], rtol, atol) +check_values(178117294.05871135, Eypml[:,:,1], rtol, atol) check_values(0.0 , Eypml[:,:,2], rtol, atol) check_values(430277101.26568377, Ezpml[:,:,0], rtol, atol) check_values(0.0 , Ezpml[:,:,1], rtol, atol) check_values(190919663.2167449 , Ezpml[:,:,2], rtol, atol) # B in PML check_values(1.0565189315366146 , Bxpml[:,:,0], rtol, atol) -check_values(0.4618191395098556 , Bxpml[:,:,1], rtol, atol) -check_values(0.6849858273929585 , Bxpml[:,:,2], rtol, atol) +check_values(0.46181913800643065, Bxpml[:,:,1], rtol, atol) +check_values(0.6849858305343736 , Bxpml[:,:,2], rtol, atol) check_values(1.7228584190213505 , Bypml[:,:,0], rtol, atol) -check_values(0.47697331996765685, Bypml[:,:,1], rtol, atol) +check_values(0.47697332248020935, Bypml[:,:,1], rtol, atol) check_values(0.0 , Bypml[:,:,2], rtol, atol) -check_values(1.5183380774611628 , Bzpml[:,:,0], rtol, atol) +check_values(1.518338068658267 , Bzpml[:,:,0], rtol, atol) check_values(0.0 , Bzpml[:,:,1], rtol, atol) check_values(0.6849858291863835 , Bzpml[:,:,2], rtol, atol) # F and G in PML diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR.json index 075fe19afb8..9b4ab5a5401 100644 --- a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR.json +++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR.json @@ -2,43 +2,43 @@ "electrons": { "particle_cpu": 32768.0, "particle_id": 1123057664.0, - "particle_momentum_x": 4.241495334522421e-20, + "particle_momentum_x": 4.245553180462133e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 4.241495334522413e-20, - "particle_position_x": 0.6553605577217855, - "particle_position_y": 0.6553605577217853, + "particle_momentum_z": 4.245553180462118e-20, + "particle_position_x": 0.6553607116882387, + "particle_position_y": 0.6553607116882386, "particle_weight": 3200000000000000.5 }, "lev=0": { "Bx": 0.0, - "By": 32.88376813880145, + "By": 35.447497788711175, "Bz": 0.0, - "Ex": 7572926260078.977, + "Ex": 7569456532381.726, "Ey": 0.0, - "Ez": 7572926260078.995, - "jx": 7298077411371068.0, + "Ez": 7569456532381.752, + "jx": 7305055873426885.0, "jy": 0.0, - "jz": 7298077411371252.0 + "jz": 7305055873426942.0 }, "lev=1": { "Bx": 0.0, - "By": 565.8504447463731, + "By": 640.5385076158955, "Bz": 0.0, - "Ex": 7546287967911.542, + "Ex": 7591603030109.285, "Ey": 0.0, - "Ez": 7546287967911.541, - "jx": 7199765368042636.0, + "Ez": 7591603030109.295, + "jx": 7225962509586529.0, "jy": 0.0, - "jz": 7199765368042662.0 + "jz": 7225962509586578.0 }, "positrons": { "particle_cpu": 32768.0, "particle_id": 3371204608.0, - "particle_momentum_x": 4.2413987484133116e-20, + "particle_momentum_x": 4.2453138721853093e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 4.241398748413304e-20, - "particle_position_x": 0.6553600696643818, - "particle_position_y": 0.6553600696643818, + "particle_momentum_z": 4.245313872185295e-20, + "particle_position_x": 0.6553599057578553, + "particle_position_y": 0.6553599057578553, "particle_weight": 3200000000000000.5 } } \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR_anisotropic.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR_anisotropic.json index c79b4dbbf77..0ec7ff20c85 100644 --- a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR_anisotropic.json +++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR_anisotropic.json @@ -2,43 +2,43 @@ "electrons": { "particle_cpu": 32768.0, "particle_id": 1123057664.0, - "particle_momentum_x": 4.2409233886523047e-20, + "particle_momentum_x": 4.2441494210115385e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 4.239636641708783e-20, - "particle_position_x": 0.6553604498033957, - "particle_position_y": 0.6553602617123965, + "particle_momentum_z": 4.2429971697106366e-20, + "particle_position_x": 0.6553605498557016, + "particle_position_y": 0.655360420612517, "particle_weight": 3200000000000000.5 }, "lev=0": { "Bx": 0.0, - "By": 29.033843714475076, + "By": 32.13899895667297, "Bz": 0.0, - "Ex": 7575734617759.523, + "Ex": 7573285869190.605, "Ey": 0.0, - "Ez": 7575399948801.469, - "jx": 7296519320465208.0, + "Ez": 7572536902181.709, + "jx": 7302040020242814.0, "jy": 0.0, - "jz": 7297090426947096.0 + "jz": 7302892615117250.0 }, "lev=1": { "Bx": 0.0, - "By": 71.46369739075365, + "By": 74.04987113513879, "Bz": 0.0, - "Ex": 4602033610493.496, + "Ex": 4601660582425.954, "Ey": 0.0, - "Ez": 7017735833493.598, - "jx": 4492590664379721.0, + "Ez": 7012131128423.561, + "jx": 4495212217937750.0, "jy": 0.0, - "jz": 6825856952745953.0 + "jz": 6848672875762947.0 }, "positrons": { "particle_cpu": 32768.0, "particle_id": 3371204608.0, - "particle_momentum_x": 4.240515207391037e-20, + "particle_momentum_x": 4.2436090244998326e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 4.2396984750798293e-20, - "particle_position_x": 0.6553601899702647, - "particle_position_y": 0.6553597467035968, + "particle_momentum_z": 4.2430343380944015e-20, + "particle_position_x": 0.6553600813503462, + "particle_position_y": 0.6553595861949815, "particle_weight": 3200000000000000.5 } } \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR_psatd.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR_psatd.json index 435afcd7372..ca13da0324b 100644 --- a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR_psatd.json +++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_MR_psatd.json @@ -2,43 +2,43 @@ "electrons": { "particle_cpu": 32768.0, "particle_id": 1123057664.0, - "particle_momentum_x": 4.2362974799075715e-20, + "particle_momentum_x": 4.237978127636407e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 4.2362974818601834e-20, - "particle_position_x": 0.6553603056279188, - "particle_position_y": 0.6553603056279399, + "particle_momentum_z": 4.2379781109485186e-20, + "particle_position_x": 0.655360346626376, + "particle_position_y": 0.655360346625584, "particle_weight": 3200000000000000.5 }, "lev=0": { "Bx": 0.0, - "By": 49.400846503194956, + "By": 49.13289051610905, "Bz": 0.0, - "Ex": 7590421951518.164, + "Ex": 7589572073848.546, "Ey": 0.0, - "Ez": 7590421950277.463, - "jx": 7279591576517866.0, + "Ez": 7589572085733.721, + "jx": 7282449902618536.0, "jy": 0.0, - "jz": 7279591579334755.0 + "jz": 7282449869844108.0 }, "lev=1": { "Bx": 0.0, - "By": 564.3980303336011, + "By": 623.630930990979, "Bz": 0.0, - "Ex": 7563548313704.213, + "Ex": 7590754010944.117, "Ey": 0.0, - "Ez": 7563548287845.812, - "jx": 6584893500547710.0, + "Ez": 7590753791450.942, + "jx": 6596452748155826.0, "jy": 0.0, - "jz": 6584893531174318.0 + "jz": 6596452562357128.0 }, "positrons": { "particle_cpu": 32768.0, "particle_id": 3371204608.0, - "particle_momentum_x": 4.236433706251501e-20, + "particle_momentum_x": 4.2381117997129824e-20, "particle_momentum_y": 0.0, - "particle_momentum_z": 4.23643370757891e-20, - "particle_position_x": 0.6553599973488732, - "particle_position_y": 0.6553599973488947, + "particle_momentum_z": 4.23811177887596e-20, + "particle_position_x": 0.6553599570597749, + "particle_position_y": 0.6553599570604656, "particle_weight": 3200000000000000.5 } } \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/LaserAccelerationMR.json b/Regression/Checksum/benchmarks_json/LaserAccelerationMR.json index cfee789f139..55c39432f07 100644 --- a/Regression/Checksum/benchmarks_json/LaserAccelerationMR.json +++ b/Regression/Checksum/benchmarks_json/LaserAccelerationMR.json @@ -33,14 +33,14 @@ }, "lev=1": { "Bx": 53303280.78531405, - "By": 213761.8564937196, + "By": 213820.3084522391, "Bz": 1321618.978929098, - "Ex": 47856565901964.84, + "Ex": 47846019388623.05, "Ey": 1.455101409934242e+16, - "Ez": 65304559419893.03, + "Ez": 65304502775874.59, "jx": 46169274324016.0, "jy": 9.274837523569297e+18, "jz": 1.926277235505425e+16, "rho": 67304284.98422323 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/LaserOnFine.json b/Regression/Checksum/benchmarks_json/LaserOnFine.json index 992ada0ce03..fbaac2ae66a 100644 --- a/Regression/Checksum/benchmarks_json/LaserOnFine.json +++ b/Regression/Checksum/benchmarks_json/LaserOnFine.json @@ -14,15 +14,15 @@ }, "lev=1": { "Bx": 0.0, - "By": 1234346.3811825265, + "By": 1234425.838348821, "Bz": 0.0, - "Ex": 369321181489775.1, + "Ex": 369306119506292.0, "Ey": 0.0, - "Ez": 17591468778544.74, + "Ez": 17592228452124.82, "divB": 0.0, "jx": 2.4549619817062897e+18, "jy": 0.0, "jz": 0.0, "part_per_cell": 102.0 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/PEC_field_mr.json b/Regression/Checksum/benchmarks_json/PEC_field_mr.json index c6d84b01123..99966cbd257 100644 --- a/Regression/Checksum/benchmarks_json/PEC_field_mr.json +++ b/Regression/Checksum/benchmarks_json/PEC_field_mr.json @@ -4,7 +4,7 @@ "Ey": 1052208601.7204809 }, "lev=1": { - "Bx": 0.36174107941906025, - "Ey": 94466060.63059175 + "Bx": 0.3714838555919295, + "Ey": 98224588.87142882 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/PlasmaAccelerationMR.json b/Regression/Checksum/benchmarks_json/PlasmaAccelerationMR.json index 13f65226f03..6b0b68729ad 100644 --- a/Regression/Checksum/benchmarks_json/PlasmaAccelerationMR.json +++ b/Regression/Checksum/benchmarks_json/PlasmaAccelerationMR.json @@ -31,12 +31,12 @@ "jz": 4235811839330752.0 }, "lev=1": { - "Bx": 4.147877077278406, - "By": 142707.7156173306, - "Bz": 2.536603360961009, - "Ex": 49284312088805.66, - "Ey": 1421188611.560146, - "Ez": 52701506458469.23, + "Bx": 4.150715838073591, + "By": 142753.1342940014, + "Bz": 2.536506757005651, + "Ex": 49292890789333.12, + "Ey": 1422005008.958596, + "Ez": 52701928262289.76, "jx": 3029834005553842.0, "jy": 92336633890.84117, "jz": 4286789664067064.0 @@ -51,4 +51,4 @@ "particle_position_y": 0.1035059304082606, "particle_weight": 823974609374999.9 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/Python_LaserAccelerationMR.json b/Regression/Checksum/benchmarks_json/Python_LaserAccelerationMR.json index cfee789f139..55c39432f07 100644 --- a/Regression/Checksum/benchmarks_json/Python_LaserAccelerationMR.json +++ b/Regression/Checksum/benchmarks_json/Python_LaserAccelerationMR.json @@ -33,14 +33,14 @@ }, "lev=1": { "Bx": 53303280.78531405, - "By": 213761.8564937196, + "By": 213820.3084522391, "Bz": 1321618.978929098, - "Ex": 47856565901964.84, + "Ex": 47846019388623.05, "Ey": 1.455101409934242e+16, - "Ez": 65304559419893.03, + "Ez": 65304502775874.59, "jx": 46169274324016.0, "jy": 9.274837523569297e+18, "jz": 1.926277235505425e+16, "rho": 67304284.98422323 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/Python_wrappers.json b/Regression/Checksum/benchmarks_json/Python_wrappers.json index 9cbe4988911..c67649c485c 100644 --- a/Regression/Checksum/benchmarks_json/Python_wrappers.json +++ b/Regression/Checksum/benchmarks_json/Python_wrappers.json @@ -2,9 +2,9 @@ "lev=0": { "Bx": 1.3641435093370002, "By": 2.064187020248739, - "Bz": 1.3641435106886308, - "Ex": 424983256.20110655, - "Ey": 159848272.83698195, - "Ez": 321466524.9277822 + "Bz": 1.364143507950999, + "Ex": 424983254.6320197, + "Ey": 159848274.5435922, + "Ez": 321466522.8953871 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/RefinedInjection.json b/Regression/Checksum/benchmarks_json/RefinedInjection.json index a13b6eb224b..975d2b8c5dc 100644 --- a/Regression/Checksum/benchmarks_json/RefinedInjection.json +++ b/Regression/Checksum/benchmarks_json/RefinedInjection.json @@ -33,14 +33,14 @@ }, "lev=1": { "Bx": 53303352.85256407, - "By": 229340.2527555132, + "By": 229398.7047140327, "Bz": 1321676.328125808, - "Ex": 52873151117961.77, + "Ex": 52862604604619.98, "Ey": 1.455090881912273e+16, - "Ez": 67855870768719.61, + "Ez": 67855814124701.18, "jx": 121460903131204.2, "jy": 9.274942706868505e+18, "jz": 1.927066428524446e+16, "rho": 76242140.83674599 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/momentum-conserving-gather.json b/Regression/Checksum/benchmarks_json/momentum-conserving-gather.json index e3862ad0949..7837d1d7436 100644 --- a/Regression/Checksum/benchmarks_json/momentum-conserving-gather.json +++ b/Regression/Checksum/benchmarks_json/momentum-conserving-gather.json @@ -31,12 +31,12 @@ "jz": 4304296344908232.0 }, "lev=1": { - "Bx": 4.06321359891687, - "By": 137867.46535533897, - "Bz": 2.398740968267233, - "Ex": 44258985973031.625, - "Ey": 1433041921.7273414, - "Ez": 51834767989855.24, + "Bx": 4.066233386167275, + "By": 137914.0920350435, + "Bz": 2.398664210322969, + "Ex": 44268789869454.29, + "Ey": 1433858315.947112, + "Ez": 51835208819539.50, "jx": 2880554704671068.0, "jy": 93913414121.91959, "jz": 4418516923648772.0 @@ -51,4 +51,4 @@ "particle_position_y": 0.1035024179589788, "particle_weight": 823974609375002.8 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json b/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json index f497b1b0cb3..6cc7803cb56 100644 --- a/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json +++ b/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json @@ -1,11 +1,11 @@ { "lev=0": { - "Bx": 2.9025658853498176e-06, - "By": 2.9007903785438975e-06, - "Bz": 2.9006439095978047e-06, - "Ex": 1639.2657294711155, - "Ey": 1648.9781131837074, - "Ez": 1649.012795508662, + "Bx": 2.916821611439007e-06, + "By": 2.905967588894673e-06, + "Bz": 2.910675767589761e-06, + "Ex": 1638.225417535941, + "Ey": 1643.469685733573, + "Ez": 1645.623696125757, "rho": 0.00017408482555034496 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/pml_x_psatd.json b/Regression/Checksum/benchmarks_json/pml_x_psatd.json index 15356580f5e..a9cf4abde80 100644 --- a/Regression/Checksum/benchmarks_json/pml_x_psatd.json +++ b/Regression/Checksum/benchmarks_json/pml_x_psatd.json @@ -1,12 +1,12 @@ { "lev=0": { - "Bx": 1.1599384673077106e-08, - "By": 1.515799977053999e-08, - "Bz": 8.995717690678642e-09, - "Ex": 3.9363123107318456, - "Ey": 4.134487508968383, - "Ez": 3.3070809317529846, - "divE": 189105.63427643626, + "Bx": 1.1599355635528331e-08, + "By": 1.5158011001988313e-08, + "Bz": 8.995705405879026e-09, + "Ex": 3.9363155460110746, + "Ey": 4.134482727292589, + "Ez": 3.3070827418135487, + "divE": 189105.46307809104, "rho": 1.6734952000013346e-06 } } \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/subcyclingMR.json b/Regression/Checksum/benchmarks_json/subcyclingMR.json index 1e1b3022632..6950d755e4a 100644 --- a/Regression/Checksum/benchmarks_json/subcyclingMR.json +++ b/Regression/Checksum/benchmarks_json/subcyclingMR.json @@ -2,7 +2,7 @@ "beam": { "particle_cpu": 0.0, "particle_id": 150005000.0, - "particle_momentum_x": 4.341970536840506e-19, + "particle_momentum_x": 4.341970536840507e-19, "particle_momentum_y": 0.0, "particle_momentum_z": 4.854125969816336e-17, "particle_position_x": 0.000631595029496105, @@ -21,44 +21,44 @@ }, "lev=0": { "Bx": 0.0, - "By": 999420.5203867452, + "By": 999414.3581034234, "Bz": 0.0, - "Ex": 388788109636517.25, + "Ex": 388778725968159.1, "Ey": 0.0, - "Ez": 475734042024292.44, - "jx": 3.8453833238655616e+17, + "Ez": 475738753763627.9, + "jx": 3.8453235309975706e+17, "jy": 0.0, - "jz": 6.358844776737537e+17 + "jz": 6.358921323286822e+17 }, "lev=1": { "Bx": 0.0, - "By": 1869002.1405830402, + "By": 1869003.9426403372, "Bz": 0.0, - "Ex": 876808143853005.0, + "Ex": 876798193970146.4, "Ey": 0.0, - "Ez": 918405586302947.2, - "jx": 1.536628245883072e+17, + "Ez": 918415238562273.8, + "jx": 1.536648937833682e+17, "jy": 0.0, - "jz": 4.129776951264712e+17 + "jz": 4.129783535034742e+17 }, "plasma_e": { "particle_cpu": 30217.0, "particle_id": 1048674286.0, - "particle_momentum_x": 1.10240793316972e-18, + "particle_momentum_x": 1.1023996684852916e-18, "particle_momentum_y": 0.0, - "particle_momentum_z": 1.430691945336632e-18, - "particle_position_x": 0.23607785290428887, - "particle_position_y": 0.290297172085528, + "particle_momentum_z": 1.430708750836071e-18, + "particle_position_x": 0.23607832427125558, + "particle_position_y": 0.29029734645671285, "particle_weight": 5532897949218750.0 }, "plasma_p": { "particle_cpu": 31872.0, "particle_id": 1058217024.0, - "particle_momentum_x": 1.8564730642732005e-18, + "particle_momentum_x": 1.856467736176977e-18, "particle_momentum_y": 0.0, - "particle_momentum_z": 2.416764596133168e-18, - "particle_position_x": 0.23904668587287342, - "particle_position_y": 0.31004658592211964, + "particle_momentum_z": 2.416775960843685e-18, + "particle_position_x": 0.2390466856387225, + "particle_position_y": 0.31004658573281996, "particle_weight": 5835937500000001.0 } } \ No newline at end of file diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 2783a6bf651..fafbd9259bf 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -42,7 +42,12 @@ struct Sigma : amrex::Gpu::DeviceVector struct SigmaBox { SigmaBox (const amrex::Box& box, const amrex::BoxArray& grids, - const amrex::Real* dx, int ncell, int delta); + const amrex::Real* dx, int ncell, int delta, const amrex::Box& regdomain); + + void define_single (const amrex::Box& regdomain, int ncell, + const amrex::Array& fac); + void define_multiple (const amrex::Box& box, const amrex::BoxArray& grids, int ncell, + const amrex::Array& fac); void ComputePMLFactorsB (const amrex::Real* dx, amrex::Real dt); void ComputePMLFactorsE (const amrex::Real* dx, amrex::Real dt); @@ -64,8 +69,9 @@ class SigmaBoxFactory : public amrex::FabFactory { public: - SigmaBoxFactory (const amrex::BoxArray& grid_ba, const amrex::Real* dx, int ncell, int delta) - : m_grids(grid_ba), m_dx(dx), m_ncell(ncell), m_delta(delta) {} + SigmaBoxFactory (const amrex::BoxArray& grid_ba, const amrex::Real* dx, int ncell, int delta, + const amrex::Box& regular_domain) + : m_grids(grid_ba), m_dx(dx), m_ncell(ncell), m_delta(delta), m_regdomain(regular_domain) {} virtual ~SigmaBoxFactory () = default; SigmaBoxFactory (const SigmaBoxFactory&) = default; @@ -77,7 +83,7 @@ public: virtual SigmaBox* create (const amrex::Box& box, int /*ncomps*/, const amrex::FabInfo& /*info*/, int /*box_index*/) const final - { return new SigmaBox(box, m_grids, m_dx, m_ncell, m_delta); } + { return new SigmaBox(box, m_grids, m_dx, m_ncell, m_delta, m_regdomain); } virtual void destroy (SigmaBox* fab) const final { delete fab; } @@ -89,6 +95,7 @@ private: const amrex::Real* m_dx; int m_ncell; int m_delta; + amrex::Box m_regdomain; }; class MultiSigmaBox @@ -96,7 +103,8 @@ class MultiSigmaBox { public: MultiSigmaBox(const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, - const amrex::BoxArray& grid_ba, const amrex::Real* dx, int ncell, int delta); + const amrex::BoxArray& grid_ba, const amrex::Real* dx, int ncell, int delta, + const amrex::Box& regular_domain); void ComputePMLFactorsB (const amrex::Real* dx, amrex::Real dt); void ComputePMLFactorsE (const amrex::Real* dx, amrex::Real dt); private: @@ -238,11 +246,24 @@ private: } #endif - static amrex::BoxArray MakeBoxArray (const amrex::Geometry& geom, + static amrex::BoxArray MakeBoxArray (bool single_box_domain, + const amrex::Box& regular_domain, + const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell, int do_pml_in_domain, - const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), - const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); + const amrex::IntVect do_pml_Lo, + const amrex::IntVect do_pml_Hi); + + static amrex::BoxArray MakeBoxArray_single (const amrex::Box& regular_domain, + const amrex::BoxArray& grid_ba, + int ncell, const amrex::IntVect do_pml_Lo, + const amrex::IntVect do_pml_Hi); + + static amrex::BoxArray MakeBoxArray_multiple (const amrex::Geometry& geom, + const amrex::BoxArray& grid_ba, + int ncell, int do_pml_in_domain, + const amrex::IntVect do_pml_Lo, + const amrex::IntVect do_pml_Hi); static void CopyToPML (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); }; diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 837f66c2e6f..34b5b68c499 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -44,6 +44,7 @@ #include #include +#include #include #include #ifdef AMREX_USE_EB @@ -54,15 +55,12 @@ using namespace amrex; namespace { - static void FillLo (int idim, Sigma& sigma, Sigma& sigma_cumsum, + static void FillLo (Sigma& sigma, Sigma& sigma_cumsum, Sigma& sigma_star, Sigma& sigma_star_cumsum, - const Box& overlap, const Box& grid, Real fac) + const int olo, const int ohi, const int glo, Real fac) { - int glo = grid.smallEnd(idim); - int olo = overlap.smallEnd(idim); - int ohi = overlap.bigEnd(idim); - int slo = sigma.m_lo; - int sslo = sigma_star.m_lo; + const int slo = sigma.m_lo; + const int sslo = sigma_star.m_lo; const int N = ohi+1-olo+1; Real* p_sigma = sigma.data(); @@ -85,15 +83,12 @@ namespace }); } - static void FillHi (int idim, Sigma& sigma, Sigma& sigma_cumsum, + static void FillHi (Sigma& sigma, Sigma& sigma_cumsum, Sigma& sigma_star, Sigma& sigma_star_cumsum, - const Box& overlap, const Box& grid, Real fac) + const int olo, const int ohi, const int ghi, Real fac) { - int ghi = grid.bigEnd(idim); - int olo = overlap.smallEnd(idim); - int ohi = overlap.bigEnd(idim); - int slo = sigma.m_lo; - int sslo = sigma_star.m_lo; + const int slo = sigma.m_lo; + const int sslo = sigma_star.m_lo; const int N = ohi+1-olo+1; Real* p_sigma = sigma.data(); @@ -115,14 +110,12 @@ namespace } #if (AMREX_SPACEDIM != 1) - static void FillZero (int idim, Sigma& sigma, Sigma& sigma_cumsum, + static void FillZero (Sigma& sigma, Sigma& sigma_cumsum, Sigma& sigma_star, Sigma& sigma_star_cumsum, - const Box& overlap) + const int olo, const int ohi) { - int olo = overlap.smallEnd(idim); - int ohi = overlap.bigEnd(idim); - int slo = sigma.m_lo; - int sslo = sigma_star.m_lo; + const int slo = sigma.m_lo; + const int sslo = sigma_star.m_lo; const int N = ohi+1-olo+1; Real* p_sigma = sigma.data(); @@ -144,7 +137,8 @@ namespace } -SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int ncell, int delta) +SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int ncell, int delta, + const amrex::Box& regdomain) { BL_ASSERT(box.cellCentered()); @@ -154,14 +148,14 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - sigma [idim].resize(sz[idim]+1); - sigma_cumsum [idim].resize(sz[idim]+1); - sigma_star [idim].resize(sz[idim]+1); - sigma_star_cumsum [idim].resize(sz[idim]+1); - sigma_fac [idim].resize(sz[idim]+1); - sigma_cumsum_fac [idim].resize(sz[idim]+1); - sigma_star_fac [idim].resize(sz[idim]+1); - sigma_star_cumsum_fac[idim].resize(sz[idim]+1); + sigma [idim].resize(sz[idim]+1,std::numeric_limits::quiet_NaN()); + sigma_cumsum [idim].resize(sz[idim]+1,std::numeric_limits::quiet_NaN()); + sigma_star [idim].resize(sz[idim]+1,std::numeric_limits::quiet_NaN()); + sigma_star_cumsum [idim].resize(sz[idim]+1,std::numeric_limits::quiet_NaN()); + sigma_fac [idim].resize(sz[idim]+1,std::numeric_limits::quiet_NaN()); + sigma_cumsum_fac [idim].resize(sz[idim]+1,std::numeric_limits::quiet_NaN()); + sigma_star_fac [idim].resize(sz[idim]+1,std::numeric_limits::quiet_NaN()); + sigma_star_cumsum_fac[idim].resize(sz[idim]+1,std::numeric_limits::quiet_NaN()); sigma [idim].m_lo = lo[idim]; sigma [idim].m_hi = hi[idim]+1; @@ -186,6 +180,58 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n fac[idim] = 4.0_rt*PhysConst::c/(dx[idim]*static_cast(delta*delta)); } + if (regdomain.ok()) { // The union of the regular grids is a single box + define_single(regdomain, ncell, fac); + } else { + define_multiple(box, grids, ncell, fac); + } +} + +void SigmaBox::define_single (const Box& regdomain, int ncell, + const Array& fac) +{ + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + const int slo = sigma[idim].lo(); + const int shi = sigma[idim].hi()-1; + const int dlo = regdomain.smallEnd(idim); + const int dhi = regdomain.bigEnd(idim); + + // Lo + int olo = std::max(slo, dlo-ncell); + int ohi = std::min(shi, dlo-1); + if (ohi >= olo) { + FillLo(sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + olo, ohi, dlo, fac[idim]); + } + +#if (AMREX_SPACEDIM != 1) + // Mid + olo = std::max(slo, dlo); + ohi = std::min(shi, dhi); + if (ohi >= olo) { + FillZero(sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + olo, ohi); + } +#endif + + // Hi + olo = std::max(slo, dhi+1); + ohi = std::min(shi, dhi+ncell); + if (ohi >= olo) { + FillHi(sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + olo, ohi, dhi, fac[idim]); + } + } + + amrex::Gpu::streamSynchronize(); +} + +void SigmaBox::define_multiple (const Box& box, const BoxArray& grids, int ncell, + const Array& fac) +{ const std::vector >& isects = grids.intersections(box, false, ncell); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) @@ -251,9 +297,10 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n #endif Box looverlap = lobox & box; if (looverlap.ok()) { - FillLo(idim, sigma[idim], sigma_cumsum[idim], + FillLo(sigma[idim], sigma_cumsum[idim], sigma_star[idim], sigma_star_cumsum[idim], - looverlap, grid_box, fac[idim]); + looverlap.smallEnd(idim), looverlap.bigEnd(idim), + grid_box.smallEnd(idim), fac[idim]); } Box hibox = amrex::adjCellHi(grid_box, idim, ncell); @@ -263,9 +310,10 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n #endif Box hioverlap = hibox & box; if (hioverlap.ok()) { - FillHi(idim, sigma[idim], sigma_cumsum[idim], + FillHi(sigma[idim], sigma_cumsum[idim], sigma_star[idim], sigma_star_cumsum[idim], - hioverlap, grid_box, fac[idim]); + hioverlap.smallEnd(idim), hioverlap.bigEnd(idim), + grid_box.bigEnd(idim), fac[idim]); } if (!looverlap.ok() && !hioverlap.ok()) { @@ -280,8 +328,9 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n const Box& grid_box = grids[gid]; const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box; if (overlap.ok()) { - FillZero(idim, sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], overlap); + FillZero(sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + overlap.smallEnd(idim), overlap.bigEnd(idim)); } else { amrex::Abort("SigmaBox::SigmaBox(): side_side_edges, how did this happen?\n"); @@ -295,17 +344,19 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n Box lobox = amrex::adjCellLo(grid_box, idim, ncell); Box looverlap = lobox.grow(jdim,ncell).grow(kdim,ncell) & box; if (looverlap.ok()) { - FillLo(idim, sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], - looverlap, grid_box, fac[idim]); + FillLo(sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + looverlap.smallEnd(idim), looverlap.bigEnd(idim), + grid_box.smallEnd(idim), fac[idim]); } Box hibox = amrex::adjCellHi(grid_box, idim, ncell); Box hioverlap = hibox.grow(jdim,ncell).grow(kdim,ncell) & box; if (hioverlap.ok()) { - FillHi(idim, sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], - hioverlap, grid_box, fac[idim]); + FillHi(sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + hioverlap.smallEnd(idim), hioverlap.bigEnd(idim), + grid_box.bigEnd(idim), fac[idim]); } if (!looverlap.ok() && !hioverlap.ok()) { @@ -324,8 +375,9 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box; #endif if (overlap.ok()) { - FillZero(idim, sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], overlap); + FillZero(sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + overlap.smallEnd(idim), overlap.bigEnd(idim)); } else { amrex::Abort("SigmaBox::SigmaBox(): side_faces, how did this happen?\n"); } @@ -339,17 +391,19 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n const Box& lobox = amrex::adjCellLo(grid_box, idim, ncell); Box looverlap = lobox & box; if (looverlap.ok()) { - FillLo(idim, sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], - looverlap, grid_box, fac[idim]); + FillLo(sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + looverlap.smallEnd(idim), looverlap.bigEnd(idim), + grid_box.smallEnd(idim), fac[idim]); } const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell); Box hioverlap = hibox & box; if (hioverlap.ok()) { - FillHi(idim, sigma[idim], sigma_cumsum[idim], - sigma_star[idim], sigma_star_cumsum[idim], - hioverlap, grid_box, fac[idim]); + FillHi(sigma[idim], sigma_cumsum[idim], + sigma_star[idim], sigma_star_cumsum[idim], + hioverlap.smallEnd(idim), hioverlap.bigEnd(idim), + grid_box.bigEnd(idim), fac[idim]); } if (!looverlap.ok() && !hioverlap.ok()) { @@ -362,7 +416,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n } } - amrex::Gpu::synchronize(); + amrex::Gpu::streamSynchronize(); } @@ -435,9 +489,10 @@ SigmaBox::ComputePMLFactorsE (const Real* a_dx, Real dt) } MultiSigmaBox::MultiSigmaBox (const BoxArray& ba, const DistributionMapping& dm, - const BoxArray& grid_ba, const Real* dx, int ncell, int delta) + const BoxArray& grid_ba, const Real* dx, int ncell, int delta, + const amrex::Box& regular_domain) : FabArray(ba,dm,1,0,MFInfo(), - SigmaBoxFactory(grid_ba,dx,ncell,delta)) + SigmaBoxFactory(grid_ba,dx,ncell,delta, regular_domain)) {} void @@ -495,19 +550,23 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri // Note that this is okay to build pml inside domain for a single patch, or joint patches // with same [min,max]. But it does not support multiple disjoint refinement patches. Box domain0 = grid_ba.minimalBox(); - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - if (do_pml_Lo[idim]){ - domain0.growLo(idim, -ncell); - } - if (do_pml_Hi[idim]){ - domain0.growHi(idim, -ncell); + if (do_pml_in_domain) { + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (do_pml_Lo[idim]){ + domain0.growLo(idim, -ncell); + } + if (do_pml_Hi[idim]){ + domain0.growHi(idim, -ncell); + } } } - const BoxArray grid_ba_reduced = BoxArray(grid_ba.boxList().intersect(domain0)); + const BoxArray grid_ba_reduced = (do_pml_in_domain) ? + BoxArray(grid_ba.boxList().intersect(domain0)) : grid_ba; + + bool is_single_box_domain = domain0.numPts() == grid_ba_reduced.numPts(); + const BoxArray& ba = MakeBoxArray(is_single_box_domain, domain0, *geom, grid_ba_reduced, ncell, + do_pml_in_domain, do_pml_Lo, do_pml_Hi); - const BoxArray& ba = (do_pml_in_domain)? - MakeBoxArray(*geom, grid_ba_reduced, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi) : - MakeBoxArray(*geom, grid_ba, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi); if (ba.empty()) { m_ok = false; return; @@ -660,12 +719,10 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri pml_G_fp->setVal(0.0); } - if (do_pml_in_domain){ - sigba_fp = std::make_unique(ba, dm, grid_ba_reduced, geom->CellSize(), ncell, delta); - } - else { - sigba_fp = std::make_unique(ba, dm, grid_ba, geom->CellSize(), ncell, delta); - } + Box single_domain_box = is_single_box_domain ? domain0 : Box(); + // Empty box (i.e., Box()) means it's not a single box domain. + sigba_fp = std::make_unique(ba, dm, grid_ba_reduced, geom->CellSize(), + ncell, delta, single_domain_box); if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) { #ifndef WARPX_USE_PSATD @@ -705,26 +762,31 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri grid_cba.coarsen(ref_ratio); // assuming that the bounding box around grid_cba is a single patch, and not disjoint patches, similar to fine patch. - amrex::Box domain1 = grid_cba.minimalBox(); - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - if (do_pml_Lo[idim]){ - // ncell is divided by refinement ratio to ensure that the - // physical width of the PML region is equal in fine and coarse patch - domain1.growLo(idim, -ncell/ref_ratio[idim]); - } - if (do_pml_Hi[idim]){ - // ncell is divided by refinement ratio to ensure that the - // physical width of the PML region is equal in fine and coarse patch - domain1.growHi(idim, -ncell/ref_ratio[idim]); + amrex::Box cdomain = grid_cba.minimalBox(); + if (do_pml_in_domain) { + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (do_pml_Lo[idim]){ + // ncell is divided by refinement ratio to ensure that the + // physical width of the PML region is equal in fine and coarse patch + cdomain.growLo(idim, -ncell/ref_ratio[idim]); + } + if (do_pml_Hi[idim]){ + // ncell is divided by refinement ratio to ensure that the + // physical width of the PML region is equal in fine and coarse patch + cdomain.growHi(idim, -ncell/ref_ratio[idim]); + } } } - const BoxArray grid_cba_reduced = BoxArray(grid_cba.boxList().intersect(domain1)); + const BoxArray grid_cba_reduced = (do_pml_in_domain) ? + BoxArray(grid_cba.boxList().intersect(cdomain)) : grid_cba; - // Assuming that refinement ratio is equal in all dimensions - const BoxArray& cba = (do_pml_in_domain) ? - MakeBoxArray(*cgeom, grid_cba_reduced, ncell/ref_ratio[0], do_pml_in_domain, do_pml_Lo, do_pml_Hi) : - MakeBoxArray(*cgeom, grid_cba, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi); +// This needs to be fixed. AMREX_ALWAYS_ASSERT(ref_ratio == ref_ratio[0]); + const int cncells = ncell/ref_ratio[0]; + const int cdelta = delta/ref_ratio[0]; + // Assuming that refinement ratio is equal in all dimensions + const BoxArray& cba = MakeBoxArray(is_single_box_domain, cdomain, *cgeom, grid_cba_reduced, + cncells, do_pml_in_domain, do_pml_Lo, do_pml_Hi); DistributionMapping cdm; if (WarpX::do_similar_dm_pml) { auto ng_sim = amrex::elemwiseMax(amrex::elemwiseMax(nge, ngb), ngf); @@ -781,12 +843,9 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri pml_j_cp[1]->setVal(0.0); pml_j_cp[2]->setVal(0.0); - if (do_pml_in_domain){ - // Note - assuming that the refinement ratio is equal in all dimensions - sigba_cp = std::make_unique(cba, cdm, grid_cba_reduced, cgeom->CellSize(), ncell/ref_ratio[0], delta/ref_ratio[0]); - } else { - sigba_cp = std::make_unique(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta); - } + single_domain_box = is_single_box_domain ? cdomain : Box(); + sigba_cp = std::make_unique(cba, cdm, grid_cba_reduced, cgeom->CellSize(), + cncells, cdelta, single_domain_box); if (WarpX::maxwell_solver_id == MaxwellSolverAlgo::PSATD) { #ifndef WARPX_USE_PSATD @@ -815,9 +874,61 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& gri } BoxArray -PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, +PML::MakeBoxArray (bool is_single_box_domain, const amrex::Box& regular_domain, + const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell, int do_pml_in_domain, const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) +{ + if (is_single_box_domain) { + return MakeBoxArray_single(regular_domain, grid_ba, ncell, do_pml_Lo, do_pml_Hi); + } else { // the union of the regular grids is *not* a single rectangular domain + return MakeBoxArray_multiple(geom, grid_ba, ncell, do_pml_in_domain, do_pml_Lo, do_pml_Hi); + } +} + +BoxArray +PML::MakeBoxArray_single (const amrex::Box& regular_domain, const amrex::BoxArray& grid_ba, + int ncell, const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) +{ + BoxList bl; + for (int i = 0, N = grid_ba.size(); i < N; ++i) { + Box const& b = grid_ba[i]; + for (OrientationIter oit; oit.isValid(); ++oit) { + // In 3d, a Box has 6 faces. This iterates over the 6 faces. + // 3 of them are on the lower side and the others are on the + // higher side. + Orientation ori = oit(); + const int idim = ori.coordDir(); // either 0 or 1 or 2 (i.e., x, y, z-direction) + bool pml_bndry = false; + if (ori.isLow() && do_pml_Lo[idim]) { // This is one of the lower side faces. + pml_bndry = b.smallEnd(idim) == regular_domain.smallEnd(idim); + } else if (ori.isHigh() && do_pml_Hi[idim]) { // This is one of the higher side faces. + pml_bndry = b.bigEnd(idim) == regular_domain.bigEnd(idim); + } + if (pml_bndry) { + Box bbox = amrex::adjCell(b, ori, ncell); + for (int jdim = 0; jdim < idim; ++jdim) { + if (do_pml_Lo[jdim] && + bbox.smallEnd(jdim) == regular_domain.smallEnd(jdim)) { + bbox.growLo(jdim, ncell); + } + if (do_pml_Hi[jdim] && + bbox.bigEnd(jdim) == regular_domain.bigEnd(jdim)) { + bbox.growHi(jdim, ncell); + } + } + bl.push_back(bbox); + } + } + } + + return BoxArray(std::move(bl)); +} + +BoxArray +PML::MakeBoxArray_multiple (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, + int ncell, int do_pml_in_domain, + const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) { Box domain = geom.Domain(); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {