ldz_main.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 "ldz_electric_field.hpp"
#include "ldz_magnetic_field.hpp"
#include "ldz_hall.hpp"
#include "ldz_gradpe.hpp"
#include "ldz_volume.hpp"
#include "fs_common.h"
#include "derivatives.hpp"
#include "fs_limiters.h"
#include "mpiconversion.h"
#include "../fieldtracing/fieldtracing.h"
#include "../logger.h"
extern Logger logFile;

bool propagateFields(fsgrids::perbspan perb,
                     fsgrids::perbspan perbdt2,
                     fsgrids::efieldspan e,
                     fsgrids::efieldspan edt2,
                     fsgrids::ehallspan ehall,
                     fsgrids::egradpespan egradpe,
                     fsgrids::egradpespan egradpedt2,
                     fsgrids::momentsspan moments,
                     fsgrids::momentsspan momentsdt2,
                     fsgrids::dperbspan dperb,
                     fsgrids::dmomentsspan dmoments,
                     fsgrids::dmomentsspan dmomentsdt2,
                     fsgrids::bgbspan bgb,
                     fsgrids::volspan vol,
                     fsgrids::technicalspan technical,
                     FieldSolverGrid &fsgrid,
                     SysBoundary& sysBoundaries,
                     creal& dt,
                     cuint subcycles) {
 
   if (subcycles == 0) {
      cerr << "Field solver subcycles cannot be 0." << endl;
      exit(1);
   }
 
   const auto* localSize = &fsgrid.getLocalSize()[0];
 
   fsgrid.parallel_for([](int timerId) -> phiprof::Timer { return phiprof::Timer{timerId}; },
                       phiprof::initializeTimer("Initialize technical.maxFsDt"), technical,
                       [=](const fsgrid::Coordinates &coordinates, const fsgrid::FsStencil& stencil, cuint sysBoundaryFlag, cuint sysBoundaryLayer) {
                          technical[stencil.ooo()].maxFsDt = std::numeric_limits<Real>::max();
                       });
 
   if (subcycles == 1) {
#ifdef FS_1ST_ORDER_TIME
      propagateMagneticFieldSimple(perb, perbdt2, bgb, e, edt2, technical, fsgrid, sysBoundaries, dt, RK_ORDER1);
      calculateDerivativesSimple(perb, moments, dperb, dmoments, technical, fsgrid, true /*doMoments*/);
      if (P::ohmGradPeTerm > 0) {
         calculateGradPeTermSimple(egradpe, egradpedt2, moments, momentsdt2, dmoments, dmomentsdt2, technical, fsgrid, sysBoundaries, RK_ORDER1);
      }
      if (P::ohmHallTerm > 0) {
         calculateHallTermSimple(
            perb,
            perbdt2,
            ehall,
            moments,
            momentsdt2,
            dperb,
            dmoments,
            dmomentsdt2,
            bgb,
            technical,
            fsgrid,
            sysBoundaries,
            RK_ORDER1,
            true // communicateMomentsDerivatives
         );
      }
      calculateUpwindedElectricFieldSimple(
         perb,
         perbdt2,
         e,
         edt2,
         ehall,
         egradpe,
         egradpedt2,
         moments,
         momentsdt2,
         dperb,
         dmoments,
         dmomentsdt2,
         bgb,
         technical,
         fsgrid,
         sysBoundaries,
         RK_ORDER1,
         true // communicateEGradPeOrMomentsDerivatives
      );
#else
      propagateMagneticFieldSimple(perb, perbdt2, bgb, e, edt2, technical, fsgrid, sysBoundaries, dt, RK_ORDER2_STEP1);
      calculateDerivativesSimple(perbdt2, momentsdt2, dperb, dmomentsdt2, technical, fsgrid, true /*doMoments*/);
      if (P::ohmGradPeTerm > 0) {
         calculateGradPeTermSimple(egradpe, egradpedt2, moments, momentsdt2, dmoments, dmomentsdt2, technical, fsgrid, sysBoundaries, RK_ORDER2_STEP1);
      }
      if (P::ohmHallTerm > 0) {
         calculateHallTermSimple(
            perb,
            perbdt2,
            ehall,
            moments,
            momentsdt2,
            dperb,
            dmoments,
            dmomentsdt2,
            bgb,
            technical,
            fsgrid,
            sysBoundaries,
            RK_ORDER2_STEP1,
            true // communicateMomentsDerivatives
         );
      }
      calculateUpwindedElectricFieldSimple(
         perb,
         perbdt2,
         e,
         edt2,
         ehall,
         egradpe,
         egradpedt2,
         moments,
         momentsdt2,
         dperb,
         dmoments,
         dmomentsdt2,
         bgb,
         technical,
         fsgrid,
         sysBoundaries,
         RK_ORDER2_STEP1,
         true // communicateEGradPeOrMomentsDerivatives
      );
 
      propagateMagneticFieldSimple(perb, perbdt2, bgb, e, edt2, technical, fsgrid, sysBoundaries, dt, RK_ORDER2_STEP2);
      calculateDerivativesSimple(perb, moments, dperb, dmoments, technical, fsgrid, true /*doMoments*/);
      if (P::ohmGradPeTerm > 0) {
         calculateGradPeTermSimple(egradpe, egradpedt2, moments, momentsdt2, dmoments, dmomentsdt2, technical, fsgrid, sysBoundaries, RK_ORDER2_STEP2);
      }
      if (P::ohmHallTerm > 0) {
         calculateHallTermSimple(
            perb,
            perbdt2,
            ehall,
            moments,
            momentsdt2,
            dperb,
            dmoments,
            dmomentsdt2,
            bgb,
            technical,
            fsgrid,
            sysBoundaries,
            RK_ORDER2_STEP2,
            true // communicateMomentsDerivatives
         );
      }
      calculateUpwindedElectricFieldSimple(
         perb,
         perbdt2,
         e,
         edt2,
         ehall,
         egradpe,
         egradpedt2,
         moments,
         momentsdt2,
         dperb,
         dmoments,
         dmomentsdt2,
         bgb,
         technical,
         fsgrid,
         sysBoundaries,
         RK_ORDER2_STEP2,
         true // communicateEGradPeOrMomentsDerivatives
      );
#endif
   } else {
      Real subcycleDt = dt / convert<Real>(subcycles);
      Real subcycleT = P::t;
      creal targetT = P::t + dt;
      uint subcycleCount = 0;
      uint maxSubcycleCount = std::numeric_limits<uint>::max();
      int myRank = fsgrid.getRank();
 
      while (subcycleCount < maxSubcycleCount) {
         // In case of subcycling, we decided to go for a blunt Runge-Kutta subcycling even though e.g. moments are not
         // going along. Result of the Summer of Debugging 2016, the behaviour in wave dispersion was much improved with
         // this.
         propagateMagneticFieldSimple(perb, perbdt2, bgb, e, edt2, technical, fsgrid, sysBoundaries, subcycleDt, RK_ORDER2_STEP1);
 
         // We need to calculate derivatives of the moments at every substep, but the moments only
         // need to be communicated in the first one.
         calculateDerivativesSimple(perbdt2, momentsdt2, dperb, dmomentsdt2, technical, fsgrid, (subcycleCount == 0) /*doMoments*/);
         if (P::ohmGradPeTerm > 0 && subcycleCount == 0) {
            calculateGradPeTermSimple(egradpe, egradpedt2, moments, momentsdt2, dmoments, dmomentsdt2, technical, fsgrid, sysBoundaries, RK_ORDER2_STEP1);
         }
         if (P::ohmHallTerm > 0) {
            calculateHallTermSimple(
               perb,
               perbdt2,
               ehall,
               moments,
               momentsdt2,
               dperb,
               dmoments,
               dmomentsdt2,
               bgb,
               technical,
               fsgrid,
               sysBoundaries,
               RK_ORDER2_STEP1,
               subcycleCount == 0 // communicateMomentsDerivatives
            );
         }
         calculateUpwindedElectricFieldSimple(
            perb,
            perbdt2,
            e,
            edt2,
            ehall,
            egradpe,
            egradpedt2,
            moments,
            momentsdt2,
            dperb,
            dmoments,
            dmomentsdt2,
            bgb,
            technical,
            fsgrid,
            sysBoundaries,
            RK_ORDER2_STEP1,
            subcycleCount == 0 // communicateEGradPeOrMomentsDerivatives
         );
 
         propagateMagneticFieldSimple(perb, perbdt2, bgb, e, edt2, technical, fsgrid, sysBoundaries, subcycleDt, RK_ORDER2_STEP2);
 
         // We need to calculate derivatives of the moments at every substep, but the moments only
         // need to be communicated in the first one.
         calculateDerivativesSimple(perb, moments, dperb, dmoments, technical, fsgrid, (subcycleCount == 0) /*doMoments*/);
         if (P::ohmGradPeTerm > 0 && subcycleCount == 0) {
            calculateGradPeTermSimple(egradpe, egradpedt2, moments, momentsdt2, dmoments, dmomentsdt2, technical, fsgrid, sysBoundaries, RK_ORDER2_STEP2);
         }
         if (P::ohmHallTerm > 0) {
            calculateHallTermSimple(
               perb,
               perbdt2,
               ehall,
               moments,
               momentsdt2,
               dperb,
               dmoments,
               dmomentsdt2,
               bgb,
               technical,
               fsgrid,
               sysBoundaries,
               RK_ORDER2_STEP2,
               subcycleCount == 0 // communicateMomentsDerivatives
            );
         }
         calculateUpwindedElectricFieldSimple(
            perb,
            perbdt2,
            e,
            edt2,
            ehall,
            egradpe,
            egradpedt2,
            moments,
            momentsdt2,
            dperb,
            dmoments,
            dmomentsdt2,
            bgb,
            technical,
            fsgrid,
            sysBoundaries,
            RK_ORDER2_STEP2,
            subcycleCount == 0 // communicateEGradPeOrMomentsDerivatives
         );
 
         phiprof::Timer subcyclingTimer{"FS subcycle stuff"};
         subcycleT += subcycleDt;
         subcycleCount++;
 
         if (subcycleT >= targetT || subcycleCount >= maxSubcycleCount) {
            // we are done
            if (subcycleT > targetT) {
               // due to roundoff we might hit this, should add delta
               std::cerr << "subcycleT > targetT, should not happen! (values: subcycleT " << subcycleT << ", subcycleDt " << subcycleDt << ", targetT " << targetT << ")" << std::endl;
            }
            break;
         }
 
         // Reassess subcycle dt
         Real dtMaxGlobal = 0.0;
         Real dtMaxLocal = fsgrid.parallel_reduction([](int timerId) ->  phiprof::Timer { return phiprof::Timer{timerId}; },
                                                   phiprof::initializeTimer("compute-subcycle-dt-reduction-loop"), technical,
                                                   [](Real a, Real b) { return std::min<Real>(a, b); },
                                                   std::numeric_limits<Real>::max(),
                                                   [=](const fsgrid::Coordinates &coordinates, const fsgrid::FsStencil& stencil, cuint sysBoundaryFlag, cuint sysBoundaryLayer, creal maximum) {
            if (sysBoundaryFlag == sysboundarytype::NOT_SYSBOUNDARY ||
               (sysBoundaryLayer == 1 && sysBoundaryFlag != sysboundarytype::NOT_SYSBOUNDARY)) {
               return technical[stencil.ooo()].maxFsDt;
            } else {
               return maximum;
            }
         });
 
         phiprof::Timer allreduceTimer{"MPI_Allreduce"};
         fsgrid.Allreduce(&(dtMaxLocal), &(dtMaxGlobal), 1, MPI_Type<Real>(), MPI_MIN);
         allreduceTimer.stop();
 
         // reduce dt if it is too high
         if (subcycleDt > dtMaxGlobal * P::fieldSolverMaxCFL) {
            creal meanFieldsCFL = 0.5 * (P::fieldSolverMaxCFL + P::fieldSolverMinCFL);
            subcycleDt = meanFieldsCFL * dtMaxGlobal;
            if (myRank == MASTER_RANK) {
               logFile << "(TIMESTEP) New field solver subcycle dt = " << subcycleDt << " computed on step " << P::tstep << " and substep " << subcycleCount << " at " << P::t << " s" << std::endl;
            }
         }
 
         // Readjust the dt to hit targetT. Try to avoid having a very
         // short delta step at the end, instead 2 more normal ones
         if (subcycleT + 1.5 * subcycleDt > targetT) {
            subcycleDt = targetT - subcycleT;
            maxSubcycleCount = subcycleCount + 1; // 1 more steps
            // check that subcyclDt has correct CFL, take 2 if not
            if (subcycleDt > dtMaxGlobal * P::fieldSolverMaxCFL) {
               subcycleDt = (targetT - subcycleT) / 2;
               maxSubcycleCount = subcycleCount + 2;
            }
         }
 
         subcyclingTimer.stop();
      }
 
      if (subcycles != subcycleCount && myRank == MASTER_RANK) {
         logFile << "Effective field solver subcycles were " << subcycleCount << " instead of " << P::fieldSolverSubcycles << " on step " << P::tstep << std::endl;
      }
   }
 
   calculateVolumeAveragedFieldsSimple(perb, e, dperb, vol, technical, fsgrid);
   calculateBVOLDerivativesSimple(vol, technical, fsgrid);
   if (FieldTracing::fieldTracingParameters.doTraceFullBox || Parameters::computeCurvature) {
      fsgrid.updateGhostCells(vol);
      calculateCurvatureSimple(vol, bgb, technical, fsgrid);
   }
   return true;
}