MantleConvection.cpp 55.1 KB
Newer Older
Nils Kohl's avatar
Nils Kohl committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
/*
 * Copyright (c) 2017-2020 Nils Kohl.
 *
 * This file is part of HyTeG
 * (see https://i10git.cs.fau.de/hyteg/hyteg).
 *
 * 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 3 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, see <http://www.gnu.org/licenses/>.
 */

21
#include <algorithm>
Nils Kohl's avatar
Nils Kohl committed
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#include <cmath>
#include <core/Environment.h>
#include <core/math/Constants.h>

#include "core/DataTypes.h"
#include "core/config/Config.h"
#include "core/mpi/MPIManager.h"

#include "hyteg/MeshQuality.hpp"
#include "hyteg/composites/StrongFreeSlipWrapper.hpp"
#include "hyteg/composites/UnsteadyDiffusion.hpp"
#include "hyteg/dataexport/SQL.hpp"
#include "hyteg/dataexport/TimingOutput.hpp"
#include "hyteg/dataexport/VTKOutput.hpp"
Nils Kohl's avatar
Nils Kohl committed
36
#include "hyteg/elementwiseoperators/P2P1ElementwiseBlendingStokesOperator.hpp"
37
#include "hyteg/functions/FunctionProperties.hpp"
Nils Kohl's avatar
Nils Kohl committed
38
39
#include "hyteg/geometry/AnnulusMap.hpp"
#include "hyteg/geometry/IcosahedralShellMap.hpp"
40
#include "hyteg/geometry/SphereTools.hpp"
Nils Kohl's avatar
Nils Kohl committed
41
42
#include "hyteg/gridtransferoperators/P2P1StokesToP2P1StokesProlongation.hpp"
#include "hyteg/gridtransferoperators/P2P1StokesToP2P1StokesRestriction.hpp"
43
#include "hyteg/memory/MemoryAllocation.hpp"
Nils Kohl's avatar
Nils Kohl committed
44
45
#include "hyteg/mesh/MeshInfo.hpp"
#include "hyteg/numerictools/CFDHelpers.hpp"
46
#include "hyteg/numerictools/SphericalHarmonicsTool.hpp"
Nils Kohl's avatar
Nils Kohl committed
47
48
49
50
51
52
53
54
55
56
57
#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
#include "hyteg/p2functionspace/P2ConstantOperator.hpp"
#include "hyteg/p2functionspace/P2Function.hpp"
#include "hyteg/p2functionspace/P2ProjectNormalOperator.hpp"
#include "hyteg/petsc/PETScBlockPreconditionedStokesSolver.hpp"
#include "hyteg/petsc/PETScLUSolver.hpp"
#include "hyteg/petsc/PETScManager.hpp"
#include "hyteg/petsc/PETScMinResSolver.hpp"
#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
#include "hyteg/primitivestorage/Visualization.hpp"
Nils Kohl's avatar
Nils Kohl committed
58
#include "hyteg/primitivestorage/loadbalancing/SimpleBalancer.hpp"
Nils Kohl's avatar
Nils Kohl committed
59
#include "hyteg/solvers/CGSolver.hpp"
60
#include "hyteg/solvers/FullMultigridSolver.hpp"
Nils Kohl's avatar
Nils Kohl committed
61
62
63
64
65
66
#include "hyteg/solvers/GaussSeidelSmoother.hpp"
#include "hyteg/solvers/GeometricMultigridSolver.hpp"
#include "hyteg/solvers/MinresSolver.hpp"
#include "hyteg/solvers/UzawaSmoother.hpp"
#include "hyteg/solvers/WeightedJacobiSmoother.hpp"
#include "hyteg/solvers/controlflow/SolverLoop.hpp"
67
#include "hyteg/solvers/preconditioners/stokes/StokesBlockDiagonalPreconditioner.hpp"
Nils Kohl's avatar
Nils Kohl committed
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
#include "hyteg/solvers/preconditioners/stokes/StokesPressureBlockPreconditioner.hpp"
#include "hyteg/solvers/preconditioners/stokes/StokesVelocityBlockBlockDiagonalPreconditioner.hpp"
#include "hyteg/solvers/solvertemplates/StokesSolverTemplates.hpp"

#include "coupling_hyteg_convection_particles/MMOCTransport.hpp"
#include "sqlite/SQLite.h"

namespace hyteg {

using walberla::int_c;
using walberla::real_t;
using walberla::math::pi;

#if 0
static std::string getDateTimeID()
{
   std::vector< char > cTimeString( 64 );
   WALBERLA_ROOT_SECTION()
   {
      std::time_t t;
      std::time( &t );
      std::strftime( cTimeString.data(), 64, "%F_%H-%M-%S", std::localtime( &t ) );
   }

   walberla::mpi::broadcastObject( cTimeString );

   std::string timeString( cTimeString.data() );
   return timeString;
}
#endif

Nils Kohl's avatar
Nils Kohl committed
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
enum class StokesSolverType : int
{
   PETSC_MUMPS         = 0,
   PETSC_MINRES_JACOBI = 1,
   PETSC_MINRES_BOOMER = 2,
   HYTEG_MINRES        = 3,
   HYTEG_MINRES_GMG    = 4,
   HYTEG_UZAWA_V       = 5,
   HYTEG_UZAWA_FMG     = 6
};

enum class DiffusionSolverType : int
{
   PETSC_MINRES = 0,
   HYTEG_CG     = 1,
   HYTEG_GMG    = 2
};

Nils Kohl's avatar
Nils Kohl committed
117
118
119
120
121
122
123
124
125
struct DomainInfo
{
   bool threeDim = false;

   real_t rMin = 0;
   real_t rMax = 0;
   uint_t nTan = 0;
   uint_t nRad = 0;

Nils Kohl's avatar
Nils Kohl committed
126
   real_t domainVolume() const
Nils Kohl's avatar
Nils Kohl committed
127
128
129
130
131
132
133
   {
      if ( !threeDim )
      {
         return walberla::math::pi * rMax * rMax - walberla::math::pi * rMin * rMin;
      }
      else
      {
Nils Kohl's avatar
Nils Kohl committed
134
         return ( 4. / 3. ) * walberla::math::pi * rMax * rMax * rMax - ( 4. / 3. ) * walberla::math::pi * rMin * rMin * rMin;
Nils Kohl's avatar
Nils Kohl committed
135
136
137
138
      }
   }
};

Nils Kohl's avatar
Nils Kohl committed
139
140
141
142
143
144
struct SolverInfo
{
   StokesSolverType stokesSolverType = StokesSolverType::PETSC_MUMPS;

   uint_t stokesMaxNumIterations           = 10;
   real_t stokesAbsoluteResidualUTolerance = 0;
145
   real_t stokesRelativeResidualUTolerance = 0;
146
147
148
   uint_t uzawaInnerIterations             = 10;
   uint_t uzawaPreSmooth                   = 6;
   uint_t uzawaPostSmooth                  = 6;
Nils Kohl's avatar
Nils Kohl committed
149
150
151
152
153

   DiffusionSolverType diffusionSolverType = DiffusionSolverType::PETSC_MINRES;

   uint_t diffusionMaxNumIterations           = 10000;
   real_t diffusionAbsoluteResidualUTolerance = 10000;
Nils Kohl's avatar
Nils Kohl committed
154
155

