Skip to content
Snippets Groups Projects
Commit 12519e6a authored by wenqing's avatar wenqing
Browse files

[HT] Modified the velocity calculation in order that it can be used for both... 

[HT] Modified the velocity calculation in order that it can be used for both the monolithic and staggered schemes.
parent cba5e249
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 152 additions and 82 deletions
ProcessLib/HT/HTFEM.h
+ 63
− 67
View file @ 12519e6a
......@@ -49,7 +49,8 @@ public:
unsigned const integration_order,
HTMaterialProperties const& material_properties,
const unsigned dof_per_node)
: _element(element),
: HTLocalAssemblerInterface(),
_element(element),
_material_properties(material_properties),
_integration_method(integration_order)
{
......@@ -86,72 +87,6 @@ public:
return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
}
std::vector<double> const& getIntPtDarcyVelocity(
const double t,
GlobalVector const& current_solution,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::vector<double>& cache) const override
{
auto const n_integration_points =
_integration_method.getNumberOfPoints();
auto const indices = NumLib::getIndices(_element.getID(), dof_table);
assert(!indices.empty());
auto const local_x = current_solution.get(indices);
cache.clear();
auto cache_mat = MathLib::createZeroedMatrix<
Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
cache, GlobalDim, n_integration_points);
SpatialPosition pos;
pos.setElementID(_element.getID());
MaterialLib::Fluid::FluidProperty::ArrayType vars;
auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
&local_x[ShapeFunction::NPOINTS], ShapeFunction::NPOINTS);
for (unsigned ip = 0; ip < n_integration_points; ++ip)
{
auto const& ip_data = _ip_data[ip];
auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
pos.setIntegrationPoint(ip);
double T_int_pt = 0.0;
double p_int_pt = 0.0;
NumLib::shapeFunctionInterpolate(local_x, N, T_int_pt, p_int_pt);
vars[static_cast<int>(
MaterialLib::Fluid::PropertyVariableType::T)] = T_int_pt;
vars[static_cast<int>(
MaterialLib::Fluid::PropertyVariableType::p)] = p_int_pt;
auto const K =
_material_properties.porous_media_properties.getIntrinsicPermeability(
t, pos).getValue(t, pos, 0.0, T_int_pt);
auto const mu = _material_properties.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
GlobalDimMatrixType const K_over_mu = K / mu;
cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
if (_material_properties.has_gravity)
{
auto const rho_w =
_material_properties.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Density, vars);
auto const b = _material_properties.specific_body_force;
// here it is assumed that the vector b is directed 'downwards'
cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
}
}
return cache;
}
protected:
MeshLib::Element const& _element;
HTMaterialProperties const& _material_properties;
......@@ -212,6 +147,67 @@ protected:
return thermal_conductivity * I + thermal_dispersivity;
}
std::vector<double> const& getIntPtDarcyVelocityLocal(
const double t, std::vector<double> const& local_p,
std::vector<double> const& local_T, std::vector<double>& cache) const
{
auto const n_integration_points =
_integration_method.getNumberOfPoints();
cache.clear();
auto cache_mat = MathLib::createZeroedMatrix<
Eigen::Matrix<double, GlobalDim, Eigen::Dynamic, Eigen::RowMajor>>(
cache, GlobalDim, n_integration_points);
SpatialPosition pos;
pos.setElementID(_element.getID());
MaterialLib::Fluid::FluidProperty::ArrayType vars;
auto const p_nodal_values = Eigen::Map<const NodalVectorType>(
&local_p[0], ShapeFunction::NPOINTS);
for (unsigned ip = 0; ip < n_integration_points; ++ip)
{
auto const& ip_data = _ip_data[ip];
auto const& N = ip_data.N;
auto const& dNdx = ip_data.dNdx;
pos.setIntegrationPoint(ip);
double T_int_pt = 0.0;
double p_int_pt = 0.0;
NumLib::shapeFunctionInterpolate(local_p, N, p_int_pt);
NumLib::shapeFunctionInterpolate(local_T, N, T_int_pt);
vars[static_cast<int>(
MaterialLib::Fluid::PropertyVariableType::T)] = T_int_pt;
vars[static_cast<int>(
MaterialLib::Fluid::PropertyVariableType::p)] = p_int_pt;
auto const K = _material_properties.porous_media_properties
.getIntrinsicPermeability(t, pos)
.getValue(t, pos, 0.0, T_int_pt);
auto const mu = _material_properties.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Viscosity, vars);
GlobalDimMatrixType const K_over_mu = K / mu;
cache_mat.col(ip).noalias() = -K_over_mu * dNdx * p_nodal_values;
if (_material_properties.has_gravity)
{
auto const rho_w =
_material_properties.fluid_properties->getValue(
MaterialLib::Fluid::FluidPropertyType::Density, vars);
auto const b = _material_properties.specific_body_force;
// here it is assumed that the vector b is directed 'downwards'
cache_mat.col(ip).noalias() += K_over_mu * rho_w * b;
}
}
return cache;
}
};
} // namespace HT
......
ProcessLib/HT/HTLocalAssemblerInterface.h
+ 21
− 0
View file @ 12519e6a
......@@ -17,6 +17,8 @@
namespace ProcessLib
{
struct CoupledSolutionsForStaggeredScheme;
namespace HT
{
template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
......@@ -42,11 +44,30 @@ class HTLocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
public NumLib::ExtrapolatableElement
{
public:
HTLocalAssemblerInterface() : _coupled_solutions(nullptr) {}
void setStaggeredCoupledSolutions(
std::size_t const /*mesh_item_id*/,
CoupledSolutionsForStaggeredScheme* const coupling_term)
{
_coupled_solutions = coupling_term;
}
virtual std::vector<double> const& getIntPtDarcyVelocity(
const double /*t*/,
GlobalVector const& /*current_solution*/,
NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
std::vector<double>& /*cache*/) const = 0;
protected:
// TODO: remove _coupled_solutions or move integration point data from local
// assembler class to a new class to make local assembler unique for each
// process.
/** Pointer to CoupledSolutionsForStaggeredScheme that is set in a member of
* Process class, setCoupledTermForTheStaggeredSchemeToLocalAssemblers.
* It is used for calculate the secondary variables like velocity for
* coupled processes.
*/
CoupledSolutionsForStaggeredScheme* _coupled_solutions;
};
} // namespace HT
......
ProcessLib/HT/HTProcess.cpp
+ 13
− 5
View file @ 12519e6a
......@@ -104,9 +104,9 @@ void HTProcess::assembleWithJacobianConcreteProcess(
}
void HTProcess::preTimestepConcreteProcess(GlobalVector const& x,
const double /*t*/,
const double /*delta_t*/,
const int process_id)
const double /*t*/,
const double /*delta_t*/,
const int process_id)
{
assert(process_id < 2);
......@@ -123,11 +123,19 @@ void HTProcess::preTimestepConcreteProcess(GlobalVector const& x,
auto& x0 = *_xs_previous_timestep[process_id];
MathLib::LinAlg::copy(x, x0);
}
auto& x0 = *_xs_previous_timestep[process_id];
MathLib::LinAlg::setLocalAccessibleVector(x0);
}
void HTProcess::setCoupledTermForTheStaggeredSchemeToLocalAssemblers()
{
DBUG("Set the coupled term for the staggered scheme to local assembers.");
GlobalExecutor::executeMemberOnDereferenced(
&HTLocalAssemblerInterface::setStaggeredCoupledSolutions,
_local_assemblers, _coupled_solutions);
}
} // namespace HT
} // namespace ProcessLib
ProcessLib/HT/HTProcess.h
+ 5
− 3
View file @ 12519e6a
......@@ -67,6 +67,8 @@ public:
return _xs_previous_timestep[process_id].get();
}
void setCoupledTermForTheStaggeredSchemeToLocalAssemblers() override;
private:
void initializeConcreteProcess(
NumLib::LocalToGlobalIndexMap const& dof_table,
......@@ -82,9 +84,9 @@ private:
const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
void preTimestepConcreteProcess(
GlobalVector const& x, double const t, double const dt,
const int process_id) override;
void preTimestepConcreteProcess(GlobalVector const& x, double const t,
double const dt,
const int process_id) override;
const std::unique_ptr<HTMaterialProperties> _material_properties;
......
ProcessLib/HT/MonolithicHTFEM.h
+ 20
− 0
View file @ 12519e6a
......@@ -183,6 +183,26 @@ public:
* in buoyancy effects */
}
}
std::vector<double> const& getIntPtDarcyVelocity(
const double t,
GlobalVector const& current_solution,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::vector<double>& cache) const override
{
auto const indices =
NumLib::getIndices(this->_element.getID(), dof_table);
assert(!indices.empty());
auto local_x = current_solution.get(indices);
std::vector<double> local_T(
std::make_move_iterator(local_x.begin() + local_x.size() / 2),
std::make_move_iterator(local_x.end()));
// only p is kept in local_x
local_x.erase(local_x.begin() + local_x.size() / 2, local_x.end());
return this->getIntPtDarcyVelocityLocal(t, local_x, local_T, cache);
}
};
} // namespace HT
......
ProcessLib/HT/StaggeredHTFEM-impl.h
+ 20
− 4
View file @ 12519e6a
......@@ -30,13 +30,13 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
{
if (coupled_term.variable_id == 0)
{
assembleHydraulicEquation(t, local_M_data, local_K_data, local_b_data,
coupled_term);
assembleHeatTransportEquation(t, local_M_data, local_K_data,
local_b_data, coupled_term);
return;
}
assembleHeatTransportEquation(t, local_M_data, local_K_data, local_b_data,
coupled_term);
assembleHydraulicEquation(t, local_M_data, local_K_data, local_b_data,
coupled_term);
}
template <typename ShapeFunction, typename IntegrationMethod,
......@@ -267,5 +267,21 @@ void StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
}
}
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
std::vector<double> const&
StaggeredHTFEM<ShapeFunction, IntegrationMethod, GlobalDim>::
getIntPtDarcyVelocity(const double t,
GlobalVector const& /*current_solution*/,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::vector<double>& cache) const
{
auto const indices = NumLib::getIndices(this->_element.getID(), dof_table);
assert(!indices.empty());
auto const local_xs =
getCurrentLocalSolutions(*(this->_coupled_solutions), indices);
return this->getIntPtDarcyVelocityLocal(t, local_xs[0], local_xs[1], cache);
}
} // namespace HT
} // namespace ProcessLib
ProcessLib/HT/StaggeredHTFEM.h
+ 6
− 0
View file @ 12519e6a
......@@ -67,6 +67,12 @@ public:
std::vector<double>& local_K_data, std::vector<double>& local_b_data,
LocalCoupledSolutions const& coupled_term) override;
std::vector<double> const& getIntPtDarcyVelocity(
const double t,
GlobalVector const& current_solution,
NumLib::LocalToGlobalIndexMap const& dof_table,
std::vector<double>& cache) const override;
private:
void assembleHydraulicEquation(double const t,
std::vector<double>& local_M_data,
......
ProcessLib/UncoupledProcessesTimeLoop.cpp
+ 4
− 3
View file @ 12519e6a
......@@ -987,11 +987,12 @@ void UncoupledProcessesTimeLoop::outputSolutions(
spd->timestepper->getTimeStep().dt(), pcs_idx);
if (is_staggered_coupling)
{
CoupledSolutionsForStaggeredScheme coupled_xs(
CoupledSolutionsForStaggeredScheme coupled_solutions(
_solutions_of_coupled_processes, 0.0, pcs_idx);
spd->process.setCoupledSolutionsForStaggeredScheme(&coupled_xs);
spd->process.setStaggeredCouplingTermToLocalAssemblers();
spd->process.setCoupledSolutionsForStaggeredScheme(
&coupled_solutions);
spd->process.setCoupledTermForTheStaggeredSchemeToLocalAssemblers();
(output_object.*output_class_member)(pcs, spd->process_output,
timestep, t, x);
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment