ldz_gradpe.cpp

Go to the documentation of this file.

/*
 * This file is part of Vlasiator.
 * Copyright 2010-2016 Finnish Meteorological Institute
 *
 * For details of usage, see the COPYING file and read the "Rules of the Road"
 * at http://www.physics.helsinki.fi/vlasiator/
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License along
 * with this program; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 */
 
#include "fs_common.h"
#include "ldz_gradpe.hpp"
#include <limits>
 
#ifdef DEBUG_VLASIATOR
   #define DEBUG_FSOLVER
#endif
 
using namespace std;

 

void calculateGradPeTerm(fsgrids::egradpespan egradpes,
                         fsgrids::constmomentsspan moments,
                         fsgrids::constdmomentsspan dmoments,
                         fsgrids::consttechnicalspan technical, const fsgrid::FsStencil& stencil,
                         const auto& gridSpacing, SysBoundary& sysBoundaries) {
#ifdef DEBUG_FSOLVER
   if (stencil.ooo() >= moments.size()) {
      cerr << "Out-of-bounds access in " << __FILE__ << ":" << __LINE__ << endl;
      exit(1);
   }
#endif
   auto& egradpe = egradpes[stencil.ooo()];
   const auto& moment = moments[stencil.ooo()];
   const auto& dmoment = dmoments[stencil.ooo()];
   const auto& tech = technical[stencil.ooo()];
 
   cuint cellSysBoundaryFlag = tech.sysBoundaryFlag;
   cuint cellSysBoundaryLayer = tech.sysBoundaryLayer;
 
   if (cellSysBoundaryFlag == sysboundarytype::DO_NOT_COMPUTE ||
       cellSysBoundaryFlag == sysboundarytype::OUTER_BOUNDARY_PADDING) {
      return;
   }
 
   if ((cellSysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY) && (cellSysBoundaryLayer != 1)) {
      sysBoundaries.getSysBoundary(cellSysBoundaryFlag)->fieldSolverBoundaryCondGradPeElectricField(egradpes, stencil, 0);
      sysBoundaries.getSysBoundary(cellSysBoundaryFlag)->fieldSolverBoundaryCondGradPeElectricField(egradpes, stencil, 1);
      sysBoundaries.getSysBoundary(cellSysBoundaryFlag)->fieldSolverBoundaryCondGradPeElectricField(egradpes, stencil, 2);
   } else {
      if (Parameters::ohmGradPeTerm == 0) {
         cerr << __FILE__ << __LINE__
              << "You shouldn't be in a electron pressure gradient term function if Parameters::ohmGradPeTerm == 0."
              << endl;
      } else if (Parameters::ohmGradPeTerm == 1) {
         auto calculateEdgeGradPeTermComponent = [&egradpe, &moment, &dmoment](size_t i, size_t j, Real spacing) {
            const Real rhoq = moment[fsgrids::moments::RHOQ];
            const Real min = Parameters::hallMinimumRhoq;
            const Real max = std::numeric_limits<Real>::max();
            const Real hallRhoq = std::clamp(rhoq, min, max);
            egradpe[i] = -dmoment[j] / (hallRhoq * spacing);
         };
         calculateEdgeGradPeTermComponent(fsgrids::egradpe::EXGRADPE, fsgrids::dmoments::dPedx, gridSpacing[0]);
         calculateEdgeGradPeTermComponent(fsgrids::egradpe::EYGRADPE, fsgrids::dmoments::dPedy, gridSpacing[1]);
         calculateEdgeGradPeTermComponent(fsgrids::egradpe::EZGRADPE, fsgrids::dmoments::dPedz, gridSpacing[2]);
      } else {
         cerr << __FILE__ << ":" << __LINE__ << "You are welcome to code higher-order Hall term correction terms." << endl;
      }
   }
}
 
void calculateGradPeTermSimple(fsgrids::egradpespan egradpe,
                               fsgrids::egradpespan egradpedt2,
                               fsgrids::momentsspan moments,
                               fsgrids::momentsspan momentsdt2,
                               fsgrids::dmomentsspan dmoments,
                               fsgrids::dmomentsspan dmomentsdt2,
                               fsgrids::technicalspan technical, FieldSolverGrid &fsgrid,
                               SysBoundary& sysBoundaries, cint& RKCase) {
   phiprof::Timer gradPeTimer{"Calculate GradPe term"};
 
   const size_t numCells = fsgrid.getNumCells();
 
   if (not(RKCase == RK_ORDER1 || RKCase == RK_ORDER2_STEP2)) {
      egradpe = egradpedt2;
      moments = momentsdt2;
      dmoments = dmomentsdt2;
   }
 
   phiprof::Timer mpiTimer{"EgradPe field update ghosts MPI", {"MPI"}};
   fsgrid.updateGhostCells(dmoments);
   mpiTimer.stop();
 
   // Calculate GradPe term
   fsgrid.parallel_for([](int timerId) -> phiprof::Timer { return phiprof::Timer{timerId}; },
                       phiprof::initializeTimer("EgradPe compute cells"), technical,
                       [=, &sysBoundaries](const fsgrid::Coordinates &coordinates, const fsgrid::FsStencil& stencil, cuint sysBoundaryFlag, cuint sysBoundaryLayer) {
                          calculateGradPeTerm(egradpe, moments, dmoments, technical, stencil, coordinates.physicalGridSpacing, sysBoundaries);
                       });
 
   gradPeTimer.stop(numCells, "Spatial Cells");
}