   real_t uzawaOmega = 0.3;
Nils Kohl's avatar
Nils Kohl committed
156
157
};

158
159
160
161
162
163
164
struct OutputInfo
{
   bool        vtk;
   bool        vtkOutputVelocity;
   bool        sphericalTemperatureSlice;
   int         sphericalTemperatureSliceNumMeridians;
   int         sphericalTemperatureSliceNumParallels;
165
   int         sphericalTemperatureSliceIcoRefinements;
166
167
168
169
170
171
172
173
   uint_t      printInterval;
   uint_t      vtkInterval;
   uint_t      sphericalTemperatureSliceInterval;
   bool        vtkOutputVertexDoFs;
   std::string outputDirectory;
   std::string outputBaseName;
};

Nils Kohl's avatar
Nils Kohl committed
174
175
176
177
178
179
180
181
182
183
184
185
/// Calculates and returns
///
///     ||u||_L2 = sqrt( u^T M u )
///
template < typename FunctionType, typename MassOperator >
real_t normL2( const FunctionType& u, const FunctionType& tmp, const MassOperator& M, const uint_t& level, const DoFType& flag )
{
   tmp.interpolate( 0, level );
   M.apply( u, tmp, level, flag );
   return std::sqrt( u.dotGlobal( tmp, level, flag ) );
}

186
187
188
189
190
191
192
193
194
195
196
197
template < typename FunctionType, typename MassOperator >
real_t normL2Squared( const FunctionType& u,
                      const FunctionType& tmp,
                      const MassOperator& M,
                      const uint_t&       level,
                      const DoFType&      flag )
{
   tmp.interpolate( 0, level );
   M.apply( u, tmp, level, flag );
   return u.dotGlobal( tmp, level, flag );
}

Nils Kohl's avatar
Nils Kohl committed
198
199
200
201
202
203
204
205
206
template < typename StokesFunction, typename StokesOperator, typename MassOperatorVelocity, typename MassOperatorPressure >
void calculateStokesResiduals( const StokesOperator&       A,
                               const MassOperatorVelocity& Mu,
                               const MassOperatorPressure& Mp,
                               const StokesFunction&       x,
                               const StokesFunction&       b,
                               uint_t                      level,
                               StokesFunction&             r,
                               StokesFunction&             tmp,
207
                               real_t&                     residual )
Nils Kohl's avatar
Nils Kohl committed
208
209
210
211
{
   tmp.interpolate( 0, level, All );
   r.interpolate( 0, level, All );
   A.apply( x, tmp, level, Inner | NeumannBoundary | FreeslipBoundary );
212
213
214
   r.assign( { 1.0, -1.0 }, { b, tmp }, level, Inner | NeumannBoundary | FreeslipBoundary );
   residual = normL2Squared( r.uvw[0], tmp.uvw[0], Mu, level, Inner | NeumannBoundary | FreeslipBoundary );
   residual += normL2Squared( r.uvw[1], tmp.uvw[1], Mu, level, Inner | NeumannBoundary | FreeslipBoundary );
Nils Kohl's avatar
Nils Kohl committed
215
216
   if ( x.getStorage()->hasGlobalCells() )
   {
217
      residual += normL2Squared( r.uvw[2], tmp.uvw[2], Mu, level, Inner | NeumannBoundary | FreeslipBoundary );
Nils Kohl's avatar
Nils Kohl committed
218
   }
219
220
   residual += normL2Squared( r.p, tmp.p, Mp, level, Inner | NeumannBoundary | FreeslipBoundary );
   residual = std::sqrt( residual );
Nils Kohl's avatar
Nils Kohl committed
221
222
}

Nils Kohl's avatar
Nils Kohl committed
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
template < typename UnsteadyDiffusion,
           typename UnsteadyDiffusionOperator,
           typename LaplaceOperator,
           typename MassOperatorVelocity,
           typename ScalarFuncionType >
void calculateDiffusionResidual( UnsteadyDiffusion&               unsteadyDiffusion,
                                 const UnsteadyDiffusionOperator& unsteadyDiffusionOperator,
                                 const LaplaceOperator&           laplacian,
                                 const MassOperatorVelocity&      Mu,
                                 const ScalarFuncionType&         c,
                                 const ScalarFuncionType&         cOld,
                                 const ScalarFuncionType&         f,
                                 ScalarFuncionType&               r,
                                 ScalarFuncionType&               tmp,
                                 uint_t                           level,
                                 real_t&                          residual )
{
   unsteadyDiffusion.calculateResidual(
       unsteadyDiffusionOperator, laplacian, Mu, c, cOld, f, f, r, level, Inner | NeumannBoundary | FreeslipBoundary );
   residual = normL2( r, tmp, Mu, level, Inner | NeumannBoundary | FreeslipBoundary );
}

Nils Kohl's avatar
Nils Kohl committed
245
246
247
template < typename StokesFunction, typename VelocityMass >
real_t velocityRMS( const StokesFunction& u, const StokesFunction& tmp, const VelocityMass& M, real_t domainVolume, uint_t level )
{
248
249
   auto norm = std::pow( normL2( u.uvw[0], tmp.uvw[0], M, level, All ), 2.0 ) +
               std::pow( normL2( u.uvw[1], tmp.uvw[1], M, level, All ), 2.0 );
Nils Kohl's avatar
Nils Kohl committed
250
251
252

   if ( u.getStorage()->hasGlobalCells() )
   {
253
      norm += std::pow( normL2( u.uvw[2], tmp.uvw[2], M, level, All ), 2.0 );
Nils Kohl's avatar
Nils Kohl committed
254
255
256
257
258
   }

   return std::sqrt( norm / domainVolume );
}

259
std::shared_ptr< SetupPrimitiveStorage > createSetupStorage( DomainInfo domainInfo )
Nils Kohl's avatar
Nils Kohl committed
260
261
262
263
264
265
266
267
268
269
270
271
272
273
{
   MeshInfo meshInfo = MeshInfo::emptyMeshInfo();
   if ( domainInfo.threeDim )
   {
      meshInfo = MeshInfo::meshSphericalShell( domainInfo.nTan, domainInfo.nRad, domainInfo.rMin, domainInfo.rMax );
   }
   else
   {
      meshInfo = MeshInfo::meshAnnulus( domainInfo.rMin, domainInfo.rMax, MeshInfo::CRISS, domainInfo.nTan, domainInfo.nRad );
   }

   auto setupStorage = std::make_shared< SetupPrimitiveStorage >(
       meshInfo, walberla::uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );

Nils Kohl's avatar
Nils Kohl committed
274
275
   loadbalancing::roundRobinVolume( *setupStorage );

276
   setupStorage->setMeshBoundaryFlagsOnBoundary( 1, 0, true );
Nils Kohl's avatar
Nils Kohl committed
277
278
279
280
281
282
283
284
285
286
287
288
289

   if ( domainInfo.threeDim )
   {
      IcosahedralShellMap::setMap( *setupStorage );
   }
   else
   {
      AnnulusMap::setMap( *setupStorage );
   }

   return setupStorage;
}

290
291
292
293
294
295
296
297
298
299
300
301
302
void runBenchmark( real_t     cflMax,
                   real_t     rayleighNumber,
                   bool       fixedTimeStep,
                   real_t     dtConstant,
                   uint_t     maxNumTimeSteps,
                   uint_t     minLevel,
                   uint_t     level,
                   SolverInfo solverInfo,
                   DomainInfo domainInfo,
                   bool       predictorCorrector,
                   real_t     simulationTime,
                   OutputInfo outputInfo,
                   bool       verbose )
Nils Kohl's avatar
Nils Kohl committed
303
304
305
{
   walberla::WcTimer localTimer;

306
   const bool outputTimingJSON = true;
Nils Kohl's avatar
Nils Kohl committed
307

Nils Kohl's avatar
Nils Kohl committed
308
309
310
311
312
   if ( verbose )
   {
      WALBERLA_LOG_INFO_ON_ROOT( "Building setup storage ..." );
   }

313
   auto setupStorage = createSetupStorage( domainInfo );
Nils Kohl's avatar
Nils Kohl committed
314
315
316
317
318
319

   if ( verbose )
   {
      WALBERLA_LOG_INFO_ON_ROOT( "Building distributed storage ..." );
   }

320
   auto storage = std::make_shared< PrimitiveStorage >( *setupStorage, 1 );
Nils Kohl's avatar
Nils Kohl committed
321

322
   if ( outputInfo.vtk )
Nils Kohl's avatar
Nils Kohl committed
323
   {
Nils Kohl's avatar
Nils Kohl committed
324
325
326
327
328
      if ( verbose )
      {
         WALBERLA_LOG_INFO_ON_ROOT( "Writing domain partitioning VTK ..." );
      }

329
      writeDomainPartitioningVTK( storage, outputInfo.outputDirectory, outputInfo.outputBaseName + "_domain" );
Nils Kohl's avatar
Nils Kohl committed
330
331
332
   }

   auto timer = storage->getTimingTree();
333

Nils Kohl's avatar
Nils Kohl committed
334
335
   timer->start( "Setup" );

Nils Kohl's avatar
Nils Kohl committed
336
337
338
339
340
341
342
343
344
   if ( verbose )
   {
      WALBERLA_LOG_INFO_ON_ROOT( "Gathering domain information ..." );
   }

   const uint_t unknownsTemperature = numberOfGlobalDoFs< P2FunctionTag >( *storage, level );
   const uint_t unknownsStokes      = numberOfGlobalDoFs< P2P1TaylorHoodFunctionTag >( *storage, level );
   const real_t hMin                = MeshQuality::getMinimalEdgeLength( storage, level );
   const real_t hMax                = MeshQuality::getMaximalEdgeLength( storage, level );
Nils Kohl's avatar
Nils Kohl committed
345

346
347
   const real_t diffusivity     = 1.0;
   const real_t internalHeating = 0.0;
348
349

   const real_t sliceEvaluationRadius = domainInfo.rMin + 0.5 * ( domainInfo.rMax - domainInfo.rMin );
Nils Kohl's avatar
Nils Kohl committed
350

Nils Kohl's avatar
Nils Kohl committed
351
   WALBERLA_LOG_INFO_ON_ROOT( "" )
352
   WALBERLA_LOG_INFO_ON_ROOT( "Benchmark name: " << outputInfo.outputBaseName )
Nils Kohl's avatar
Nils Kohl committed
353
354
355
356
357
358
359
360
361
362
363
364
   WALBERLA_LOG_INFO_ON_ROOT( " - time discretization: " )
   WALBERLA_LOG_INFO_ON_ROOT( "   + fixed time step:                              " << fixedTimeStep )
   if ( fixedTimeStep )
   {
      WALBERLA_LOG_INFO_ON_ROOT( "   + dt (constant):                                " << dtConstant )
   }
   else
   {
      WALBERLA_LOG_INFO_ON_ROOT( "   + max CFL:                                      " << cflMax )
   }
   WALBERLA_LOG_INFO_ON_ROOT( "   + simulation time:                              " << simulationTime )
   WALBERLA_LOG_INFO_ON_ROOT( " - space discretization: " )
Nils Kohl's avatar
Nils Kohl committed
365
366
367
368
   WALBERLA_LOG_INFO_ON_ROOT( "   + rMin:                                         " << domainInfo.rMin )
   WALBERLA_LOG_INFO_ON_ROOT( "   + rMax:                                         " << domainInfo.rMax )
   WALBERLA_LOG_INFO_ON_ROOT( "   + nTan:                                         " << domainInfo.nTan )
   WALBERLA_LOG_INFO_ON_ROOT( "   + nRad:                                         " << domainInfo.nRad )
Nils Kohl's avatar
Nils Kohl committed
369
370
   WALBERLA_LOG_INFO_ON_ROOT( "   + dimensions:                                   " << ( storage->hasGlobalCells() ? "3" : "2" ) )
   WALBERLA_LOG_INFO_ON_ROOT( "   + level:                                        " << level )
Nils Kohl's avatar
Nils Kohl committed
371
   WALBERLA_LOG_INFO_ON_ROOT( "   + unknowns temperature, including boundary:     " << unknownsTemperature )
372
373
374
375
376
   for ( uint_t l = minLevel; l <= level; l++ )
   {
      const uint_t unknownsStokesLevel = numberOfGlobalDoFs< P2P1TaylorHoodFunctionTag >( *storage, l );
      WALBERLA_LOG_INFO_ON_ROOT( "   + unknowns Stokes, including boundary, level " << l << ": " << unknownsStokesLevel )
   }
Nils Kohl's avatar
Nils Kohl committed
377
378
379
380
381
382
   WALBERLA_LOG_INFO_ON_ROOT( "   + h_min:                                        " << hMin )
   WALBERLA_LOG_INFO_ON_ROOT( "   + h_max:                                        " << hMax )
   WALBERLA_LOG_INFO_ON_ROOT( " - benchmark settings: " )
   WALBERLA_LOG_INFO_ON_ROOT( "   + Rayleigh number:                              " << rayleighNumber )
   WALBERLA_LOG_INFO_ON_ROOT( "   + internal heating:                             " << internalHeating )
   WALBERLA_LOG_INFO_ON_ROOT( " - app settings: " )
Nils Kohl's avatar
Nils Kohl committed
383
   WALBERLA_LOG_INFO_ON_ROOT( "   + Stokes solver:                                " << int_c( solverInfo.stokesSolverType ) )
Nils Kohl's avatar
Nils Kohl committed
384
   WALBERLA_LOG_INFO_ON_ROOT( "   + Uzawa omega:                                  " << solverInfo.uzawaOmega )
Nils Kohl's avatar
Nils Kohl committed
385
   WALBERLA_LOG_INFO_ON_ROOT( "   + diffusion solver:                             " << int_c( solverInfo.diffusionSolverType ) )
Nils Kohl's avatar
Nils Kohl committed
386

387
388
389
390
391
   if ( outputInfo.vtk )
   {
      WALBERLA_LOG_INFO_ON_ROOT( "   + VTK interval:                                 " << outputInfo.vtkInterval )
   }
   if ( outputInfo.sphericalTemperatureSlice )
Nils Kohl's avatar
Nils Kohl committed
392
   {
393
394
      WALBERLA_LOG_INFO_ON_ROOT(
          "   + spherical slice output interval:              " << outputInfo.sphericalTemperatureSliceInterval )
Nils Kohl's avatar
Nils Kohl committed
395
   }
396
397
398
   WALBERLA_LOG_INFO_ON_ROOT( "   + print interval:                               " << outputInfo.printInterval )
   WALBERLA_LOG_INFO_ON_ROOT( "   + output directory:                             " << outputInfo.outputDirectory )
   WALBERLA_LOG_INFO_ON_ROOT( "   + output base name:                             " << outputInfo.outputBaseName )
Nils Kohl's avatar
Nils Kohl committed
399
400
   WALBERLA_LOG_INFO_ON_ROOT( "" )

Nils Kohl's avatar
Nils Kohl committed
401
402
403
404
   auto storageInfo = storage->getGlobalInfo();
   WALBERLA_LOG_INFO_ON_ROOT( storageInfo );
   WALBERLA_LOG_INFO_ON_ROOT( "" );

405
   FixedSizeSQLDB db( outputInfo.outputDirectory + "/" + outputInfo.outputBaseName + ".db" );
Nils Kohl's avatar
Nils Kohl committed
406

Nils Kohl's avatar
Nils Kohl committed
407
   db.setConstantEntry( "ra", rayleighNumber );
408
409
410
411
   db.setConstantEntry( "ntan", domainInfo.nTan );
   db.setConstantEntry( "nrad", domainInfo.nRad );
   db.setConstantEntry( "rmin", domainInfo.rMin );
   db.setConstantEntry( "rmax", domainInfo.rMax );
Nils Kohl's avatar
Nils Kohl committed
412
413
414
415
416
   db.setConstantEntry( "fixed_time_step", fixedTimeStep );
   db.setConstantEntry( "cfl_max", cflMax );
   db.setConstantEntry( "dt_constant", dtConstant );
   db.setConstantEntry( "simulation_time", simulationTime );
   db.setConstantEntry( "level", uint_c( level ) );
Nils Kohl's avatar
Nils Kohl committed
417
418
   db.setConstantEntry( "unknowns_temperature", uint_c( unknownsTemperature ) );
   db.setConstantEntry( "unknowns_stokes", uint_c( unknownsStokes ) );
Nils Kohl's avatar
Nils Kohl committed
419
420
421
422
423
424
425
426
427
   db.setConstantEntry( "h_min", hMin );
   db.setConstantEntry( "h_max", hMax );
   db.setConstantEntry( "num_macro_cells", uint_c( setupStorage->getNumberOfCells() ) );
   db.setConstantEntry( "num_macro_faces", uint_c( setupStorage->getNumberOfFaces() ) );
   db.setConstantEntry( "num_macro_edges", uint_c( setupStorage->getNumberOfEdges() ) );
   db.setConstantEntry( "num_macro_vertices", uint_c( setupStorage->getNumberOfVertices() ) );
   db.setConstantEntry( "num_macro_primitives", uint_c( setupStorage->getNumberOfPrimitives() ) );
   db.setConstantEntry( "diffusivity", diffusivity );

428
   db.setConstantEntry( "stokes_solver_type", int64_c( solverInfo.stokesSolverType ) );
429
430
431
432
433
   db.setConstantEntry( "uzawa_omega", solverInfo.uzawaOmega );
   db.setConstantEntry( "uzawa_inner_smooth", solverInfo.uzawaInnerIterations );
   db.setConstantEntry( "uzawa_pre_smooth", solverInfo.uzawaPreSmooth );
   db.setConstantEntry( "uzawa_post_smooth", solverInfo.uzawaPostSmooth );

434
435
436
   db.setConstantEntry( "stokes_absolute_residual_threshold", solverInfo.stokesAbsoluteResidualUTolerance );
   db.setConstantEntry( "stokes_relative_residual_threshold", solverInfo.stokesRelativeResidualUTolerance );

437
438
439
440
441
   typedef P2P1TaylorHoodFunction< real_t >       StokesFunction;
   typedef P2Function< real_t >                   ScalarFunction;
   typedef P2P1ElementwiseBlendingStokesOperator  StokesOperator;
   typedef P2ElementwiseBlendingLaplaceOperator   LaplaceOperator;
   typedef P2ElementwiseBlendingMassOperator      MassOperatorVelocity;
Nils Kohl's avatar
Nils Kohl committed
442
   typedef P1ConstantMassOperator                 MassOperatorPressure;
443
   typedef P2ElementwiseUnsteadyDiffusionOperator UnsteadyDiffusionOperator;
Nils Kohl's avatar
Nils Kohl committed
444
445
446
447

   BoundaryCondition bcVelocity;
   BoundaryCondition bcTemperature;

448
449
450
451
   bcVelocity.createDirichletBC( "outerBoundaryVelocity", 1 );
   bcVelocity.createDirichletBC( "innerBoundaryVelocity", 1 );
   bcTemperature.createDirichletBC( "outerBoundaryTemperature", 1 );
   bcTemperature.createDirichletBC( "innerBoundaryTemperature", 1 );
Nils Kohl's avatar
Nils Kohl committed
452

Nils Kohl's avatar
Nils Kohl committed
453
454
455
456
457
   if ( verbose )
   {
      WALBERLA_LOG_INFO_ON_ROOT( "Allocating functions ..." );
   }

458
459
460
461
462
463
   ScalarFunction c( "c", storage, minLevel, level, bcTemperature );
   ScalarFunction cPr( "cPr", storage, minLevel, level, bcTemperature );
   ScalarFunction cOld( "cOld", storage, minLevel, level, bcTemperature );
   ScalarFunction cTmp( "cTmp", storage, minLevel, level, bcTemperature );
   ScalarFunction cTmp2( "cTmp2", storage, minLevel, level, bcTemperature );
   ScalarFunction q( "q", storage, minLevel, level, bcTemperature );
Nils Kohl's avatar
Nils Kohl committed
464

465
466
467
468
469
470
   StokesFunction u( "u", storage, minLevel, level, bcVelocity );
   StokesFunction uLast( "uLast", storage, minLevel, level, bcVelocity );
   StokesFunction f( "f", storage, minLevel, level, bcVelocity );
   StokesFunction outwardNormal( "outwardNormal", storage, minLevel, level, bcVelocity );
   StokesFunction stokesTmp( "stokesTmp", storage, minLevel, level, bcVelocity );
   StokesFunction stokesResidual( "stokesResidual", storage, minLevel, level, bcVelocity );
Nils Kohl's avatar
Nils Kohl committed
471

Nils Kohl's avatar
Nils Kohl committed
472
473
   ScalarFunction uMagnitudeSquared( "uMagnitudeSquared", storage, minLevel, level, bcVelocity );

474
475
   ScalarFunction uTmp( "uTmp", storage, minLevel, level, bcVelocity );
   ScalarFunction uTmp2( "uTmp2", storage, minLevel, level, bcVelocity );
Nils Kohl's avatar
Nils Kohl committed
476

477
478
   SphericalHarmonicsTool sphTool( 6 );

Nils Kohl's avatar
Nils Kohl committed
479
480
   std::function< real_t( const Point3D& ) > initialTemperature = [&]( const Point3D& x ) {
      auto radius = std::sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] );
481
482
483
484
485
486
487
488
489
490
491

