The TS-GitLab will have to shut down towards the end of the year — please think about migrating your projects to GitLab.com or GitHub.
(This is still a very early message, meant to keep you informed. There will be more urgent ones in the future.)

transport.cc 16.1 KB
Newer Older
1
#include <dune/dorie/model/transport/transport.hh>
2

3 4 5 6
namespace Dune{
namespace Dorie{

template<typename Traits>
7
ModelTransport<Traits>::ModelTransport(
8 9
  Dune::ParameterTree&              _inifile,
  const GridCreator&                _grid_creator,
10
  Dune::MPIHelper&                  _helper
11
) : ModelBase(log_transport,
12 13
              _inifile.get<std::string>("output.logLevel"),
              _helper)
14
  , helper(_helper)
15
  , inifile(_inifile)
16
  , time_step_ctrl(inifile, this->_log)
17 18 19
  , output_type(_inifile.get<bool>("output.asciiVtk") ? 
      Dune::VTK::OutputType::ascii : Dune::VTK::OutputType::base64)
  , grid(_grid_creator.grid())
20
  , gv(grid->leafGridView())
21 22
  , mbe_slop(estimate_mbe_entries<typename MBE::size_type>(Traits::dim,Traits::GridGeometryType))
  , mbe_tlop(1)
23
  , time_steps(0)
24
  , _courant_number(inifile.get<RF>("numerics.courant"))
25
  , enable_fluxrc(_inifile.get<bool>("fluxReconstruction.enable"))
26
{
27 28 29
  // initialize time
  this->_time = time_step_ctrl.get_time_begin();

30 31 32 33 34 35 36 37 38 39 40 41 42 43
  // --- Output Policy ---
  std::string output_policy_str 
      = inifile.get<std::string>("output.policy");

  if (output_policy_str == "endOfRichardsStep") {
    this->set_policy(OutputPolicy::EndOfRichardsStep);
  } else if (output_policy_str == "endOfTransportStep") {
    this->set_policy(OutputPolicy::EndOfTransportStep);
  } else if (output_policy_str == "none") {
    this->set_policy(OutputPolicy::None);
  } else {
    this->_log->error("Invalid output policy '{}'", output_policy_str);
    DUNE_THROW(NotImplemented,"Invalid output policy: " << output_policy_str);
  }
44

45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
  // --- Solution Check Policy --- //
  const auto check_negative_policy_str
    = inifile.get<std::string>("solutionCheck.negativeConcentration");
  if (check_negative_policy_str == "pass") {
    check_negative_policy = CheckNegativePolicy::Pass;
  }
  else if (check_negative_policy_str == "warn") {
    check_negative_policy = CheckNegativePolicy::Warn;
  }
  else if (check_negative_policy_str == "retry") {
    check_negative_policy = CheckNegativePolicy::Retry;
  }
  else {
    this->_log->error("Invalid policy for negative solute concentrations: {}",
                      check_negative_policy_str);
    DUNE_THROW(IOError, "Invalid policy for negative concentrations");
  }

63 64 65
  // initialize time
  this->_time = time_step_ctrl.get_time_begin();

66
  // --- Grid Function Space ---
67
  this->_log->trace("Setting up GridFunctionSpace");
68 69 70 71 72
  gfs = std::make_shared<GFS>(GFSHelper::create(gv));
  gfs->name("solute");
  cc = std::make_unique<CC>();
  Dune::PDELab::constraints(*gfs,*cc,false);

73
  // --- Create the new parameter class
74
  auto element_map = _grid_creator.element_index_map();
75 76
  sparam = std::make_shared<SoluteParameters>(inifile, grid, element_map);

77
  // --- Operator Helper Classes ---
78 79 80 81
  const auto& boundary_index_map = _grid_creator.boundary_index_map();
  sboundary = std::make_shared<SoluteBoundary>(inifile,
                                               boundary_index_map,
                                               this->_log);
82
  sinitial = SoluteInitialFactory::create(inifile, gv, this->_log);
83

84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
	// Read local operator settings
	namespace OP = Dune::Dorie::Operator;
	const auto settings = OP::read_transport_operator_settings(inifile);

  // --- Local Operators ---  
  if constexpr (order>0)
  {
    this->_log->debug("Setting up local grid operators: DG method");

    namespace OP = Dune::Dorie::Operator;
    const auto method = std::get<OP::TransportDGMethod>(settings);
    const auto upwinding = std::get<OP::TransportUpwinding>(settings);
    const auto weights = std::get<OP::TransportDGWeights>(settings);

    slop = std::make_unique<SLOP>(inifile, sparam, sboundary, 
                                  method, upwinding, 
                                  weights);
  }
  else {
    this->_log->debug("Setting up local grid operators: FV method");
    slop = std::make_unique<SLOP>(sparam, sboundary/*, upwinding*/);
  }

  tlop = std::make_unique<TLOP>();

109 110 111 112 113 114
  // --- Solution Vectors and Initial Condition ---
  const RF ini_value = inifile.get<RF>("initial.value",0.);
  u = std::make_shared<U>(*gfs,ini_value);
  
  Dune::PDELab::interpolate(*sinitial,*gfs,*u);

115
  // Select time stepping method
116 117 118 119 120 121 122 123 124 125
  std::string ts_method = _inifile.get<std::string>("numerics.timestepMethod");
  if (ts_method== "explicit_euler")
    ts_param = std::make_unique<Dune::PDELab::ExplicitEulerParameter<RF>>();
  else if (ts_method== "implicit_euler")
    ts_param = std::make_unique<Dune::PDELab::ImplicitEulerParameter<RF>>();
  else if (ts_method== "alex2")
    ts_param = std::make_unique<Dune::PDELab::Alexander2Parameter<RF>>();
  else
    DUNE_THROW(Dune::NotImplemented,"Time method not supported!");

126 127 128 129 130 131 132 133 134 135 136 137
  if (not ts_param->implicit())
  {
    if (FloatCmp::le(_courant_number, 0.0)) {
      this->_log->error("Courant number must be larger than zero");
      DUNE_THROW(IOError, "Invalid courant number");
    }
    else if (FloatCmp::gt(_courant_number, 1.0)) {
      this->_log->warn("Courant number larger than 1.0! Solver might not "
                       "converge, leading to instabilities in the solution.");
    }
  }

138
  // --- Writer Setup --- //
139 140
  if (output_policy() != OutputPolicy::None)
  {
141 142 143 144 145 146 147
    const int subsamling_lvl =
      _inifile.get<int>("output.subsamplingLevel", 0);
    const auto output_path = inifile.get<std::string>("output.outputPath");
		this->_log->debug("Creating a VTK writer with subsampling level: {}",
						  subsamling_lvl);

    Dorie::create_directories(output_path);
148 149 150
    const auto subsamling_intervals = Dune::refinementLevels(subsamling_lvl);
    auto sub_vtk = std::make_shared<Dune::SubsamplingVTKWriter<GV>>(gv,
      subsamling_intervals);
151

152 153
    vtkwriter = std::make_unique<Writer>(sub_vtk,
              inifile.get<std::string>("output.fileName"),
154
              output_path,
155
              "./");
156 157
  }

158 159
  // --- Operator Setup ---
  operator_setup();
160

161
  this->_log->info("Setup complete");
162 163 164 165
}


template<typename Traits>
166
void ModelTransport<Traits>::operator_setup()
167
{
168 169 170 171 172 173 174 175
  auto dof_count = gfs->globalSize();
  this->_log->debug("Setting up grid operators and solvers");
  if (helper.size() > 1) {
    this->_log->debug("  DOF of this process: {}", dof_count);
  }
  dof_count = gv.comm().sum(dof_count);
  this->_log->debug("  Total number of DOF: {}", dof_count);

176
  // --- Grid Operators ---
177
  go0 = std::make_shared<GO0>(*gfs,*cc,*gfs,*cc,*slop,mbe_slop);
178
  go1 = std::make_unique<GO1>(*gfs,*cc,*gfs,*cc,*tlop,mbe_tlop);
179 180 181 182
  if (ts_param->implicit())
    implicit_igo = std::make_unique<ImplicitIGO>(*go0,*go1);
  else
    explicit_igo = std::make_unique<ExplicitIGO>(*go0,*go1);
183 184

  // --- Solvers ---
185 186 187 188 189 190 191
  // Initialize helper spaces for DG linear solver only
  if constexpr (order > 0)
  {
    lsgfs = std::make_unique<LSGFS>(LSGFSHelper::create(gv));
    lscc = std::make_unique<LSCC>();
  }

192
  if (ts_param->implicit())
193 194 195 196 197 198 199 200 201
  {
    if constexpr (order == 0)
      implicit_ls = std::make_unique<ImplicitLS>(*gfs, 1000, 0, true, true);
    else
    {
      implicit_ls = std::make_unique<ImplicitLS>(
        *implicit_igo, *cc, *lsgfs, *lscc, 1000, 0, true, true);
    }
  }
202
  else
203 204 205 206 207 208 209 210 211
  {
    if constexpr (order == 0)
      explicit_ls = std::make_unique<ExplicitLS>(*gfs, 1000, 0, true, true);
    else
    {
      explicit_ls = std::make_unique<ExplicitLS>(
        *explicit_igo, *cc, *lsgfs, *lscc, 1000, 0, true, true);
    }
  }
212 213

  // --- Time Step Operators ---
214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231
  if (ts_param->implicit()){
    // Set linear problem solver verbosity to zero
    inifile["solverParameters.verbosity"] = "0";

    implicit_slps = std::make_unique<ImplicitSLPS>(
      *implicit_igo, *implicit_ls, inifile.sub("solverParameters")
    );
    implicit_osm = std::make_unique<ImplicitOSM>(
      *ts_param, *implicit_igo, *implicit_slps
    );
    implicit_osm->setVerbosityLevel(0);
  }
  else {
    explicit_osm = std::make_unique<ExplicitOSM>(
      *ts_param, *explicit_igo, *explicit_ls
    );
    explicit_osm->setVerbosityLevel(0);
  }
232 233 234 235 236

  gfs->update();
}

template<typename Traits>
237
void ModelTransport<Traits>::step()
238
{
239 240
  // long time step log for level 'debug'
  if (this->_log->should_log(spdlog::level::debug)) {
241
    this->_log->info("Time Step {}:",time_steps);
242
  }
243 244 245 246 247 248 249 250

  if (!water_flux_gf) {
    this->_log->error("Water flux grid function required");
    DUNE_THROW(Dune::InvalidStateException,
               "Pointer to water_flux_gf is invalid!");
  }
  if (!water_content_gf) {
    this->_log->error("Water content grid function required");
251
    DUNE_THROW(Dune::InvalidStateException,
252 253
               "Pointer to water_content_gf is invalid!"); 
  }
254

255
  // Suggest times step for explicit methods. 
256 257 258 259 260 261 262 263
  // NOTE: Do not futher increase time step on failure!
  if (not ts_param->implicit())
  {
    const auto cfl = Dune::Dorie::cfl_condition(*water_flux_gf,
                                                *water_content_gf,
                                                order);
    this->_log->trace("CFL condition: {}s", cfl);
    time_step_ctrl.suggest_time_step(_courant_number * cfl);
264
  }
265

266
  // get time step suggestion from boundary conditions
267 268
  sboundary->set_time(this->_time);
  const auto dt_suggestion = sboundary->max_time_step(this->_time);
269 270 271 272
  time_step_ctrl.suggest_time_step(dt_suggestion);

  // Allocate DOF vector for next time step
  auto unext = std::make_shared<U>(*u);
273

274
  // Obtain time variables
275
  // NOTE: Communicate minimal time step between processors!
276 277
  const auto time = this->_time;
  auto dt = time_step_ctrl.get_time_step(time);
278
  dt = gv.comm().min(dt);
279 280 281 282

  // Solver loop
  // NOTE: Runs infinitely until solution is computed and valid (causes break)
  //       or exception is thrown.
283
  while (true)
284
  {
285 286 287 288 289
    // Long time step log for level 'debug'
    if (this->_log->should_log(spdlog::level::debug)) {
      this->_log->debug("  Time {:.2e} + {:.2e} -> {:.2e}",
                time, dt, time+dt);
    }
290

291 292 293
    try
    {
      // Apply the solver
294
      dwarn.push(false);
295 296 297 298
      if (ts_param->implicit())
        implicit_osm->apply(time, dt, *u, *unext);
      else
        explicit_osm->apply(time, dt, *u, *unext);
299 300
      dwarn.pop();

301 302 303 304 305 306 307 308 309 310 311
      // Check the solution
      if (auto [valid, msg] = check_solution(*u, *unext);
          valid)
      {
        // Solution computed and valid: leave loop
        break;
      }
      else
      {
        this->_log->debug("  Invalid solution");
        this->_log->trace("  Error during solution check: {}", msg);
312
      }
313 314
    }
    catch (Dune::ISTLError &e){
315 316
      this->_log->debug("  Solver not converged");
      this->_log->trace("  Error in linear solver: {}", e.what());
317 318
    }
    catch (Dune::Exception &e){
319
      this->_log->error("  Unexpected error in solver: {}", e.what());
320
      DUNE_THROW(Dune::Exception, "Critical error in solver");
321
    }
322 323 324 325
    // rethrow anything else
    catch (...) {
      throw;
    }
326

327
    // abort if minimal time step is already reached
328
    if (Dune::FloatCmp::le(dt, time_step_ctrl.get_time_step_min())) {
329
      this->_log->error("Computing a solution for the minimal time "
330
                        "step failed");
331 332
      DUNE_THROW(Exception, "No solution for minimal time step");
    }
333

334
    // Reduce time step and try again
335
    // NOTE: Communicate minimal time step between processes!
336
    time_step_ctrl.decrease_time_step();
337
    dt = time_step_ctrl.get_time_step(time);
338
    dt = gv.comm().min(dt);
339
  } // while (true)
340

341 342 343 344 345 346 347 348 349
  // Short time step log for level 'info' (after success)
  if (not this->_log->should_log(spdlog::level::debug)) {
    this->_log->info("Time Step {}: {:.2e} + {:.2e} -> {:.2e}",
                     time_steps,
                     time,
                     dt,
                     time+dt);
  }

350 351 352 353 354 355 356 357 358 359 360 361
  // Print computing times in trace log
  if (ts_param->implicit())
  {
    auto result = implicit_osm->result();
    this->_log->trace("  Matrix assembly: {:.2e}s (total), {:.2e}s (success)",
                      result.total.assembler_time,
                      result.successful.assembler_time);
    this->_log->trace("  Linear solver:   {:.2e}s (total), {:.2e}s (success)",
                      result.total.linear_solver_time,
                      result.successful.linear_solver_time);
  }

362 363 364 365 366 367
  // Overwrite last solution with new solution
  u = unext;

  // Update time and time step
  this->_time += dt;
  time_step_ctrl.increase_time_step();
368
  time_steps++;
369

370 371
  cache_solute_flux_gf_rt.reset();

372
  if (this->output_policy() == OutputPolicy::EndOfTransportStep)
373 374 375 376
    write_data();
}

template<typename Traits>
377
void ModelTransport<Traits>::write_data () const
378
{
379 380
  if (vtkwriter)
  {   
381 382
    auto peclet = std::make_shared<const GFPeclet>(gv, sparam, water_flux_gf, water_content_gf);
    auto d_hd = std::make_shared<const GFEffectiveHydrodynamicDispersion>(gv, sparam, water_flux_gf, water_content_gf);
383

384
    if (inifile.get<bool>("output.vertexData")) {
385 386 387
      vtkwriter->addVertexData(get_solute(),"solute");
      vtkwriter->addVertexData(get_total_solute(),"total_solute");
      vtkwriter->addVertexData(peclet,"micro_peclet");
388 389
      
      if (inifile.get<bool>("output.writeDispersionTensor"))
390
        vtkwriter->addVertexData(d_hd,"eff_hd_dispersion");
391

392
      if constexpr (enable_rt_engine)
393 394
        if (enable_fluxrc) {
          auto RT_name = "flux_RT" + std::to_string(flux_order);
395
          vtkwriter->addVertexData(get_solute_flux_reconstructed(),RT_name);
396
        }
397
    } else {
398 399 400
      vtkwriter->addCellData(get_solute(),"solute");
      vtkwriter->addCellData(get_total_solute(),"total_solute");
      vtkwriter->addCellData(peclet,"micro_peclet");
401 402
      
      if (inifile.get<bool>("output.writeDispersionTensor"))
403
        vtkwriter->addCellData(d_hd,"eff_hd_dispersion");
404
      
405
      if constexpr (enable_rt_engine)
406 407
        if (enable_fluxrc) {
          auto RT_name = "flux_RT" + std::to_string(flux_order);
408
          vtkwriter->addCellData(get_solute_flux_reconstructed(),RT_name);
409
        }
410 411 412
    }

    try{
413 414
      this->_log->trace("Writing solution at time {:.2e}", this->_time);
      vtkwriter->write(this->_time, output_type);
415
    }
416 417 418
    catch (Dune::Exception& e) {
      this->_log->error("Writing VTK output failed: {}", e.what());
      DUNE_THROW(Dune::IOError, "Cannot write VTK output!");
419 420
    }
    catch(...){
421 422
      this->_log->error("Writing VTK output failed for unknown reason");
      DUNE_THROW(Dune::IOError, "Cannot write VTK output!");
423 424
    }
    vtkwriter->clear();
425 426 427 428
  } else {
    this->_log->error("Calling 'write_data' on object without VTK writer");
    DUNE_THROW(Dune::InvalidStateException, "No vtk writer configured!");
  }
429 430
}

431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451
template<typename Traits>
void ModelTransport<Traits>::mark_grid()
{
  DUNE_THROW(NotImplemented,
    "Dorie does not implement adaptivity for this simulation object!");
}

template<typename Traits>
void ModelTransport<Traits>::adapt_grid()
{
  DUNE_THROW(NotImplemented,
   "Dorie does not implement adaptivity for this simulation object!");
}


template<typename Traits>
void ModelTransport<Traits>::post_adapt_grid()
{
  operator_setup();
}

452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494
template<typename Traits>
std::pair<bool, std::string> ModelTransport<Traits>::check_solution(
  const U& u_before,
  const U& u_after
) const
{
  // Check for negative values
  if (check_negative_policy != CheckNegativePolicy::Pass)
  {
    // Compute tolerance from largest solute concentration and user input
    const auto max = *std::max_element(u_after.begin(), u_after.end());
    const auto eps = max * inifile.get<RF>("solutionCheck.negativeTol");

    // Check if values are negative and record minimum value encountered
    auto min = std::numeric_limits<RF>::max();
    auto check_negative = [eps, &min](const auto val)
                          {
                            min = std::min(min, val);
                            return val < -eps;
                          };

    // Perform the check
    if (std::any_of(u_after.begin(), u_after.end(), check_negative))
    {
      // NOTE: fmt library included via spdlog
      const std::string msg = fmt::format("Negative value in solution: {:.2e}",
                                          min);

      // Issue warning, but pass check
      if (check_negative_policy == CheckNegativePolicy::Warn) {
        this->_log->warn(msg);
      }
      // Do not pass check
      else {
        return std::make_pair(false, msg);
      }
    }
  }

  // Default: Pass without message
  return std::make_pair(true, "");
}

495
} // namespace Dorie
496
} // namespace Dune