lpt_props.h

Go to the documentation of this file.

// Copyright (C) 2024 The m-AIA AUTHORS
//
// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
//
// SPDX-License-Identifier: LGPL-3.0-only
 
#include <cmath>
#include "UTIL/functions.h"
#include "lpt.h"
 
template <MInt nDim>
void LPT<nDim>::readModelProps() {
  TRACE();
 
  m_log << "reading model properties...";
 
  m_restart = false;
  m_restart = Context::getSolverProperty<MBool>("restartFile", solverId(), AT_, &m_restart);
 
  m_restartTimeStep = Context::getSolverProperty<MInt>("restartTimeStep", solverId(), AT_);
 
  m_maxNoParticles = 1E6;
  m_maxNoParticles = Context::getSolverProperty<MInt>("maxNoParticles", solverId(), AT_, &m_maxNoParticles);
 
  m_exchangeBufferSize =
      Context::getSolverProperty<MInt>("particleExchangeBufferSize", solverId(), AT_, &m_exchangeBufferSize);
 
 
  m_sizeLimit = 1E-12;
  m_sizeLimit = Context::getSolverProperty<MFloat>("particleSizeLimit", solverId(), AT_, &m_sizeLimit);
 
  m_spawnParticles = false;
  m_spawnParticles = Context::getSolverProperty<MBool>("spawnParticles", solverId(), AT_, &m_spawnParticles);
 
  if(m_spawnParticles) {
    m_log << "Particle spawning has been activated." << std::endl;
  }
 
  m_activePrimaryBUp = false;
  m_activePrimaryBUp = Context::getSolverProperty<MBool>("sprayPrimaryBUp", solverId(), AT_, &m_activePrimaryBUp);
  m_activeSecondaryBUp = false;
  m_activeSecondaryBUp = Context::getSolverProperty<MBool>("spraySecondaryBUp", solverId(), AT_, &m_activeSecondaryBUp);
  m_ellipsoids = false;
  m_ellipsoids = Context::getSolverProperty<MBool>("particleIncludeEllipsoids", solverId(), AT_, &m_ellipsoids);
 
  m_collisions = Context::getSolverProperty<MInt>("particleCollisions", solverId(), AT_, &m_collisions);
 
  if(m_collisions > 0) {
    m_collisionModel = std::unique_ptr<ParticleCollision<nDim>>(
        new ParticleCollision<nDim>(this, m_collisions, solverId(), domainId(), noDomains(), mpiComm()));
  }
 
  m_wallCollisions = Context::getSolverProperty<MBool>("particleWallCollisions", solverId(), AT_, &m_wallCollisions);
 
  m_initializationMethod =
      Context::getSolverProperty<MInt>("particleInitializationMethod", solverId(), AT_, &m_initializationMethod);
 
  m_dragModelType = Context::getSolverProperty<MInt>("particleDrag", solverId(), AT_, &m_dragModelType);
 
  m_liftModelType = Context::getSolverProperty<MInt>("particleLift", solverId(), AT_, &m_liftModelType);
 
  m_torqueModelType = Context::getSolverProperty<MInt>("particleTorque", solverId(), AT_, &m_torqueModelType);
 
  m_motionEquationType = Context::getSolverProperty<MInt>("motionEquation", solverId(), AT_, &m_motionEquationType);
 
  m_respawn = Context::getSolverProperty<MBool>("particleRespawn", solverId(), AT_, &m_respawn);
 
  if(m_respawn) {
    m_respawnPlane = Context::getSolverProperty<MFloat>("particleRespawnPlane", solverId(), AT_);
 
    // create a vector of cells which are used for the "respawning" of particles
    MFloat halfLength = NAN;
    for(MInt i = 0; i < noInternalCells(); i++) {
      if(!c_isLeafCell(i)) continue;
      MInt level = c_level(i);
      halfLength = c_cellLengthAtLevel(level + 1);
      if(fabs(c_coordinate(i, 0) - m_respawnPlane) < halfLength) {
        auto weightedTimes = static_cast<MInt>(pow(pow(2, nDim), (maxRefinementLevel() - level)));
        for(MInt count = 0; count < weightedTimes; ++count) {
          m_respawnCells.push_back(i);
        }
      }
    }
    m_log << m_respawnCells.size() << " respawn cells found in this solver" << std::endl;
  }
 
  m_xCutOff = Context::getSolverProperty<MFloat>("particlePartOffset", solverId(), AT_, &m_xCutOff);
 
  m_momentumCoupling = false;
  m_momentumCoupling =
      Context::getSolverProperty<MBool>("particleMomentumCoupling", solverId(), AT_, &m_momentumCoupling);
  if(m_momentumCoupling) {
    m_log << "Momentum coupling has been activated." << std::endl;
  }
 
  m_heatCoupling = false;
  m_heatCoupling = Context::getSolverProperty<MBool>("particleHeatCoupling", solverId(), AT_, &m_heatCoupling);
 
  if(m_heatCoupling) {
    m_log << "Heat coupling has been activated." << std::endl;
  }
 
  m_massCoupling = false;
  m_massCoupling = Context::getSolverProperty<MBool>("particleMassCoupling", solverId(), AT_, &m_massCoupling);
  if(m_massCoupling) {
    m_log << "Mass coupling has been activated." << std::endl;
  }
 
  m_evaporation = false;
  m_evaporation = Context::getSolverProperty<MBool>("particleEvaporation", solverId(), AT_, &m_evaporation);
 
  if(m_evaporation) {
    m_log << "Evaporation modeling has been activated." << std::endl;
  }
 
  m_maxNoBndryCells = Context::getSolverProperty<MInt>("maxNoBndryCells", solverId(), AT_);
 
 
  MInt range[2] = {2, 0};
  for(MInt i = 0; i < 2; i++) {
    range[i] = Context::getSolverProperty<MInt>("particleAdapRange", solverId(), AT_, &range[i], i);
  }
 
  mAlloc(m_bandWidth, maxRefinementLevel() + 1, "m_bandWidth", 0, AT_);
  m_bandWidth[maxRefinementLevel() - 1] = range[0];
  for(MInt i = maxRefinementLevel() - 2; i >= 0; i--) {
    m_bandWidth[i] = (m_bandWidth[i + 1] / 2) + 1 + range[1];
  }
 
  m_innerBound = Context::getSolverProperty<MInt>("innerBound", solverId(), AT_, &m_innerBound);
 
  MFloat referenceTemperature = 273.15;
  referenceTemperature =
      Context::getSolverProperty<MFloat>("referenceTemperature", m_solverId, AT_, &referenceTemperature);
 
  m_sutherlandConstant = 110.4;
  m_sutherlandConstant =
      Context::getSolverProperty<MFloat>("sutherlandConstant", solverId(), AT_, &m_sutherlandConstant);
 
  m_sutherlandConstant /= referenceTemperature;
  m_sutherlandPlusOne = m_sutherlandConstant + F1;
 
  m_cfl = 1.0;
  m_cfl = Context::getSolverProperty<MFloat>("cfl", solverId(), AT_, &m_cfl);
 
  m_timeStepComputationInterval = m_restart ? -1 : 0;
  m_timeStepComputationInterval =
      Context::getSolverProperty<MInt>("timeStepComputationInterval", solverId(), AT_, &m_timeStepComputationInterval);
 
 
  m_nonBlockingComm = Context::getSolverProperty<MBool>("nonBlockingComm", solverId(), AT_, &m_nonBlockingComm);
 
  if(m_nonBlockingComm) {
    m_nonBlockingStage = 0;
    if(m_ellipsoids) {
      cerr0 << "LPT-Warning: nonBlocking currently only works with either spherical or ellipsoidal particles "
            << "--> disabling exchange of spherical particles" << std::endl;
    }
  }
 
  m_weightBaseCell = 0.0;
  m_weightBaseCell = Context::getSolverProperty<MFloat>("weightBaseCell", solverId(), AT_, &m_weightBaseCell);
  m_weightLeafCell = 0.0;
  m_weightLeafCell = Context::getSolverProperty<MFloat>("weightLeafCell", solverId(), AT_, &m_weightLeafCell);
  m_weightParticleCell = 0.1;
  m_weightParticleCell =
      Context::getSolverProperty<MFloat>("weightParticleCell", solverId(), AT_, &m_weightParticleCell);
  m_weightParticle = 1.0;
  m_weightParticle = Context::getSolverProperty<MFloat>("weightParticle", solverId(), AT_, &m_weightParticle);
  m_weightSpawnCell = 0.0;
  m_weightSpawnCell = Context::getSolverProperty<MFloat>("weightSpawnCell", solverId(), AT_, &m_weightSpawnCell);
  m_weightMulitSolverFactor = 1.0;
  m_weightMulitSolverFactor =
      Context::getSolverProperty<MFloat>("weightMulitSolverFactor", solverId(), AT_, &m_weightMulitSolverFactor);
  m_limitWeights = false;
  m_limitWeights = Context::getSolverProperty<MBool>("limitDLBWeights", solverId(), AT_, &m_limitWeights);
 
  m_weightSourceCells = false;
  m_weightSourceCells =
      Context::getSolverProperty<MBool>("weightLPTSourceCells", solverId(), AT_, &m_weightSourceCells);
 
  m_domainIdOutput = false;
  m_domainIdOutput = Context::getSolverProperty<MBool>("domainIdOutput", m_solverId, AT_, &m_domainIdOutput);
 
  // LPT sleep time in milli-seconds
  m_sleepLPT = -1;
  m_sleepLPT = Context::getSolverProperty<MInt>("sleepTimeLPT", m_solverId, AT_, &m_sleepLPT);
 
  m_skipLPT = false;
  m_skipLPT = Context::getSolverProperty<MBool>("skipLPT", m_solverId, AT_, &m_skipLPT);
 
  m_log << "finished" << std::endl;
}
 