      const auto linearSlope =
          std::max( real_c( 0 ), 1 - ( ( radius - domainInfo.rMin ) / ( domainInfo.rMax - domainInfo.rMin ) ) );

      const auto xn = x / x.norm();

      // sph ~ in (-2, 2)
      const auto powerPertubation = sphTool.shconvert_eval( 4, 3, xn[0], xn[1], xn[2] );

      const auto powerSlope = std::pow( linearSlope, 4 + powerPertubation );
      return powerSlope;
Nils Kohl's avatar
Nils Kohl committed
492
493
   };

Nils Kohl's avatar
Nils Kohl committed
494
495
496
497
498
   if ( verbose )
   {
      WALBERLA_LOG_INFO_ON_ROOT( "Interpolating initial temperature and internal heating ..." );
   }

499
   for ( uint_t l = minLevel; l <= level; l++ )
500
501
502
503
   {
      c.interpolate( initialTemperature, l, All );
      q.interpolate( internalHeating, l, All );
   }
Nils Kohl's avatar
Nils Kohl committed
504
505
506
507
508

   std::function< real_t( const Point3D& ) > normalX = []( const Point3D& x ) { return x[0] / x.norm(); };
   std::function< real_t( const Point3D& ) > normalY = []( const Point3D& x ) { return x[1] / x.norm(); };
   std::function< real_t( const Point3D& ) > normalZ = []( const Point3D& x ) { return x[2] / x.norm(); };

