[WIP] 1D CSR Model #627

xbojanic · 2024-05-31T22:38:08Z

This PR implements a 1D CSR model. Currently, this is the ultrarelativistic steady-state wakefield model (eg, eq 19 of
E. L. Saldin et al, NIMA 398, p. 373-394 (1997), DOI:10.1016/S0168-9002(97)00822-X).

This includes:

wake function calculation
beam profile convolution
apply the force "kick" to the particle momentum.
benchmark example

Close #452

for more information, see https://pre-commit.ci

examples/fodo/Wakefield_Functions_FODO_copy.py

for more information, see https://pre-commit.ci

cemitch99

Congratulations on posting your first ImpactX pull request! The changes to the examples (your local tests) probably don't need to appear in this request. We can add a new test that includes CSR, once the wakefield kick is fully implemented.

for more information, see https://pre-commit.ci

Not part of this PR.

…es into it. Created 2 new files for the wake force momentum kick.

…ge_Deposit_Step Created the directory wakefields for all new wakefield scripts.

for more information, see https://pre-commit.ci

…ge_Deposit_Step Deleted copy files and added paths to new wakefields directory for compiling

for more information, see https://pre-commit.ci

src/particles/wakefields/WakeConvolution.H

src/particles/wakefields/ChargeBinning.H

Co-authored-by: Chad Mitchell <46825199+cemitch99@users.noreply.github.com>

Merge remote-tracking branch 'xbojanic/Charge_Deposit_Step' into Charge_Deposit_Step

for more information, see https://pre-commit.ci

cemitch99 · 2024-06-11T23:11:39Z

src/ImpactX.cpp

+                        /*
+                        else if constexpr (std::is_same_v<std::decay_t<decltype(element)>, ExactSbend>) //Currently internal calculation for m_rc
+                        {
+                            R = element.m_rc;


This case (the ExactSBend) will cause a problem because m_rc is not a member variable. Instead, rc is computed internally from other parameters. For now, we could add a comment that this case will be addressed later. We can make a simple change to ExactSBend to address this.

…est.

Merge remote-tracking branch 'xbojanic/Charge_Deposit_Step' into Charge_Deposit_Step

for more information, see https://pre-commit.ci

Modernize FFT build.

Guard FFTW logic and throw runtime errors if used but not built.

…it_Step

for more information, see https://pre-commit.ci

src/ImpactX.cpp

for more information, see https://pre-commit.ci

src/ImpactX.cpp

+                    /*bool element_has_csr = false; // Updates to true for example with bend element
+                    double R = 0.0; // Updates for bend element rc
+
+                    //Define lambda function inside of std::visit
+                    std::visit([&R, &element_has_csr](auto &&element)
+                    {
+                        if constexpr (std::is_same_v<std::decay_t<decltype(element)>, Sbend>)
+                        {
+                            R = element.m_rc;
+                            std::cout << "My radius of curvature is:" << R << std::endl;
+                            element_has_csr = true;
+                        }
+                        else if constexpr (std::is_same_v<std::decay_t<decltype(element)>, CFbend>)
+                        {
+                            R = element.m_rc;
+                            std::cout << "My radius of curvature is:" << R << std::endl;
+                            element_has_csr = true;
+                        }
+
+                        else if constexpr (std::is_same_v<std::decay_t<decltype(element)>, ExactSbend>) //Currently internal calculation for m_rc
+                        {
+                            R = element.m_rc;
+                            std::cout << "My radius of curvature is:" << R << std::endl;
+                            element_has_csr = true;
+                        }
+
+                    }, element_variant);
+
+                    Enter loop if lattice has bend element
+                    if (element_has_csr)
+                    {
+                        // Measure beam size, extract the min, max of particle positions
+                        auto const [x_min, y_min, t_min, x_max, y_max, t_max] =
+                            amr_data->m_particle_container->MinAndMaxPositions();
+
+                        using amrex::Real;
+
+                        // Set parameters for charge deposition
+                        bool is_unity_particle_weight = true;
+                        bool GetNumberDensity = true;
+
+                        int num_bins = 100;  // Set resolution
+                        Real bin_min = t_min;
+                        Real bin_max = t_max;
+                        Real bin_size = (bin_max - bin_min) / num_bins;
+
+                        // Allocate memory for the charge profile
+                        Real* dptr_data = new Real[num_bins]();
+                        auto& particle_container = *(amr_data->m_particle_container);
+
+                        //Call charge deposition function
+                        DepositCharge1D(particle_container, dptr_data, num_bins, bin_min, bin_size, is_unity_particle_weight);
+
+                        // Call charge density derivative function
+                        std::vector<double> charge_distribution(dptr_data, dptr_data + num_bins);
+                        std::vector<double> slopes(num_bins - 1);
+                        DerivativeCharge1D(charge_distribution.data(), slopes.data(), num_bins, bin_size, GetNumberDensity); //Use number derivatives for convolution with CSR
+
+                        // Call wake function
+
+                        // Read in external variable bunch_charge
+                        std::cout << "My beam charge is:" << bunch_charge << std::endl;
+
+                        std::vector<double> wake_function(num_bins);
+                        for (int i = 0; i < num_bins; ++i)
+                        {
+                            double s = bin_min + i * bin_size;
+                            wake_function[i] = w_l_csr(s, R, bunch_charge);
+                        }
+
+                        // Call convolution function
+                        std::vector<double> convoluted_wakefield(2 * num_bins - 1);
+                        convolve_fft(slopes.data(), wake_function.data(), slopes.size(), wake_function.size(), bin_size, convoluted_wakefield.data(), 1);
+
+                        //Check convolution
+                        std::cout << "Convoluted wakefield: ";
+                        std::ofstream outfile("convoluted_wakefield.txt");
+                        for (int i = 0; i < convoluted_wakefield.size(); ++i)
+                        {
+                            std::cout << convoluted_wakefield[i] << " ";
+                            outfile << convoluted_wakefield[i] << std::endl;
+                        }
+                        std::cout << std::endl;
+                        outfile.close();
+                        delete[] dptr_data;
+
+                        // Kick particles with wake
+                        impactx::wakepush::WakePush(particle_container, convoluted_wakefield, bin_size);
+                    }*/