template <MInt nDim>
void LPT<nDim>::readSpawnProps() {
  m_spawnSeed = 5489U;
  if(Context::propertyExists("spawnSeed", solverId())) {
    m_spawnSeed = (MLong)Context::getSolverProperty<MInt>("spawnSeed", solverId(), AT_);
  }
  // initialise the psydo-randon number generator with the given seed
  // this is done identically on all ranks and for all random numbers!
  m_PRNGSpawn.seed(m_spawnSeed);
  m_PRNGRespawn.seed(m_spawnSeed);
 
  m_spawnDiameter = 0.00001;
  m_spawnDiameter = Context::getSolverProperty<MFloat>("spawnDiameter", solverId(), AT_, &m_spawnDiameter);
 
  m_spawnVelocity = Context::getSolverProperty<MFloat>("spawnParticlesInitVelo", solverId(), AT_, &m_spawnVelocity);
 
  m_spawnParticlesCount =
      Context::getSolverProperty<MInt>("spawnParticlesCount", solverId(), AT_, &m_spawnParticlesCount);
 
  m_spawnDistSigmaCoeff =
      Context::getSolverProperty<MFloat>("spawnDistSigmaCoeff", solverId(), AT_, &m_spawnDistSigmaCoeff);
 
  m_spawnParticlesConeAngle =
      Context::getSolverProperty<MFloat>("spawnParticlesConeAngle", solverId(), AT_, &m_spawnParticlesConeAngle);
  for(MInt i = 0; i < nDim; i++) {
    m_spawnDir[i] = 0;
    if(i == 0) m_spawnDir[i] = 1.0;
  }
 
  if(Context::propertyExists("spawnDir", solverId())) {
    if(Context::propertyLength("spawnDir", solverId()) != nDim) {
      TERMM(1, "Need to give a Coordinate for every dimension");
    }
 
    for(MInt i = 0; i < nDim; i++) {
      m_spawnDir[i] = Context::getSolverProperty<MFloat>("spawnDir", solverId(), AT_, i);
    }
  }
 
  m_engineSetup = false;
  m_engineSetup = Context::getSolverProperty<MBool>("engineSetup", solverId(), AT_, &m_engineSetup);
}
 