Nils Kohl's avatar
Nils Kohl committed
509
510
511
512
   if ( verbose )
   {
      WALBERLA_LOG_INFO_ON_ROOT( "Interpolating normals ..." );
   }
Nils Kohl's avatar
Nils Kohl committed
513

514
   for ( uint_t l = minLevel; l <= level; l++ )
Nils Kohl's avatar
Nils Kohl committed
515
   {
516
      outwardNormal.uvw.interpolate( { normalX, normalY, normalZ }, l );
Nils Kohl's avatar
Nils Kohl committed
517
518
   }

Nils Kohl's avatar
Nils Kohl committed
519
520
521
522
523
   if ( verbose )
   {
      WALBERLA_LOG_INFO_ON_ROOT( "Preparing operators ..." );
   }

524
525
526
527
528
529
   UnsteadyDiffusionOperator diffusionOperator(
       storage, minLevel, level, 1.0, diffusivity, DiffusionTimeIntegrator::ImplicitEuler );
   auto                            A = std::make_shared< StokesOperator >( storage, minLevel, level );
   LaplaceOperator                 L( storage, minLevel, level );
   MassOperatorVelocity            MVelocity( storage, minLevel, level );
   MassOperatorPressure            MPressure( storage, minLevel, level );
530
   MMOCTransport< ScalarFunction > transport( storage, minLevel, level, TimeSteppingScheme::RK4 );
Nils Kohl's avatar
Nils Kohl committed
531

532
533
534
535
536
537
538
   if ( storage->hasGlobalCells() )
   {
      A->computeAndStoreLocalElementMatrices();
      L.computeAndStoreLocalElementMatrices();
      MVelocity.computeAndStoreLocalElementMatrices();
   }

Nils Kohl's avatar
Nils Kohl committed
539
540
541
542
543
   if ( verbose )
   {
      WALBERLA_LOG_INFO_ON_ROOT( "Preparing solvers ..." );
   }

544
545
   std::shared_ptr< Solver< StokesOperator > > stokesSolver;
   std::shared_ptr< Solver< StokesOperator > > stokesSolverBlockPrecMinRes;
Nils Kohl's avatar
Nils Kohl committed
546

Nils Kohl's avatar
Nils Kohl committed
547
548
549
550
   real_t initialResiudalU          = 0;
   real_t vCycleResidualULast       = 0;
   real_t vCycleResidualCLast       = 0;
   uint_t numVCycles                = 0;
Nils Kohl's avatar
Nils Kohl committed
551
   real_t averageResidualReductionU = 0;
Nils Kohl's avatar
Nils Kohl committed
552

Nils Kohl's avatar
Nils Kohl committed
553
   if ( solverInfo.stokesSolverType == StokesSolverType::PETSC_MUMPS )
554
555
556
557
   {
      stokesSolver = std::make_shared< PETScLUSolver< StokesOperator > >( storage, level );
   }

Nils Kohl's avatar
Nils Kohl committed
558
559
   else if ( solverInfo.stokesSolverType == StokesSolverType::PETSC_MINRES_JACOBI ||
             solverInfo.stokesSolverType == StokesSolverType::PETSC_MINRES_BOOMER )
Nils Kohl's avatar
Nils Kohl committed
560
   {
561
      auto velocityPreconditionerType = 1;
Nils Kohl's avatar
Nils Kohl committed
562
      if ( solverInfo.stokesSolverType == StokesSolverType::PETSC_MINRES_BOOMER )
563
564
565
      {
         velocityPreconditionerType = 3;
      }
Nils Kohl's avatar
Nils Kohl committed
566
567
568
569
570
571
572
      auto stokesSolverTmp =
          std::make_shared< PETScBlockPreconditionedStokesSolver< StokesOperator > >( storage,
                                                                                      level,
                                                                                      solverInfo.stokesAbsoluteResidualUTolerance,
                                                                                      solverInfo.stokesMaxNumIterations,
                                                                                      velocityPreconditionerType,
                                                                                      1 );
Nils Kohl's avatar
Nils Kohl committed
573
      stokesSolverTmp->reassembleMatrix( false );
574
      stokesSolverTmp->setVerbose( verbose );
Nils Kohl's avatar
Nils Kohl committed
575
576
577
      stokesSolver = stokesSolverTmp;
   }

