1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 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
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
// *****************************************************************************
/*!
  \file      src/PDE/MultiMat/Problem/UserDefined.cpp
  \copyright 2012-2015 J. Bakosi,
             2016-2018 Los Alamos National Security, LLC.,
             2019-2021 Triad National Security, LLC.
             All rights reserved. See the LICENSE file for details.
  \brief     Problem configuration for the compressible flow equations
  \details   This file defines a Problem policy class for the compressible flow
    equations, defined in PDE/MultiMat/MultiMat.h. See PDE/MultiMat/Problem.h
    for general requirements on Problem policy classes for MultiMat.
*/
// *****************************************************************************

#include <limits>

#include "UserDefined.hpp"
#include "Inciter/InputDeck/InputDeck.hpp"
#include "FieldOutput.hpp"
#include "MultiMat/MultiMatIndexing.hpp"

namespace inciter {

extern ctr::InputDeck g_inputdeck;

} // ::inciter

using inciter::MultiMatProblemUserDefined;

tk::InitializeFn::result_type
MultiMatProblemUserDefined::initialize( ncomp_t ncomp,
                                        const std::vector< EOS >& mat_blk,
                                        tk::real,
                                        tk::real,
                                        tk::real,
                                        tk::real )
// *****************************************************************************
//! Set initial conditions
//! \param[in] ncomp Number of scalar components in this PDE system
//! \return Values of all components
//! \note The function signature must follow tk::InitializeFn
// *****************************************************************************
{
  tk::InitializeFn::result_type s( ncomp, 0.0 );

  auto nmat = g_inputdeck.get< tag::param, eq, tag::nmat >();<--- Variable 'nmat' is assigned a value that is never used.
  const auto& solidx = g_inputdeck.get< tag::param, tag::multimat,
    tag::matidxmap >().template get< tag::solidx >();

  // Set background ICs
  const auto& ic = g_inputdeck.get< tag::param, eq, tag::ic >();
  const auto& bgmatid = ic.get< tag::materialid >();
  const auto& bgvelic = ic.get< tag::velocity >();
  const auto& bgpreic = ic.get< tag::pressure >();
  const auto& bgtempic = ic.get< tag::temperature >();

  Assert( bgtempic.size() > 0, "No background temperature IC" );
  Assert( bgpreic.size() > 0, "No background pressure IC" );

  auto alphamin = 1.0e-12;

  // initialize background material states
  for (std::size_t k=0; k<nmat; ++k) {
    if (k == bgmatid.at(0)-1) {
      s[volfracIdx(nmat,k)] = 1.0 - (static_cast< tk::real >(nmat-1))*alphamin;
    }
    else {
      s[volfracIdx(nmat,k)] = alphamin;
    }
  }

  tk::real u = bgvelic[0];
  tk::real v = bgvelic[1];
  tk::real w = bgvelic[2];

  auto rb = 0.0;
  for (std::size_t k=0; k<nmat; ++k) {
    // density
    auto rhok = mat_blk[k].compute< EOS::density >(bgpreic[0], bgtempic[0]);
    // partial density
    s[densityIdx(nmat,k)] = s[volfracIdx(nmat,k)] * rhok;
    // deformation gradients
    std::array< std::array< tk::real, 3 >, 3 > g;
    if (solidx[k] > 0) {
      for (std::size_t i=0; i<3; ++i) {
        for (std::size_t j=0; j<3; ++j) {
          if (i==j) g[i][j] = 1.0;
          else g[i][j] = 0.0;
          s[deformIdx(nmat,solidx[k],i,j)] = s[volfracIdx(nmat,k)]*g[i][j];
        }
      }
    }
    else {
      g = {{}};
    }
    // total specific energy
    s[energyIdx(nmat,k)] = s[volfracIdx(nmat,k)] *
      mat_blk[k].compute< EOS::totalenergy >(rhok, u, v, w, bgpreic[0],
      g);
    // bulk density
    rb += s[densityIdx(nmat,k)];
  }

  // bulk momentum
  s[momentumIdx(nmat,0)] = rb * u;
  s[momentumIdx(nmat,1)] = rb * v;
  s[momentumIdx(nmat,2)] = rb * w;

  if (bgpreic.empty() || bgtempic.empty())
    Throw("User must specify background pressure and temperature in IC.");

  return s;
}