template <MInt nDim>
void LPT<nDim>::readEllipsoidProps() {
  m_ellipsoidRandomOrientation = Context::getSolverProperty<MInt>("particleEllipsoidRandomOrientation", solverId(), AT_,
                                                                  &m_ellipsoidRandomOrientation);
  m_ellipsoidRandomOrientationSeed = 4865U;
  if(Context::propertyExists("particleEllipsoidRandomOrientationSeed", solverId())) {
    m_ellipsoidRandomOrientationSeed =
        (MLong)Context::getSolverProperty<MInt>("particleEllipsoidRandomOrientationSeed", solverId(), AT_);
  }
}
 
template <MInt nDim>
void LPT<nDim>::readMomentumCouplingProps() {
  m_couplingRedist = false;
  m_couplingRedist = Context::getSolverProperty<MBool>("particleCouplingRedist", solverId(), AT_, &m_couplingRedist);
 
  if(m_couplingRedist) {
    m_noRedistLayer = Context::getSolverProperty<MInt>("particleNoRedistLayer", solverId(), AT_, &m_noRedistLayer);
 
    if(m_noRedistLayer > grid().noHaloLayers() && noDomains() > 1) {
      mTerm(-1, AT_, "ERROR: Number of halo layers needs to be higher than the redistribution area!");
    }
 
    m_log << "Momentum redistribution is active with a redistribution area of " << m_noRedistLayer << std::endl;
  }
}
 