Nils Kohl's avatar
Nils Kohl committed
578
   else if ( solverInfo.stokesSolverType == StokesSolverType::HYTEG_MINRES )
579
580
   {
      stokesSolver = solvertemplates::stokesMinResSolver< StokesOperator >(
Nils Kohl's avatar
Nils Kohl committed
581
          storage, level, solverInfo.stokesAbsoluteResidualUTolerance, solverInfo.stokesMaxNumIterations, verbose );
582
   }
Nils Kohl's avatar
Nils Kohl committed
583

Nils Kohl's avatar
Nils Kohl committed
584
   else if ( solverInfo.stokesSolverType == StokesSolverType::HYTEG_MINRES_GMG )
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
   {
      typedef CGSolver< P2ElementwiseBlendingLaplaceOperator >                 CoarseGridSolver_T;
      typedef GeometricMultigridSolver< P2ElementwiseBlendingLaplaceOperator > GMGSolver_T;

      auto coarseGridSolver = std::make_shared< CoarseGridSolver_T >( storage, minLevel, level );
      auto smoother =
          std::make_shared< WeightedJacobiSmoother< P2ElementwiseBlendingLaplaceOperator > >( storage, minLevel, level, 0.66 );
      auto prolongationOperator = std::make_shared< P2toP2QuadraticProlongation >();
      auto restrictionOperator  = std::make_shared< P2toP2QuadraticRestriction >();
      auto gmgSolver            = std::make_shared< GMGSolver_T >(
          storage, smoother, coarseGridSolver, restrictionOperator, prolongationOperator, minLevel, level, 3, 3 );

      auto lumpedInvPressureMass = std::make_shared< P1LumpedInvMassOperator >( storage, minLevel, level );

      typedef StokesBlockDiagonalPreconditioner< StokesOperator, P1LumpedInvMassOperator > Preconditioner_T;
      auto preconditioner = std::make_shared< Preconditioner_T >( storage, minLevel, level, 3, gmgSolver );

Nils Kohl's avatar
Nils Kohl committed
602
603
      auto stokesSolverTmp = std::make_shared< MinResSolver< StokesOperator > >(
          storage, minLevel, level, solverInfo.stokesMaxNumIterations, 1e-30, preconditioner );
604
      stokesSolverTmp->setPrintInfo( verbose );
Nils Kohl's avatar
Nils Kohl committed
605
606
      stokesSolverTmp->setAbsoluteTolerance( solverInfo.stokesAbsoluteResidualUTolerance );
      if ( solverInfo.stokesSolverType == StokesSolverType::HYTEG_MINRES_GMG )
607
608
609
610
611
612
613
      {
         stokesSolver = stokesSolverTmp;
      }
      else
      {
         stokesSolverBlockPrecMinRes = stokesSolverTmp;
      }
Nils Kohl's avatar
Nils Kohl committed
614
615
   }

Nils Kohl's avatar
Nils Kohl committed
616
617
   else if ( solverInfo.stokesSolverType == StokesSolverType::HYTEG_UZAWA_V ||
             solverInfo.stokesSolverType == StokesSolverType::HYTEG_UZAWA_FMG )
618
619
620
621
622
623
624
625
626
627
628
   {
      auto prolongationOperator = std::make_shared< P2P1StokesToP2P1StokesProlongation >();
      auto restrictionOperator  = std::make_shared< P2P1StokesToP2P1StokesRestriction >( true );

      std::shared_ptr< Solver< typename StokesOperator::VelocityOperator_T > > smoother =
          std::make_shared< WeightedJacobiSmoother< typename StokesOperator::VelocityOperator_T > >(
              storage, minLevel, level, 0.66 );

      auto uzawaVelocityPreconditioner =
          std::make_shared< StokesVelocityBlockBlockDiagonalPreconditioner< StokesOperator > >( storage, smoother );

Nils Kohl's avatar
Nils Kohl committed
629
630
631
632
633
634
635
      auto uzawaSmoother = std::make_shared< UzawaSmoother< StokesOperator > >( storage,
                                                                                uzawaVelocityPreconditioner,
                                                                                minLevel,
                                                                                level,
                                                                                solverInfo.uzawaOmega,
                                                                                Inner | NeumannBoundary,
                                                                                solverInfo.uzawaInnerIterations );
636
637
638
639

      std::shared_ptr< Solver< StokesOperator > > coarseGridSolverInternal;

      auto petscSolverInternalTmp = std::make_shared< PETScBlockPreconditionedStokesSolver< StokesOperator > >(
640
          storage, minLevel, solverInfo.stokesAbsoluteResidualUTolerance, 1000, 1 );
641
642
643
644
645
646
647
648
649
650
      petscSolverInternalTmp->setVerbose( verbose );
      auto coarseGridSolver = petscSolverInternalTmp;

      auto multigridSolver = std::make_shared< GeometricMultigridSolver< StokesOperator > >( storage,
                                                                                             uzawaSmoother,
                                                                                             coarseGridSolver,
                                                                                             restrictionOperator,
                                                                                             prolongationOperator,
                                                                                             minLevel,
                                                                                             level,
651
652
                                                                                             solverInfo.uzawaPreSmooth,
                                                                                             solverInfo.uzawaPostSmooth,
653
654
655
                                                                                             2,
                                                                                             CycleType::VCYCLE );

Nils Kohl's avatar
Nils Kohl committed
656
      if ( solverInfo.stokesSolverType == StokesSolverType::HYTEG_UZAWA_V )
657
      {
Nils Kohl's avatar
Nils Kohl committed
658
659
660
         auto stopIterationCallback =
             [&]( const StokesOperator& _A, const StokesFunction& _u, const StokesFunction& _b, const uint_t _level ) {
                real_t r_u;
661

662
663
664
665
                calculateStokesResiduals( _A, MVelocity, MPressure, _u, _b, _level, stokesResidual, stokesTmp, r_u );

                if ( numVCycles == 0 )
                {
Nils Kohl's avatar
Nils Kohl committed
666
667
                   WALBERLA_LOG_INFO_ON_ROOT(
                       walberla::format( "[Uzawa] iter %3d | residual: %10.5e | initial ", 0, vCycleResidualULast ) );
668
                }
Nils Kohl's avatar
Nils Kohl committed
669

Nils Kohl's avatar
Nils Kohl committed
670
                auto reductionRateU = r_u / vCycleResidualULast;
Nils Kohl's avatar
Nils Kohl committed
671

Nils Kohl's avatar
Nils Kohl committed
672
                vCycleResidualULast = r_u;
673

Nils Kohl's avatar
Nils Kohl committed
674
675
                numVCycles++;
                averageResidualReductionU += reductionRateU;
Nils Kohl's avatar
Nils Kohl committed
676

Nils Kohl's avatar
Nils Kohl committed
677
678
679
                if ( verbose )
                {
                   WALBERLA_LOG_INFO_ON_ROOT( walberla::format(
680
                       "[Uzawa] iter %3d | residual: %10.5e | reduction: %10.5e ", numVCycles, r_u, reductionRateU ) );
Nils Kohl's avatar
Nils Kohl committed
681
                }
682

683
684
685
686
687
688
                if ( r_u / initialResiudalU < solverInfo.stokesRelativeResidualUTolerance )
                {
                   WALBERLA_LOG_INFO_ON_ROOT( "[Uzawa] reached relative residual threshold" )
                   return true;
                }

Nils Kohl's avatar
Nils Kohl committed
689
690
                if ( r_u < solverInfo.stokesAbsoluteResidualUTolerance )
                {
691
                   WALBERLA_LOG_INFO_ON_ROOT( "[Uzawa] reached absolute residual threshold" )
Nils Kohl's avatar
Nils Kohl committed
692
693
                   return true;
                }
694

Nils Kohl's avatar
Nils Kohl committed
695
696
                if ( reductionRateU > 0.8 )
                {
697
                   WALBERLA_LOG_INFO_ON_ROOT( "[Uzawa] reached convergence rate threshold" )
Nils Kohl's avatar
Nils Kohl committed
698
699
                   return true;
                }
700

Nils Kohl's avatar
Nils Kohl committed
701
702
                return false;
             };
703

Nils Kohl's avatar
Nils Kohl committed
704
705
         stokesSolver = std::make_shared< SolverLoop< StokesOperator > >(
             multigridSolver, solverInfo.stokesMaxNumIterations, stopIterationCallback );
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
      }
      else
      {
         auto fmgSolver = std::make_shared< FullMultigridSolver< StokesOperator > >(
             storage, multigridSolver, prolongationOperator, minLevel, level, 3 );
         stokesSolver = fmgSolver;
      }
   }
   else
   {
      WALBERLA_ABORT( "Invalid solver type." );
   }

   std::shared_ptr< Solver< UnsteadyDiffusionOperator > > diffusionLinearSolver;

