Skip to content

Commit

Permalink
Nuclear fusion algorithm modifications (ECP-WarpX#5133)
Browse files Browse the repository at this point in the history 
* using new method for fusion reaction probability.

* weights of colliding particles are reduced locally for fusion.

* simplifying changes to NuclearFusionFunc.H

* updating checksum.

* removing DSMC check lines from ParticleCreationFunc.H

* removed unused variables.

* updaing 3D PB11 checksum.

* alter comment. pushing again to test unexepected CI test failure.

* fixed mistake in checksum update.

* reducing duplicate coding.
  • Loading branch information
@JustinRayAngus
JustinRayAngus authored Aug 18, 2024
1 parent 90e858d commit bb7d243
Show file tree
Hide file tree
Showing 4 changed files with 96 additions and 104 deletions.
46 changes: 23 additions & 23 deletions Regression/Checksum/benchmarks_json/Proton_Boron_Fusion_2D.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,14 +16,14 @@
"particle_momentum_z": 9.279446607709548e-15,
"particle_position_x": 4096178.1664224654,
"particle_position_y": 8192499.7060386725,
"particle_weight": 6.399486212205656e+30
"particle_weight": 6.399499907503547e+30
},
"alpha5": {
"particle_momentum_x": 2.3827075626988002e-14,
"particle_momentum_y": 2.388035811027714e-14,
"particle_momentum_z": 2.4040151761604078e-14,
"particle_position_x": 2457556.8571638423,
"particle_position_y": 4914659.635379322,
"particle_momentum_x": 2.3859593638404282e-14,
"particle_momentum_y": 2.38840514351394e-14,
"particle_momentum_z": 2.3991308560597987e-14,
"particle_position_x": 2457556.8571638414,
"particle_position_y": 4914659.635379324,
"particle_weight": 3.839999999999998e-19
},
"boron2": {
Expand Down Expand Up @@ -51,12 +51,12 @@
"particle_weight": 1.2800000000000004e+37
},
"proton3": {
"particle_momentum_x": 1.6847290893972251e-15,
"particle_momentum_y": 1.6827074502304075e-15,
"particle_momentum_z": 1.6802489646490975e-15,
"particle_position_x": 2457270.6999197667,
"particle_position_y": 4914315.665267942,
"particle_weight": 1.279486212205663e+30
"particle_momentum_x": 1.6842001035542442e-15,
"particle_momentum_y": 1.6822700479972656e-15,
"particle_momentum_z": 1.6797740090740159e-15,
"particle_position_x": 2456930.3827292607,
"particle_position_y": 4913715.277730561,
"particle_weight": 1.2794999075035503e+30
},
"alpha2": {
"particle_momentum_x": 4.140109871048478e-15,
Expand Down Expand Up @@ -91,11 +91,11 @@
"particle_weight": 127.99999999999999
},
"alpha4": {
"particle_momentum_x": 2.3854178157605226e-14,
"particle_momentum_y": 2.3882457581100904e-14,
"particle_momentum_z": 2.403451456378969e-14,
"particle_position_x": 2457367.4582781536,
"particle_position_y": 4915112.044373056,
"particle_momentum_x": 2.3844872070315362e-14,
"particle_momentum_y": 2.3880317964250153e-14,
"particle_momentum_z": 2.4010930856446718e-14,
"particle_position_x": 2457367.4582781526,
"particle_position_y": 4915112.044373058,
"particle_weight": 384.0000000000002
},
"proton2": {
Expand All @@ -107,11 +107,11 @@
"particle_weight": 1.2885512788807322e+28
},
"alpha3": {
"particle_momentum_x": 5.079136896829917e-16,
"particle_momentum_y": 5.075922501147803e-16,
"particle_momentum_z": 4.962151258214859e-16,
"particle_position_x": 52678.192400911765,
"particle_position_y": 105483.59950020742,
"particle_weight": 1.5413633830148085e+27
"particle_momentum_x": 4.758972533179803e-16,
"particle_momentum_y": 4.744688612812961e-16,
"particle_momentum_z": 4.703298089277572e-16,
"particle_position_x": 49191.438625552386,
"particle_position_y": 98893.17267714937,
"particle_weight": 1.500277489351144e+27
}
}
52 changes: 26 additions & 26 deletions Regression/Checksum/benchmarks_json/Proton_Boron_Fusion_3D.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@
"particle_position_x": 4096178.1664224654,
"particle_position_y": 4096499.7060386725,
"particle_position_z": 8191465.586938233,
"particle_weight": 5.1196455431551905e+31
"particle_weight": 5.119677303391259e+31
},
"alpha2": {
"particle_momentum_x": 4.172264972648105e-15,
Expand Down Expand Up @@ -58,13 +58,13 @@
"particle_weight": 1023.9999999999999
},
"proton3": {
"particle_momentum_x": 1.6847282386883186e-15,
"particle_momentum_y": 1.6828065767793222e-15,
"particle_momentum_z": 1.6803456707569493e-15,
"particle_position_x": 2457343.371083716,
"particle_position_y": 2457033.3891170574,
"particle_position_z": 4914529.855222688,
"particle_weight": 1.023645543155193e+31
"particle_momentum_x": 1.6843593058271961e-15,
"particle_momentum_y": 1.6823892851456668e-15,
"particle_momentum_z": 1.6800222030208115e-15,
"particle_position_x": 2457072.7669634065,
"particle_position_y": 2456771.196542574,
"particle_position_z": 4914021.8271464575,
"particle_weight": 1.0236773033912632e+31
},
"proton4": {
"particle_momentum_x": 0.0,
Expand All @@ -76,29 +76,29 @@
"particle_weight": 1.0240000000000003e+38
},
"alpha3": {
"particle_momentum_x": 4.822525301678111e-16,
"particle_momentum_y": 4.78793407068675e-16,
"particle_momentum_z": 4.710411155693663e-16,
"particle_position_x": 50987.442011759704,
"particle_position_y": 48999.674675246955,
"particle_position_z": 101142.57224226737,
"particle_weight": 1.0633705344227063e+28
"particle_momentum_x": 4.529369912319074e-16,
"particle_momentum_y": 4.55762819982174e-16,
"particle_momentum_z": 4.501223434274883e-16,
"particle_position_x": 48732.63289890166,
"particle_position_y": 46262.85664926968,
"particle_position_z": 95593.72900482587,
"particle_weight": 9.680898262134895e+27
},
"alpha5": {
"particle_momentum_x": 2.3856918226484716e-14,
"particle_momentum_y": 2.3827536503955812e-14,
"particle_momentum_z": 2.405465166862308e-14,
"particle_position_x": 2457556.857163843,
"particle_position_y": 2457059.6353793247,
"particle_position_z": 4915847.043341331,
"particle_momentum_x": 2.3875154301762454e-14,
"particle_momentum_y": 2.3859525730749002e-14,
"particle_momentum_z": 2.4029799910201597e-14,
"particle_position_x": 2457556.8571638423,
"particle_position_y": 2457059.6353793237,
"particle_position_z": 4915847.043341332,
"particle_weight": 3.0719999999999984e-18
},
"alpha4": {
"particle_momentum_x": 2.386438069023501e-14,
"particle_momentum_y": 2.3877492448396915e-14,
"particle_momentum_z": 2.4013650859080414e-14,
"particle_position_x": 2457367.458278153,
"particle_position_y": 2457512.0443730573,
"particle_momentum_x": 2.3842314895045605e-14,
"particle_momentum_y": 2.386245392919704e-14,
"particle_momentum_z": 2.4007606987094537e-14,
"particle_position_x": 2457367.4582781545,
"particle_position_y": 2457512.044373057,
"particle_position_z": 4914475.7765130745,
"particle_weight": 3072.000000000002
},
Expand Down
91 changes: 47 additions & 44 deletions Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,11 +145,13 @@ public:
amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux];
amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy];
amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz];
uint64_t* AMREX_RESTRICT idcpu1 = soa_1.m_idcpu;