template <MInt nDim>
void LPT<nDim>::findSpawnCellId() {
  const MBool redetermination = m_spawnCellId >= 0;
  const MInt lastSpawnCellId = m_spawnCellId;
  const MInt lastSpawnDomainId = m_spawnDomainId;
 
  m_spawnCellId = -1;
  if(Context::propertyExists("spawnCoordinates", solverId())) {
    if(Context::propertyLength("spawnCoordinates", solverId()) != nDim) {
      TERMM(1, "Need to give a Coordinate for every dimension");
    }
 
    for(MInt i = 0; i < nDim; i++) {
      m_spawnCoord[i] = Context::getSolverProperty<MFloat>("spawnCoordinates", solverId(), AT_, i);
    }
 
    m_spawnCellId = grid().findContainingLeafCell(&m_spawnCoord[0]);
 
    // If we are running parallel we need to find the cell which is
    // closest to the given coordinates on all domains.
    if(noDomains() > 1) {
      MInt spawnDomain = domainId();
      MInt minDomain = -1;
      if(m_spawnCellId < 0 || a_isHalo(m_spawnCellId)) {
        spawnDomain = std::numeric_limits<MInt>::max();
      }
 
      // If more than one domain has cells with the same distance only
      // spawn on the domain with the lowest Id
      MPI_Allreduce(&spawnDomain, &minDomain, 1, MPI_INT, MPI_MIN, mpiComm(), AT_, "spawnDomain", "minDomain");
 
      if(minDomain != domainId()) {
        m_spawnCellId = -1;
      }
 
      m_spawnDomainId = minDomain;
 
      m_log << " LPT selected spawn-cell " << m_spawnCellId << " on domain " << m_spawnDomainId << std::endl;
      if(domainId() == m_spawnDomainId) {
        std::cerr << " LPT selected spawn-cell " << m_spawnCellId << " on domain " << m_spawnDomainId
                  << " with coordinates " << c_coordinate(m_spawnCellId, 0) << " " << c_coordinate(m_spawnCellId, 1)
                  << " " << c_coordinate(m_spawnCellId, nDim - 1) << std::endl;
      }
    }
 
    if(redetermination) {
      std::cerr << " LPT spawn-cellId was redetermined";
      if(m_spawnCellId < 0) {
        std::cerr << " and spawnDomain changed from " << domainId() << " to " << m_spawnDomainId;
        std::cerr << " left-over mass was : " << m_particleResiduum << std::endl;
      } else {
        std::cerr << std::endl;
      }
    }
    if(lastSpawnCellId < 0 && m_spawnCellId > 0 && lastSpawnDomainId > 0) {
      std::cerr << "Revieved left-over mass is : " << m_particleResiduum << std::endl;
    }
    if((m_spawnDomainId + 1 > noDomains() && noDomains() > 1) || (m_spawnCellId < 0 && noDomains() == 1)) {
      TERMM(-1, "Invalid spawn coordinates, values are probably outside the fluid-domain/box?");
    }
  }
}