Nils Kohl's avatar
Nils Kohl committed
721
722
723
724
   UnsteadyDiffusion< ScalarFunction, UnsteadyDiffusionOperator, LaplaceOperator, MassOperatorVelocity > diffusionSolver(
       storage, minLevel, level, diffusionLinearSolver );

   if ( solverInfo.diffusionSolverType == DiffusionSolverType::PETSC_MINRES )
725
   {
Nils Kohl's avatar
Nils Kohl committed
726
      auto internalDiffusionSolver = std::make_shared< PETScMinResSolver< UnsteadyDiffusionOperator > >(
727
          storage, level, 1e-30, solverInfo.diffusionAbsoluteResidualUTolerance, solverInfo.diffusionMaxNumIterations );
728
729
730
      internalDiffusionSolver->reassembleMatrix( true );
      diffusionLinearSolver = internalDiffusionSolver;
   }
Nils Kohl's avatar
Nils Kohl committed
731
   else if ( solverInfo.diffusionSolverType == DiffusionSolverType::HYTEG_CG )
732
   {
Nils Kohl's avatar
Nils Kohl committed
733
734
      auto internalDiffusionSolver = std::make_shared< CGSolver< UnsteadyDiffusionOperator > >(
          storage, minLevel, level, solverInfo.diffusionMaxNumIterations, solverInfo.diffusionAbsoluteResidualUTolerance );
735
736
737
      internalDiffusionSolver->setPrintInfo( verbose );
      diffusionLinearSolver = internalDiffusionSolver;
   }
Nils Kohl's avatar
Nils Kohl committed
738
739
740
741
742
   else if ( solverInfo.diffusionSolverType == DiffusionSolverType::HYTEG_GMG )
   {
      WALBERLA_ABORT( "Somethings not working here for 3D" );
      typedef CGSolver< UnsteadyDiffusionOperator >                 CoarseGridSolver_T;
      typedef GeometricMultigridSolver< UnsteadyDiffusionOperator > GMGSolver_T;
Nils Kohl's avatar
Nils Kohl committed
743

Nils Kohl's avatar
Nils Kohl committed
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
      auto coarseGridSolver = std::make_shared< CoarseGridSolver_T >( storage, minLevel, level );
      auto smoother = std::make_shared< WeightedJacobiSmoother< UnsteadyDiffusionOperator > >( storage, minLevel, level, 0.66 );
      auto prolongationOperator = std::make_shared< P2toP2QuadraticProlongation >();
      auto restrictionOperator  = std::make_shared< P2toP2QuadraticRestriction >();
      auto gmgSolver            = std::make_shared< GMGSolver_T >(
          storage, smoother, coarseGridSolver, restrictionOperator, prolongationOperator, minLevel, level, 3, 3 );

      auto stopIterationCallback = [&]( const UnsteadyDiffusionOperator&,
                                        const ScalarFunction& _u,
                                        const ScalarFunction& _b,
                                        const uint_t ) {
         real_t r_c;

         calculateDiffusionResidual( diffusionSolver, diffusionOperator, L, MVelocity, _u, cOld, _b, cTmp, cTmp2, level, r_c );

         auto reductionRateC = r_c / vCycleResidualCLast;

         vCycleResidualCLast = r_c;

         if ( verbose )
         {
            WALBERLA_LOG_INFO_ON_ROOT(
                walberla::format( "[Diffusion GMG] residual: %10.5e | reduction: %10.5e", r_c, reductionRateC ) );
         }

         if ( r_c < solverInfo.diffusionAbsoluteResidualUTolerance )
         {
            return true;
         }

         if ( reductionRateC > 0.8 )
         {
            return true;
         }

         return false;
      };

      diffusionLinearSolver = std::make_shared< SolverLoop< UnsteadyDiffusionOperator > >( gmgSolver, 20, stopIterationCallback );
   }

   diffusionSolver.setSolver( diffusionLinearSolver );
Nils Kohl's avatar
Nils Kohl committed
786
787
788
789
790
791
792
793
794
795

   real_t timeTotal = 0;
   real_t vMax      = 0;
   real_t vRms      = 0;
   real_t residualU = 0;

   real_t timeStepTotal = 0;
   real_t timeStokes    = 0;
   real_t timeMMOC      = 0;
   real_t timeDiffusion = 0;
796
   real_t timeVTK       = 0;
Nils Kohl's avatar
Nils Kohl committed
797

798
   hyteg::VTKOutput vtkOutput( outputInfo.outputDirectory, outputInfo.outputBaseName, storage, outputInfo.vtkInterval );
799
   vtkOutput.setVTKDataFormat( vtk::DataFormat::BINARY );
Nils Kohl's avatar
Nils Kohl committed
800

801
   if ( outputInfo.vtkOutputVertexDoFs )
Nils Kohl's avatar
Nils Kohl committed
802
803
804
805
806
807
808
809
   {
      vtkOutput.add( c.getVertexDoFFunction() );
   }
   else
   {
      vtkOutput.add( c );
   }

810
   if ( outputInfo.vtkOutputVelocity )
811
   {
812
      if ( outputInfo.vtkOutputVertexDoFs )
Nils Kohl's avatar
Nils Kohl committed
813
      {
814
815
816
         vtkOutput.add( u.uvw[0].getVertexDoFFunction() );
         vtkOutput.add( u.uvw[1].getVertexDoFFunction() );
         vtkOutput.add( u.uvw[2].getVertexDoFFunction() );
Nils Kohl's avatar
Nils Kohl committed
817
818
819
820
821
      }
      else
      {
         vtkOutput.add( u );
      }
822
   }
Nils Kohl's avatar
Nils Kohl committed
823
824
   else
   {
825
      if ( outputInfo.vtkOutputVertexDoFs )
Nils Kohl's avatar
Nils Kohl committed
826
827
828
829
830
831
832
      {
         vtkOutput.add( uMagnitudeSquared.getVertexDoFFunction() );
      }
      else
      {
         vtkOutput.add( uMagnitudeSquared );
      }
Nils Kohl's avatar
Nils Kohl committed
833
   }
Nils Kohl's avatar
Nils Kohl committed
834
835
836
837
838
839

   WALBERLA_LOG_INFO_ON_ROOT( "" );
   printFunctionAllocationInfo( *storage, 1 );
   WALBERLA_LOG_INFO_ON_ROOT( "" );
   printCurrentMemoryUsage();
   WALBERLA_LOG_INFO_ON_ROOT( "" );
Nils Kohl's avatar
Nils Kohl committed
840
841
842
843
844

   timer->stop( "Setup" );

   timer->start( "Simulation" );

Nils Kohl's avatar
Nils Kohl committed
845
846
847
848
849
850
851
852
853
   uint_t            timeStep = 0;
   walberla::WcTimer timeStepTimer;

   if ( verbose )
   {
      WALBERLA_LOG_INFO_ON_ROOT( "Initial Stokes solve ..." );
   }

   timeStepTimer.start();
Nils Kohl's avatar
Nils Kohl committed
854

855
   for ( uint_t l = minLevel; l <= level; l++ )
Nils Kohl's avatar
Nils Kohl committed
856
   {
857
858
      MVelocity.apply( c, f.uvw[0], l, All );
      MVelocity.apply( c, f.uvw[1], l, All );
859
860
      if ( storage->hasGlobalCells() )
      {
861
         MVelocity.apply( c, f.uvw[2], l, All );
862
863
      }

864
865
      f.uvw.multElementwise( { f.uvw, outwardNormal.uvw }, l );
      f.uvw.assign( { rayleighNumber }, { f.uvw }, l, All );
Nils Kohl's avatar
Nils Kohl committed
866
867
   }

868
   calculateStokesResiduals( *A, MVelocity, MPressure, u, f, level, stokesResidual, stokesTmp, residualU );
869

Nils Kohl's avatar
Nils Kohl committed
870
   initialResiudalU    = residualU;
871
872
873
874
875
876
   vCycleResidualULast = residualU;

   localTimer.start();
   stokesSolver->solve( *A, u, f, level );
   localTimer.end();
   timeStokes = localTimer.last();
