File diffusion.h

File List > code_source > cosmolattice > include > CosmoInterface > definitions > defectsmodule > diffusion.h

Go to the documentation of this file

#ifndef COSMOINTERFACE_DEFINITIONS_DEFECTSMODULE_DIFFUSION_H
#define COSMOINTERFACE_DEFINITIONS_DEFECTSMODULE_DIFFUSION_H

/* This file is part of CosmoLattice, available at www.cosmolattice.net .
 *   Copyright Daniel G. Figueroa, Adrien Florio, Francisco Torrenti and Wessel Valkenburg.
 *   Released under the MIT license, see LICENSE.md. */

// File info: Main contributor(s): Jorge Baeza-Ballesteros,  Year: 2026


#include "TempLat/util/exception.h"
#include "CosmoInterface/evolvers/evolver.h"
#include "CosmoInterface/evolvers/rk2nstoragefields.h"

namespace TempLat
{

    MakeException(DiffusionEvolverInvalid);

    class Diffusion
    {
    public:
        // Put public methods here. These should change very little over time.

        Diffusion() = delete;

        template<typename Model, typename T>
        static void diffuse(Model &model, RunParameters<T> &rPar, FilesManager<Model::NDim> &filesManager, ExtraFields<Model> extraFields)
        {

            bool iAmRoot  = model.getToolBox()->amIRoot();
            typename Model::FloatType auxg;

            if constexpr (Model::DefectsModel && Model::Ns == 1) {
                auxg = model.g;
                model.g = 0;
            }

            if (!RK2NStorageParameters<T>::isRK2n(rPar.diffType)) throw(DiffusionEvolverInvalid("You need to choose a RK integrator to diffuse initial conditions for defect simulations. Abort."));
            if (iAmRoot) say << "Diffusing initial defect conditions";

            auto diffuser = RK2NStorageFields(model, rPar.diffType);

            std::optional<DefectsMeasurer<T>> measurerS;
            std::optional<ComplexScalarMeasurer<T>> measurerCS;
            std::optional<U1Measurer<T>> measurerU1;

            if constexpr (Model::Ns > 0)  measurerS.emplace(model, filesManager, rPar, rPar.appendMode, "_diffusion", true);
            if constexpr (Model::NCs > 0) measurerCS.emplace(model, filesManager, rPar, rPar.appendMode, "_diffusion", false);
            if constexpr (Model::NU1 > 0) measurerU1.emplace(model, filesManager, rPar, rPar.appendMode, "_diffusion", false, false);
            EnergiesMeasurer<T> measurerE(model, filesManager, rPar, rPar.appendMode, "_diffusion", false);
            EnergySnapshotsMeasurer<Model> snapshots(model, rPar, filesManager, rPar.energySnapshotMeasDiffusion, "_diffusion");

            model.dt = rPar.dtdiff;
            T tdiff = 0.;
            ptrdiff_t outputFreqDiffusion = static_cast<int>(round(rPar.tOutFreqDiff / rPar.dtdiff));
            ptrdiff_t outputRareFreqDiffusion = static_cast<int>(round(rPar.tOutRareFreqDiff / rPar.dtdiff));

            for (int i = 0; tdiff < rPar.tmaxdiff + rPar.dtdiff * T(0.5); i++) {
                if ( i % outputFreqDiffusion == 0) {
                    ForEachConjugateMomenta(Model, fld, n, model.getField(fld)(n) = Kernels::get(fld, model, n, KernelsTypes::Diffusion<Model>()) ); //Sets the conjugate momenta fields, which are used as auxiliary variables for the evolution, to zero.
                    if (iAmRoot) say << "Diffusion step " << i << " completed. Current diffusion time " << tdiff << "\n";
                    if constexpr (Model::Ns > 0) measurerS->measure(model, tdiff);
                    if constexpr (Model::NCs > 0) measurerCS->measureStandard(model, tdiff);
                    if constexpr (Model::NU1 > 0) measurerU1->measureStandard(model, tdiff);
                    measurerE.measure(model, tdiff, i == 0);
                }
                if ( i % outputRareFreqDiffusion == 0) snapshots.measure(model, tdiff);

                diffuser.evolve(model, tdiff, KernelsTypes::Diffusion<Model>());
                tdiff += rPar.dtdiff;
            }

            ForLoop(n, 0, Model::Ns - 1, model.piS(n) = 0.; );
            ForLoop(n, 0, Model::NCs - 1, model.piCS(n) = Complexify(0.,0.); );
            ForLoop(n, 0, Model::NU1 - 1, ForLoop(i, 1, Model::NDim, model.piU1(n)(i) = 0.;) );

            model.dt = rPar.dt;

            if constexpr (Model::DefectsModel && Model::Ns == 1) {
                model.g = auxg;
            }
        }

    };


} // namespace TempLat

#endif