amrex::ParticleReal * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w];
amrex::ParticleReal * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux];
amrex::ParticleReal * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy];
amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz];
uint64_t* AMREX_RESTRICT idcpu2 = soa_2.m_idcpu;

// Number of macroparticles of each species
const index_type NI1 = I1e - I1s;
Expand All @@ -159,18 +161,9 @@ public:

index_type pair_index = cell_start_pair + coll_idx;

// Because the number of particles of each species is not always equal (NI1 != NI2
// in general), some macroparticles will be paired with multiple macroparticles of the
// other species and we need to decrease their weight accordingly.
// c1 corresponds to the minimum number of times a particle of species 1 will be paired
// with a particle of species 2. Same for c2.
// index_type(1): https://github.com/AMReX-Codes/amrex/pull/3684
const index_type c1 = amrex::max(NI2/NI1, index_type(1));
const index_type c2 = amrex::max(NI1/NI2, index_type(1));

// multiplier ratio to take into account unsampled pairs
const auto multiplier_ratio = static_cast<int>(
m_isSameSpecies ? 2*max_N - 1 : max_N);
m_isSameSpecies ? min_N + max_N - 1 : min_N);

#if (defined WARPX_DIM_RZ)
amrex::ParticleReal * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta];
Expand All @@ -182,50 +175,60 @@ public:
// stride (smaller set size) until we do all collisions (larger set size)
for (index_type k = coll_idx; k < max_N; k += min_N)
{
// c1k : how many times the current particle of species 1 is paired with a particle
// of species 2. Same for c2k.
const index_type c1k = (k%NI1 < max_N%NI1) ? c1 + 1: c1;
const index_type c2k = (k%NI2 < max_N%NI2) ? c2 + 1: c2;

// do not check for collision if a particle's weight was
// reduced to zero from a previous collision
if (idcpu1[ I1[i1] ]==amrex::ParticleIdCpus::Invalid ||
idcpu2[ I2[i2] ]==amrex::ParticleIdCpus::Invalid ) {
p_mask[pair_index] = false;
}
else {

#if (defined WARPX_DIM_RZ)
/* In RZ geometry, macroparticles can collide with other macroparticles
* in the same *cylindrical* cell. For this reason, collisions between macroparticles
* are actually not local in space. In this case, the underlying assumption is that
* particles within the same cylindrical cell represent a cylindrically-symmetry
* momentum distribution function. Therefore, here, we temporarily rotate the
* momentum of one of the macroparticles in agreement with this cylindrical symmetry.
* (This is technically only valid if we use only the m=0 azimuthal mode in the simulation;
* there is a corresponding assert statement at initialization.) */
amrex::ParticleReal const theta = theta2[I2[i2]]-theta1[I1[i1]];
amrex::ParticleReal const u1xbuf = u1x[I1[i1]];
u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta);
u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta);
/* In RZ geometry, macroparticles can collide with other macroparticles
* in the same *cylindrical* cell. For this reason, collisions between macroparticles
* are actually not local in space. In this case, the underlying assumption is that
* particles within the same cylindrical cell represent a cylindrically-symmetry
* momentum distribution function. Therefore, here, we temporarily rotate the
* momentum of one of the macroparticles in agreement with this cylindrical symmetry.
* (This is technically only valid if we use only the m=0 azimuthal mode in the simulation;
* there is a corresponding assert statement at initialization.) */
amrex::ParticleReal const theta = theta2[I2[i2]]-theta1[I1[i1]];
amrex::ParticleReal const u1xbuf = u1x[I1[i1]];
u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta);
u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta);
#endif