Nils Kohl's avatar
Nils Kohl committed
877

878
   calculateStokesResiduals( *A, MVelocity, MPressure, u, f, level, stokesResidual, stokesTmp, residualU );
Nils Kohl's avatar
Nils Kohl committed
879
880
881

   if ( storage->hasGlobalCells() )
   {
882
      vMax = velocityMaxMagnitude( u.uvw[0], u.uvw[1], u.uvw[2], uTmp, uMagnitudeSquared, level, All );
Nils Kohl's avatar
Nils Kohl committed
883
884
885
   }
   else
   {
886
      vMax = velocityMaxMagnitude( u.uvw[0], u.uvw[1], uTmp, uMagnitudeSquared, level, All );
Nils Kohl's avatar
Nils Kohl committed
887
888
   }

889
   localTimer.start();
890
   if ( outputInfo.vtk )
891
   {
Nils Kohl's avatar
Nils Kohl committed
892
      vtkOutput.write( level );
893
894
895
   }
   localTimer.end();
   timeVTK = localTimer.last();
Nils Kohl's avatar
Nils Kohl committed
896

897
898
899
900
   if ( outputInfo.sphericalTemperatureSlice && timeStep % outputInfo.sphericalTemperatureSliceInterval == 0 )
   {
      if ( verbose )
      {
901
         WALBERLA_LOG_INFO_ON_ROOT( "Evaluating spherical slices ..." );
902
      }
903
904

      evaluateSphericalSliceUV( sliceEvaluationRadius,
905
906
907
908
                                outputInfo.sphericalTemperatureSliceNumMeridians,
                                outputInfo.sphericalTemperatureSliceNumParallels,
                                c,
                                level,
909
                                outputInfo.outputDirectory + "/" + outputInfo.outputBaseName + "_temp_slice_uv_" +
910
                                    std::to_string( timeStep ) + ".csv" );
911
912
913
914
915
916
917

      evaluateSphericalSliceIco( sliceEvaluationRadius,
                                 outputInfo.sphericalTemperatureSliceIcoRefinements,
                                 c,
                                 level,
                                 outputInfo.outputDirectory + "/" + outputInfo.outputBaseName + "_temp_slice_ico_" +
                                     std::to_string( timeStep ) + ".csv" );
918
919
920
921
922
923
      if ( verbose )
      {
         WALBERLA_LOG_INFO_ON_ROOT( "Done." );
      }
   }

Nils Kohl's avatar
Nils Kohl committed
924
925
   timeStepTimer.end();
   timeStepTotal = timeStepTimer.last();
Nils Kohl's avatar
Nils Kohl committed
926

927
   WALBERLA_LOG_INFO_ON_ROOT(
928
       " timestep |           dt |   time total | velocity RMS | velocity max magnitude |   residual u |  total | Stokes |   MMOC |   diff |    VTK |" )
929
   WALBERLA_LOG_INFO_ON_ROOT(
930
       "----------+--------------+--------------+--------------+------------------------+--------------+--------+--------+--------+--------+--------+" )
931
   WALBERLA_LOG_INFO_ON_ROOT(
932
       walberla::format( " %8s | %12s | %12.8f | %12.4f | %22.4f | %12.5e | %6.2f | %6.2f | %6.2f | %6.2f | %6.2f |",
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
                         "initial",
                         "-",
                         timeTotal,
                         vRms,
                         vMax,
                         residualU,
                         timeStepTotal,
                         timeStokes,
                         timeMMOC,
                         timeDiffusion,
                         timeVTK ) )

   db.setVariableEntry( "ts", uint_c( timeStep ) );
   db.setVariableEntry( "sim_time", timeTotal );
   db.setVariableEntry( "run_time_advection", timeMMOC );
   db.setVariableEntry( "run_time_diffusion", timeDiffusion );
   db.setVariableEntry( "run_time_stokes", timeStokes );
   db.setVariableEntry( "run_time_time_step", timeStepTotal );
   db.setVariableEntry( "run_time_vtk", timeVTK );
   db.setVariableEntry( "v_max", vMax );
   db.setVariableEntry( "v_rms", vRms );
   db.setVariableEntry( "dt", real_c( 0 ) );

Nils Kohl's avatar
Nils Kohl committed
956
957
958
959
960
961
962
963
964
   db.setVariableEntry( "initial_residual_u_predictor", initialResiudalU );
   db.setVariableEntry( "num_v_cycles_predictor", numVCycles );
   db.setVariableEntry( "avg_residual_reduction_u_predictor", averageResidualReductionU / real_c( numVCycles ) );
   db.setVariableEntry( "final_residual_u_predictor", vCycleResidualULast );

   db.setVariableEntry( "initial_residual_u_corrector", real_c( 0 ) );
   db.setVariableEntry( "num_v_cycles_corrector", real_c( 0 ) );
   db.setVariableEntry( "avg_residual_reduction_u_corrector", real_c( 0 ) );
   db.setVariableEntry( "final_residual_u_corrector", real_c( 0 ) );
Nils Kohl's avatar
Nils Kohl committed
965

966
967
   db.writeRowOnRoot();

968
969
   timer->stop( "Simulation" );

Nils Kohl's avatar
Nils Kohl committed
970
   while ( timeTotal < simulationTime && timeStep < maxNumTimeSteps )
