TESLocalAssemblerInner-impl.h
/**
* \file
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
* The code of this file is used to decouple the evaluation of matrix elements
* from the rest of OGS6,
* not all of OGS6 has to be recompiled every time a small change is done.
*/
#pragma once
#include "NumLib/Function/Interpolation.h"
#include "TESLocalAssemblerInner-fwd.h"
#include "TESOGS5MaterialModels.h"
#include "TESReactionAdaptor.h"
namespace ProcessLib
{
namespace TES
{
template <typename Traits>
TESLocalAssemblerInner<Traits>::TESLocalAssemblerInner(
const AssemblyParams& ap, const unsigned num_int_pts,
const unsigned dimension)
: _d{ap, num_int_pts, dimension}
{
}
template <typename Traits>
Eigen::Matrix3d TESLocalAssemblerInner<Traits>::getMassCoeffMatrix(
const unsigned int_pt)
{
// TODO: Dalton's law property
const double dxn_dxm = Adsorption::AdsorptionReaction::dMolarFraction(
_d.vapour_mass_fraction, _d.ap.M_react, _d.ap.M_inert);
const double R = MaterialLib::PhysicalConstant::IdealGasConstant;
const double M_pp = _d.ap.poro / _d.p * _d.rho_GR;
const double M_pT = -_d.ap.poro / _d.T * _d.rho_GR;
const double M_px = (_d.ap.M_react - _d.ap.M_inert) * _d.p / (R * _d.T) *
dxn_dxm * _d.ap.poro;
const double M_Tp = -_d.ap.poro;
const double M_TT =
_d.ap.poro * _d.rho_GR * _d.ap.cpG // TODO: vapour heat capacity
+
(1.0 - _d.ap.poro) * _d.solid_density[int_pt] *
_d.ap.cpS; // TODO: adsorbate heat capacity
const double M_Tx = 0.0;
const double M_xp = 0.0;
const double M_xT = 0.0;
const double M_xx = _d.ap.poro * _d.rho_GR;
Eigen::Matrix3d M;
M << M_pp, M_pT, M_px, M_Tp, M_TT, M_Tx, M_xp, M_xT, M_xx;
return M;
}
template <typename Traits>
typename Traits::LaplaceMatrix
TESLocalAssemblerInner<Traits>::getLaplaceCoeffMatrix(const unsigned /*int_pt*/,
const unsigned dim)
{
const double eta_GR = fluid_viscosity(_d.p, _d.T, _d.vapour_mass_fraction);
const double lambda_F =
fluid_heat_conductivity(_d.p, _d.T, _d.vapour_mass_fraction);
const double lambda_S = _d.ap.solid_heat_cond;
using Mat = typename Traits::MatrixDimDim;
typename Traits::LaplaceMatrix L =
Traits::LaplaceMatrix::Zero(dim * NODAL_DOF, dim * NODAL_DOF);
// TODO: k_rel
// L_pp
Traits::blockDimDim(L, 0, 0, dim, dim) =
Traits::blockDimDim(_d.ap.solid_perm_tensor, 0, 0, dim, dim) *
_d.rho_GR / eta_GR;
// TODO: add zeolite part
// L_TT
Traits::blockDimDim(L, dim, dim, dim, dim) =
Mat::Identity(dim, dim) *
(_d.ap.poro * lambda_F + (1.0 - _d.ap.poro) * lambda_S);
// L_xx
Traits::blockDimDim(L, 2 * dim, 2 * dim, dim, dim) =
Mat::Identity(dim, dim) * (_d.ap.tortuosity * _d.ap.poro * _d.rho_GR *
_d.ap.diffusion_coefficient_component);
return L;
}
template <typename Traits>
Eigen::Matrix3d TESLocalAssemblerInner<Traits>::getAdvectionCoeffMatrix(
const unsigned /*int_pt*/)
{
const double A_pp = 0.0;
const double A_pT = 0.0;
const double A_px = 0.0;
const double A_Tp = 0.0;
const double A_TT = _d.rho_GR * _d.ap.cpG; // porosity?
const double A_Tx = 0.0;
const double A_xp = 0.0;
const double A_xT = 0.0;
const double A_xx = _d.rho_GR; // porosity?
Eigen::Matrix3d A;
A << A_pp, A_pT, A_px, A_Tp, A_TT, A_Tx, A_xp, A_xT, A_xx;
return A;
}
template <typename Traits>
Eigen::Matrix3d TESLocalAssemblerInner<Traits>::getContentCoeffMatrix(
const unsigned /*int_pt*/)
{
const double C_pp = 0.0;
const double C_pT = 0.0;
const double C_px = 0.0;
const double C_Tp = 0.0;
const double C_TT = 0.0;
const double C_Tx = 0.0;
const double C_xp = 0.0;
const double C_xT = 0.0;
const double C_xx = (_d.ap.poro - 1.0) * _d.qR;
Eigen::Matrix3d C;
C << C_pp, C_pT, C_px, C_Tp, C_TT, C_Tx, C_xp, C_xT, C_xx;
return C;
}
template <typename Traits>
Eigen::Vector3d TESLocalAssemblerInner<Traits>::getRHSCoeffVector(
const unsigned int_pt)
{
const double reaction_enthalpy =
_d.ap.react_sys->getEnthalpy(_d.p_V, _d.T, _d.ap.M_react);
const double rhs_p =
(_d.ap.poro - 1.0) * _d.qR; // TODO [CL] body force term
const double rhs_T =
_d.rho_GR * _d.ap.poro * _d.ap.fluid_specific_heat_source +
(1.0 - _d.ap.poro) * _d.qR * reaction_enthalpy +
_d.solid_density[int_pt] * (1.0 - _d.ap.poro) *
_d.ap.solid_specific_heat_source;
// TODO [CL] momentum production term
const double rhs_x =
(_d.ap.poro - 1.0) * _d.qR; // TODO [CL] what if x < 0.0
Eigen::Vector3d rhs;
rhs << rhs_p, rhs_T, rhs_x;
return rhs;
}
template <typename Traits>
void TESLocalAssemblerInner<Traits>::initReaction(const unsigned int_pt)
{
auto const& rate = _d.reaction_adaptor->initReaction(int_pt);
_d.qR = rate.reaction_rate;
_d.reaction_rate[int_pt] = rate.reaction_rate;
_d.solid_density[int_pt] = rate.solid_density;
}
template <typename Traits>
void TESLocalAssemblerInner<Traits>::preEachAssembleIntegrationPoint(
const unsigned int_pt,
const std::vector<double>& localX,
typename Traits::ShapeMatrices const& sm)
{
#ifndef NDEBUG
// fill local data with garbage to aid in debugging
_d.p = _d.T = _d.vapour_mass_fraction = _d.p_V = _d.rho_GR = _d.qR =
std::numeric_limits<double>::quiet_NaN();
#endif
NumLib::shapeFunctionInterpolate(localX, sm.N, _d.p, _d.T,
_d.vapour_mass_fraction);
// pre-compute certain properties
_d.p_V = _d.p * Adsorption::AdsorptionReaction::getMolarFraction(
_d.vapour_mass_fraction, _d.ap.M_react, _d.ap.M_inert);
initReaction(int_pt);
assert(_d.p > 0.0);
assert(_d.T > 0.0);
assert(0.0 <= _d.vapour_mass_fraction && _d.vapour_mass_fraction <= 1.0);
_d.rho_GR = fluid_density(_d.p, _d.T, _d.vapour_mass_fraction);
}
template <typename Traits>
void TESLocalAssemblerInner<Traits>::assembleIntegrationPoint(
unsigned integration_point,
std::vector<double> const& localX,
typename Traits::ShapeMatrices const& sm,
const double weight,
Eigen::Map<typename Traits::LocalMatrix>& local_M,
Eigen::Map<typename Traits::LocalMatrix>& local_K,
Eigen::Map<typename Traits::LocalVector>& local_b)
{
preEachAssembleIntegrationPoint(integration_point, localX, sm);
auto const N =
static_cast<unsigned>(sm.dNdx.cols()); // number of integration points
auto const D =
static_cast<unsigned>(sm.dNdx.rows()); // global dimension: 1, 2 or 3
assert(N * NODAL_DOF == local_M.cols());
auto const laplaceCoeffMat = getLaplaceCoeffMatrix(integration_point, D);
assert(laplaceCoeffMat.cols() == D * NODAL_DOF);
auto const massCoeffMat = getMassCoeffMatrix(integration_point);
auto const advCoeffMat = getAdvectionCoeffMatrix(integration_point);
auto const contentCoeffMat = getContentCoeffMatrix(integration_point);
// calculate velocity
assert((unsigned)sm.dNdx.rows() == _d.velocity.size() &&
(unsigned)sm.dNdx.cols() == _d.velocity[0].size());
auto const velocity =
(Traits::blockDimDim(laplaceCoeffMat, 0, 0, D, D) *
(sm.dNdx * Eigen::Map<const typename Traits::Vector1Comp>(localX.data(),
N) // grad_p
/
-_d.rho_GR))
.eval();
assert(velocity.size() == D);
for (unsigned d = 0; d < D; ++d)
{
_d.velocity[d][integration_point] = velocity[d];
}
auto const detJ_w_im_NT =
(sm.detJ * weight * sm.integralMeasure * sm.N.transpose()).eval();
auto const detJ_w_im_NT_N = (detJ_w_im_NT * sm.N).eval();
assert(detJ_w_im_NT_N.rows() == N && detJ_w_im_NT_N.cols() == N);
auto const detJ_w_im_NT_vT_dNdx =
(detJ_w_im_NT * velocity.transpose() * sm.dNdx).eval();
assert(detJ_w_im_NT_vT_dNdx.rows() == N && detJ_w_im_NT_vT_dNdx.cols() == N);
for (unsigned r = 0; r < NODAL_DOF; ++r)
{
for (unsigned c = 0; c < NODAL_DOF; ++c)
{
Traits::blockShpShp(local_K, N * r, N * c, N, N).noalias() +=
sm.detJ * weight * sm.integralMeasure * sm.dNdx.transpose() *
Traits::blockDimDim(laplaceCoeffMat, D * r, D * c, D, D) *
sm.dNdx // end Laplacian part
+ detJ_w_im_NT_N * contentCoeffMat(r, c) +
detJ_w_im_NT_vT_dNdx * advCoeffMat(r, c);
Traits::blockShpShp(local_M, N * r, N * c, N, N).noalias() +=
detJ_w_im_NT_N * massCoeffMat(r, c);
}
}
auto const rhsCoeffVector = getRHSCoeffVector(integration_point);
for (unsigned r = 0; r < NODAL_DOF; ++r)
{
Traits::blockShp(local_b, N * r, N).noalias() +=
rhsCoeffVector(r) * sm.N.transpose() * sm.detJ * weight *
sm.integralMeasure;
}
}
template <typename Traits>
void TESLocalAssemblerInner<Traits>::preEachAssemble()
{
if (_d.ap.iteration_in_current_timestep == 1)
{
if (_d.ap.number_of_try_of_iteration ==
1) // TODO has to hold if the above holds.
{
_d.solid_density_prev_ts = _d.solid_density;
_d.reaction_rate_prev_ts = _d.reaction_rate;
_d.reaction_adaptor->preZerothTryAssemble();
}
else
{
_d.solid_density = _d.solid_density_prev_ts;
}
}
}
} // namespace TES
} // namespace ProcessLib