CD_ItoSolverImplem.H

Go to the documentation of this file.

/*
 * SPDX-FileCopyrightText: 2021-2026 SINTEF Energy Research
 *
 * SPDX-License-Identifier: GPL-3.0-or-later
 */
 
#ifndef CD_ITOSOLVERIMPLEM_H
#define CD_ITOSOLVERIMPLEM_H
 
// Std includes
#include <chrono>
 
// Chombo includes
#include <EBAlias.H>
#include <PolyGeom.H>
 
// Our includes
#include <CD_ItoSolver.H>
#include <CD_DataOps.H>
#include <CD_Random.H>
#include <CD_NamespaceHeader.H>
 
inline RealVect
ItoSolver::randomGaussian() const
{
  // TLDR: We draw a random number from a Gaussian distribution for each coordinate, and truncate the distribution at m_normalDistributionTruncation.
 
  auto sign = [](const Real& a) -> Real {
    return (a > 0) - (a < 0);
  };
 
  RealVect r = RealVect::Zero;
  for (int i = 0; i < SpaceDim; i++) {
    r[i] = Random::getNormal01();
 
    r[i] = sign(r[i]) * std::min(std::abs(r[i]), m_normalDistributionTruncation);
  }
 
  return r;
}
 
template <ItoSolver::WhichContainer C>
void
ItoSolver::addParticles(ListBox<ItoParticle>& a_inputParticles,
                        const int             a_lvl,
                        const DataIndex&      a_dit,
                        const bool            a_destructive)
{
  CH_TIME("ItoSolver::addParticles");
  if (m_verbosity > 5) {
    pout() << m_name + "::addParticles" << endl;
  }
 
  ParticleContainer<ItoParticle>& particles = m_particleContainers.at(C);
 
  ListBox<ItoParticle>& my_particles = particles[a_lvl][a_dit];
  if (a_destructive) {
    my_particles.addItemsDestructive(a_inputParticles.listItems());
  }
  else {
    my_particles.addItems(a_inputParticles.listItems());
  }
}
 
template <class P, class Ret, Ret (P ::*MemberFunc)() const>
void
ItoSolver::depositParticles(EBAMRCellData& a_phi, ParticleContainer<P>& a_particles) const
{
  CH_TIME("ItoSolver::depositParticles");
  if (m_verbosity > 5) {
    pout() << m_name + "::depositParticles" << endl;
  }
 
  this->depositParticles<P, Ret, MemberFunc>(a_phi, a_particles, m_deposition, m_coarseFineDeposition);
}
 
template <class P, class Ret, Ret (P ::*MemberFunc)() const>
void
ItoSolver::depositParticlesNGP(LevelData<EBCellFAB>&       a_output,
                               const ParticleContainer<P>& a_particles,
                               const int                   a_level) const noexcept
{
  CH_TIME("ItoSolver::depositParticlesNGP");
  if (m_verbosity > 5) {
    pout() << m_name + "::depositParticlesNGP" << endl;
  }
 
  CH_assert(a_level >= 0);
  CH_assert(a_level <= m_amr->getFinestLevel());
 
  const ProblemDomain&     domain = m_amr->getDomains()[a_level];
  const DisjointBoxLayout& dbl    = m_amr->getGrids(a_particles.getRealm())[a_level];
  const DataIterator&      dit    = dbl.dataIterator();
  const EBISLayout&        ebisl  = m_amr->getEBISLayout(a_particles.getRealm(), m_phase)[a_level];
  const Real               dx     = m_amr->getDx()[a_level];
  const RealVect           probLo = m_amr->getProbLo();
 
  CH_assert(a_output.disjointBoxLayout() == dbl);
 
  const int nbox = dit.size();
 
#pragma omp parallel for schedule(runtime)
  for (int mybox = 0; mybox < nbox; mybox++) {
    const DataIndex& din = dit[mybox];
 
    const Box      cellBox = dbl[din];
    const EBISBox& ebisbox = ebisl[din];
 
    EBParticleMesh particleMesh(domain, cellBox, ebisbox, dx * RealVect::Unit, probLo);
 
    EBCellFAB&     output    = a_output[din];
    const List<P>& particles = a_particles[a_level][din].listItems();
 
    particleMesh.deposit<P, Ret, MemberFunc>(output, particles, DepositionType::NGP, 1.0, true);
  }
}
 
template <class P, class Ret, Ret (P ::*MemberFunc)() const>
void
ItoSolver::depositParticles(EBAMRCellData&             a_phi,
                            ParticleContainer<P>&      a_particles,
                            const DepositionType       a_deposition,
                            const CoarseFineDeposition a_coarseFineDeposition) const
{
  CH_TIME("ItoSolver::depositParticles");
  if (m_verbosity > 5) {
    pout() << m_name + "::depositParticles" << endl;
  }
 
  CH_assert(a_phi[0]->nComp() == 1);
  CH_assert(!a_particles.isOrganizedByCell());
 
  // TLDR: First, deposit onto the mesh as usual (as if the EB wasn't there). If the user asks for it, he can redistribute mass in order to
  //       conserve total mass (if that is important). But the corresponding scheme will be O(1) accurate.
  this->depositKappaConservative<P, Ret, MemberFunc>(a_phi, a_particles, a_deposition, a_coarseFineDeposition);
 
  // Redistribution magic.
  this->redistributeAMR(a_phi);
 
  // Average down and interpolate
  m_amr->conservativeAverage(a_phi, m_realm, m_phase);
  m_amr->interpGhost(a_phi, m_realm, m_phase);
}
 
template <class P, class Ret, Ret (P ::*MemberFunc)() const>
void
ItoSolver::depositKappaConservative(EBAMRCellData&             a_phi,
                                    ParticleContainer<P>&      a_particles,
                                    const DepositionType       a_deposition,
                                    const CoarseFineDeposition a_coarseFineDeposition) const
{
  CH_TIME("ItoSolver::depositKappaConservative");
  if (m_verbosity > 5) {
    pout() << m_name + "::depositKappaConservative" << endl;
  }
 
  CH_assert(a_phi[0]->nComp() == 1);
 
  m_amr->depositParticles<P, Ret, MemberFunc>(a_phi,
                                              m_realm,
                                              m_phase,
                                              a_particles,
                                              a_deposition,
                                              a_coarseFineDeposition,
                                              m_forceIrregDepositionNGP);
}
 
#include <CD_NamespaceFooter.H>
 
#endif