Nils Kohl's avatar
Nils Kohl committed
971
   {
972
973
      timer->start( "Simulation" );

Nils Kohl's avatar
Nils Kohl committed
974
975
976
977
978
979
980
981
      timeStepTimer.start();

      timeStep++;

      // new time step size

      if ( storage->hasGlobalCells() )
      {
982
         vMax = velocityMaxMagnitude( u.uvw[0], u.uvw[1], u.uvw[2], uTmp, uTmp2, level, All );
Nils Kohl's avatar
Nils Kohl committed
983
984
985
      }
      else
      {
986
         vMax = velocityMaxMagnitude( u.uvw[0], u.uvw[1], uTmp, uTmp2, level, All );
Nils Kohl's avatar
Nils Kohl committed
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
      }

      real_t dt;
      if ( fixedTimeStep )
      {
         dt = dtConstant;
      }
      else
      {
         dt = ( cflMax / vMax ) * hMin;
      }

      // energy

      // advection

      // start value for predictor
1004
      cPr.assign( { 1.0 }, { c }, level, All );
Nils Kohl's avatar
Nils Kohl committed
1005
1006

      // let's just use the current velocity for the prediction
1007
      uLast.assign( { 1.0 }, { u }, level, All );
Nils Kohl's avatar
Nils Kohl committed
1008
1009
1010
1011
1012

      // diffusion

      cPr.interpolate( initialTemperature, level, DirichletBoundary );

1013
      cOld.assign( { 1.0 }, { cPr }, level, All );
Nils Kohl's avatar
Nils Kohl committed
1014
1015
1016

      diffusionOperator.setDt( 0.5 * dt );

1017
1018
1019
1020
1021
      if ( storage->hasGlobalCells() )
      {
         diffusionOperator.getOperator().computeAndStoreLocalElementMatrices();
      }

Nils Kohl's avatar
Nils Kohl committed
1022
1023
1024
      calculateDiffusionResidual(
          diffusionSolver, diffusionOperator, L, MVelocity, cPr, cOld, q, cTmp, cTmp2, level, vCycleResidualCLast );

Nils Kohl's avatar
Nils Kohl committed
1025
1026
1027
1028
1029
1030
      localTimer.start();
      diffusionSolver.step( diffusionOperator, L, MVelocity, cPr, cOld, q, q, level, Inner | NeumannBoundary | FreeslipBoundary );
      localTimer.end();
      timeDiffusion = localTimer.last();

      localTimer.start();
Marcus Mohr's avatar
Marcus Mohr committed
1031
      transport.step( cPr, u.uvw, uLast.uvw, level, All, dt, 1, true );
Nils Kohl's avatar
Nils Kohl committed
1032
1033
1034
1035
1036
1037
1038
      localTimer.end();
      timeMMOC = localTimer.last();

      // diffusion

      cPr.interpolate( initialTemperature, level, DirichletBoundary );

1039
      cOld.assign( { 1.0 }, { cPr }, level, All );
Nils Kohl's avatar
Nils Kohl committed
1040

Nils Kohl's avatar
Nils Kohl committed
1041
1042
1043
      calculateDiffusionResidual(
          diffusionSolver, diffusionOperator, L, MVelocity, cPr, cOld, q, cTmp, cTmp2, level, vCycleResidualCLast );

Nils Kohl's avatar
Nils Kohl committed
1044
1045
1046
1047
1048
1049
1050
      localTimer.start();
      diffusionSolver.step( diffusionOperator, L, MVelocity, cPr, cOld, q, q, level, Inner | NeumannBoundary | FreeslipBoundary );
      localTimer.end();
      timeDiffusion += localTimer.last();

      // Stokes

1051
      for ( uint_t l = minLevel; l <= level; l++ )
Nils Kohl's avatar
Nils Kohl committed
1052
      {
1053
1054
         MVelocity.apply( cPr, f.uvw[0], l, All );
         MVelocity.apply( cPr, f.uvw[1], l, All );
1055
1056
         if ( storage->hasGlobalCells() )
         {
1057
            MVelocity.apply( cPr, f.uvw[2], l, All );
1058
1059
         }

1060
1061
         f.uvw.multElementwise( { f.uvw, outwardNormal.uvw }, l );
         f.uvw.assign( { rayleighNumber }, { f.uvw }, l, All );
Nils Kohl's avatar
Nils Kohl committed
1062
      }
1063

1064
      calculateStokesResiduals( *A, MVelocity, MPressure, u, f, level, stokesResidual, stokesTmp, residualU );
1065

Nils Kohl's avatar
Nils Kohl committed
1066
1067
      vCycleResidualULast       = residualU;
      initialResiudalU          = residualU;
Nils Kohl's avatar
Nils Kohl committed
1068
      averageResidualReductionU = real_c( 0 );
Nils Kohl's avatar
Nils Kohl committed
1069
      numVCycles                = 0;
Nils Kohl's avatar
Nils Kohl committed
1070
1071

      localTimer.start();
1072
      stokesSolver->solve( *A, u, f, level );
Nils Kohl's avatar
Nils Kohl committed
1073
1074
1075
      localTimer.end();
      timeStokes = localTimer.last();

1076
      calculateStokesResiduals( *A, MVelocity, MPressure, u, f, level, stokesResidual, stokesTmp, residualU );
Nils Kohl's avatar
Nils Kohl committed
1077

Nils Kohl's avatar
Nils Kohl committed
1078
1079
1080
1081
1082
      db.setVariableEntry( "initial_residual_u_predictor", initialResiudalU );
      db.setVariableEntry( "num_v_cycles_predictor", numVCycles );
      db.setVariableEntry( "avg_residual_reduction_u_predictor", averageResidualReductionU / real_c( numVCycles ) );
      db.setVariableEntry( "final_residual_u_predictor", vCycleResidualULast );

Nils Kohl's avatar
Nils Kohl committed
1083
1084
1085
1086
1087
1088
1089
1090
      if ( predictorCorrector )
      {
         // energy

         // diffusion

         c.interpolate( initialTemperature, level, DirichletBoundary );

1091
         cOld.assign( { 1.0 }, { c }, level, All );
Nils Kohl's avatar
Nils Kohl committed
1092

Nils Kohl's avatar
Nils Kohl committed
1093
1094
1095
         calculateDiffusionResidual(
             diffusionSolver, diffusionOperator, L, MVelocity, c, cOld, q, cTmp, cTmp2, level, vCycleResidualCLast );

Nils Kohl's avatar
Nils Kohl committed
1096
1097
1098
1099
1100
1101
1102
1103
1104
         localTimer.start();
         diffusionSolver.step(
             diffusionOperator, L, MVelocity, c, cOld, q, q, level, Inner | NeumannBoundary | FreeslipBoundary );
         localTimer.end();
         timeDiffusion += localTimer.last();

         // advection

         localTimer.start();
Marcus Mohr's avatar
Marcus Mohr committed
1105
         transport.step( c, u.uvw, uLast.uvw, level, All, dt, 1, true );
Nils Kohl's avatar
Nils Kohl committed
1106
1107
1108
1109
1110
1111
1112
         localTimer.end();
         timeMMOC += localTimer.last();

         // diffusion

         c.interpolate( initialTemperature, level, DirichletBoundary );

1113
         cOld.assign( { 1.0 }, { c }, level, All );
Nils Kohl's avatar
Nils Kohl committed
1114

Nils Kohl's avatar
Nils Kohl committed
1115
1116
1117
         calculateDiffusionResidual(
             diffusionSolver, diffusionOperator, L, MVelocity, c, cOld, q, cTmp, cTmp2, level, vCycleResidualCLast );

Nils Kohl's avatar
Nils Kohl committed
1118
1119
1120
1121
1122
1123
1124
1125
         localTimer.start();
         diffusionSolver.step(
             diffusionOperator, L, MVelocity, c, cOld, q, q, level, Inner | NeumannBoundary | FreeslipBoundary );
         localTimer.end();
         timeDiffusion += localTimer.last();

         // Stokes

1126
         for ( uint_t l = minLevel; l <= level; l++ )
Nils Kohl's avatar
Nils Kohl committed
1127
         {
1128
1129
            MVelocity.apply( c, f.uvw[0], l, All );
            MVelocity.apply( c, f.uvw[1], l, All );
1130
1131
            if ( storage->hasGlobalCells() )
            {
1132
               MVelocity.apply( c, f.uvw[2], l, All );
1133
1134
            }

1135
1136
            f.uvw.multElementwise( { f.uvw, outwardNormal.uvw }, l );
            f.uvw.assign( { rayleighNumber }, { f.uvw }, l, All );
Nils Kohl's avatar
Nils Kohl committed
1137
         }
1138

1139
         calculateStokesResiduals( *A, MVelocity, MPressure, u, f, level, stokesResidual, stokesTmp, residualU );
1140

Nils Kohl's avatar
Nils Kohl committed
1141
1142
         vCycleResidualULast       = residualU;
         initialResiudalU          = residualU;
Nils Kohl's avatar
Nils Kohl committed
1143
         averageResidualReductionU = real_c( 0 );
Nils Kohl's avatar
Nils Kohl committed
1144
         numVCycles                = 0;
Nils Kohl's avatar
Nils Kohl committed
1145
1146

         localTimer.start();
1147
         stokesSolver->solve( *A, u, f, level );
Nils Kohl's avatar
Nils Kohl committed
1148
1149
1150
         localTimer.end();
         timeStokes += localTimer.last();

1151
         calculateStokesResiduals( *A, MVelocity, MPressure, u, f, level, stokesResidual, stokesTmp, residualU );
Nils Kohl's avatar
Nils Kohl committed
1152
1153
1154
1155
1156

         db.setVariableEntry( "initial_residual_u_corrector", initialResiudalU );
         db.setVariableEntry( "num_v_cycles_corrector", numVCycles );
         db.setVariableEntry( "avg_residual_reduction_u_corrector", averageResidualReductionU / real_c( numVCycles ) );
         db.setVariableEntry( "final_residual_u_corrector", vCycleResidualULast );
Nils Kohl's avatar
Nils Kohl committed
1157
1158
1159
1160
      }
      else
      {
         // use predicted value
1161
         c.assign( { 1.0 }, { cPr }, level, All );
Nils Kohl's avatar
Nils Kohl committed
1162
1163
1164
1165
1166

         db.setVariableEntry( "initial_residual_u_corrector", real_c( 0 ) );
         db.setVariableEntry( "num_v_cycles_corrector", real_c( 0 ) );
         db.setVariableEntry( "avg_residual_reduction_u_corrector", real_c( 0 ) );
         db.setVariableEntry( "final_residual_u_corrector", real_c( 0 ) );
Nils Kohl's avatar
Nils Kohl committed
1167
1168
1169
1170
      }

      timeTotal += dt;

Nils Kohl's avatar
Nils Kohl committed
1171
      vRms = velocityRMS( u, stokesTmp, MVelocity, domainInfo.domainVolume(), level );
Nils Kohl's avatar
Nils Kohl committed
1172

Nils Kohl's avatar
Nils Kohl committed
1173
1174
      if ( storage->hasGlobalCells() )
      {
1175
         vMax = velocityMaxMagnitude( u.uvw[0], u.uvw[1], u.uvw[2], uTmp, uMagnitudeSquared, level, All );
Nils Kohl's avatar
Nils Kohl committed
1176
1177
1178
      }
      else
      {
1179
         vMax = velocityMaxMagnitude( u.uvw[0], u.uvw[1], uTmp, uMagnitudeSquared, level, All );
Nils Kohl's avatar
Nils Kohl committed
1180
1181
      }

1182
      localTimer.start();
1183
      if ( outputInfo.vtk )
Nils Kohl's avatar
Nils Kohl committed
1184
1185
1186