SingleNuclearFusionEvent(
u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ],
u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ],
m1, m2, w1[ I1[i1] ]/c1k, w2[ I2[i2] ]/c2k,
dt, dV, static_cast<int>(pair_index), p_mask, p_pair_reaction_weight,
m_fusion_multiplier, multiplier_ratio,
m_probability_threshold,
m_probability_target_value,
m_fusion_type, engine);
SingleNuclearFusionEvent(
u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ],
u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ],
m1, m2, w1[ I1[i1] ], w2[ I2[i2] ],
dt, dV, static_cast<int>(pair_index), p_mask, p_pair_reaction_weight,
m_fusion_multiplier, multiplier_ratio,
m_probability_threshold,
m_probability_target_value,
m_fusion_type, engine);

#if (defined WARPX_DIM_RZ)
amrex::ParticleReal const u1xbuf_new = u1x[I1[i1]];
u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta);
u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta);
amrex::ParticleReal const u1xbuf_new = u1x[I1[i1]];
u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta);
u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta);
#endif

// Remove pair reaction weight from the colliding particles' weights
if (p_mask[pair_index]) {
BinaryCollisionUtils::remove_weight_from_colliding_particle(
w1[ I1[i1] ], idcpu1[ I1[i1] ], p_pair_reaction_weight[pair_index]);
BinaryCollisionUtils::remove_weight_from_colliding_particle(
w2[ I2[i2] ], idcpu2[ I2[i2] ], p_pair_reaction_weight[pair_index]);
}

}

p_pair_indices_1[pair_index] = I1[i1];
p_pair_indices_2[pair_index] = I2[i2];
if (max_N == NI1) {
i1 += min_N;
}
if (max_N == NI2) {
i2 += min_N;
}
if (max_N == NI1) { i1 += min_N; }
if (max_N == NI2) { i2 += min_N; }
pair_index += min_N;
}
}
Expand Down
11 changes: 0 additions & 11 deletions Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,11 +131,6 @@ public:
const auto soa_1 = ptile1.getParticleTileData();
const auto soa_2 = ptile2.getParticleTileData();

amrex::ParticleReal* AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
amrex::ParticleReal* AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w];
uint64_t* AMREX_RESTRICT idcpu1 = soa_1.m_idcpu;
uint64_t* AMREX_RESTRICT idcpu2 = soa_2.m_idcpu;

// Create necessary GPU vectors, that will be used in the kernel below
amrex::Vector<SoaData_type> soa_products;
for (int i = 0; i < m_num_product_species; i++)
Expand Down Expand Up @@ -197,12 +192,6 @@ public:
}
}

// Remove p_pair_reaction_weight[i] from the colliding particles' weights
BinaryCollisionUtils::remove_weight_from_colliding_particle(
w1[p_pair_indices_1[i]], idcpu1[p_pair_indices_1[i]], p_pair_reaction_weight[i]);
BinaryCollisionUtils::remove_weight_from_colliding_particle(
w2[p_pair_indices_2[i]], idcpu2[p_pair_indices_2[i]], p_pair_reaction_weight[i]);

// Initialize the product particles' momentum, using a function depending on the
// specific collision type
if (t_collision_type == CollisionType::ProtonBoronToAlphasFusion)
Expand Down

0 comments on commit bb7d243

Please sign in to comment.