From b8e22cac39dd2e271e227a84a73b88b70df0aa99 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Wed, 5 Oct 2016 16:13:31 +0200
Subject: [PATCH] [PCS] Added the computation of the secondary variables after
during the un-coupling
---
.../LiquidFlowLocalAssembler-impl.h | 62 ++++++++++++++++++-
.../LiquidFlow/LiquidFlowLocalAssembler.h | 2 +-
ProcessLib/LiquidFlow/LiquidFlowProcess.cpp | 8 +++
ProcessLib/LiquidFlow/LiquidFlowProcess.h | 2 +
ProcessLib/LocalAssemblerInterface.cpp | 10 +--
ProcessLib/LocalAssemblerInterface.h | 7 ++-
ProcessLib/Process.cpp | 5 ++
ProcessLib/Process.h | 5 ++
ProcessLib/UncoupledProcessesTimeLoop.cpp | 2 +
9 files changed, 95 insertions(+), 8 deletions(-)
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 388d69d93b9..6886303a067 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -154,9 +154,67 @@ void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
template <typename ShapeFunction, typename IntegrationMethod,
unsigned GlobalDim>
void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
- computeSecondaryVariable(std::vector<double> const& local_x)
+ computeSecondaryVariableConcrete(std::vector<double> const& local_x)
{
- (void) local_x;
+ auto const local_matrix_size = local_x.size();
+ assert(local_matrix_size == ShapeFunction::NPOINTS);
+
+ const auto local_p_vec =
+ MathLib::toVector<NodalVectorType>(local_x, local_matrix_size);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ SpatialPosition pos;
+ pos.setElementID(_element.getID());
+ _material_properties.setMaterialID(pos);
+ const Eigen::MatrixXd& perm =
+ _material_properties.getPermeability(t, pos, _element.getDimension());
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& sm = _shape_matrices[ip];
+ double p = 0.;
+ NumLib::shapeFunctionInterpolate(local_x, sm.N, p);
+ // TODO : compute _temperature from the heat transport pcs
+
+ const double rho_g =
+ _material_properties.getLiquidDensity(p, _temperature) *
+ _gravitational_acceleration;
+ // Compute viscosity:
+ const double mu = _material_properties.getViscosity(p, _temperature);
+
+ // Assemble Laplacian, K, and RHS by the gravitational term
+ if (perm.size() == 1) // Save time for isotropic permeability.
+ {
+ // Use scalar number for isotropic permeability
+ // to save the computation time.
+ const double K = perm(0, 0) / mu;
+ // Compute the velocity
+ GlobalDimVectorType darcy_velocity = -K * sm.dNdx * local_p_vec;
+ // gravity term
+ if (_gravitational_axis_id >= 0)
+ darcy_velocity[GlobalDim - 1] -= K * rho_g;
+ for (unsigned d = 0; d < GlobalDim; ++d)
+ {
+ _darcy_velocities[d][ip] = darcy_velocity[d];
+ }
+ }
+ else
+ {
+ // Compute the velocity
+ GlobalDimVectorType darcy_velocity = -perm * sm.dNdx * local_p_vec / mu;
+ if (_gravitational_axis_id >= 0)
+ {
+ darcy_velocity.noalias() -= rho_g * perm.col(GlobalDim - 1) / mu;
+ }
+ for (unsigned d = 0; d < GlobalDim; ++d)
+ {
+ _darcy_velocities[d][ip] = darcy_velocity[d];
+ }
+ }
+ }
+
}
} // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index 39a20ecf653..b75a21fecd3 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -85,7 +85,7 @@ public:
std::vector<double>& local_K_data,
std::vector<double>& local_b_data) override;
- void computeSecondaryVariable(std::vector<double> const& local_x) override;
+ void computeSecondaryVariableConcrete(std::vector<double> const& local_x) override;
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
index e06ff0db8b1..ef4ca9047bc 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -113,5 +113,13 @@ void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
dx_dx, M, K, b, Jac);
}
+void LiquidFlowProcess::computeSecondaryVariableConcrete(GlobalVector const& x)
+{
+ DBUG("Compute the velocity for LiquidFlowProcess.");
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LiquidFlowLocalAssemblerInterface::computeSecondaryVariable,
+ _local_assemblers, *_local_to_global_index_map, x);
+}
+
} // end of namespace
} // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
index ff1a5aac205..550d6038f16 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowProcess.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -69,6 +69,8 @@ public:
double const gravitational_acceleration,
BaseLib::ConfigTree const& config);
+ void computeSecondaryVariableConcrete(GlobalVector const& x) override;
+
bool isLinear() const override { return true; }
private:
void initializeConcreteProcess(
diff --git a/ProcessLib/LocalAssemblerInterface.cpp b/ProcessLib/LocalAssemblerInterface.cpp
index 6ce33e70073..e756ae1d92e 100644
--- a/ProcessLib/LocalAssemblerInterface.cpp
+++ b/ProcessLib/LocalAssemblerInterface.cpp
@@ -26,11 +26,13 @@ void LocalAssemblerInterface::assembleWithJacobian(
"assembler.");
}
-void LocalAssemblerInterface::computeSecondaryVariable(std::vector<double> const& /*local_x*/)
+void LocalAssemblerInterface::computeSecondaryVariable(std::size_t const mesh_item_id,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ GlobalVector const& x)
{
- OGS_FATAL(
- "computeSecondaryVariable(...) function is not implemented in the local "
- "assembler.");
+ auto const indices = NumLib::getIndices(mesh_item_id, dof_table);
+ auto const local_x = x.get(indices);
+ computeSecondaryVariableConcrete(local_x);
}
void LocalAssemblerInterface::preTimestep(
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index 83084e12bea..372b6101fc2 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -45,7 +45,9 @@ public:
std::vector<double>& local_b_data,
std::vector<double>& local_Jac_data);
- virtual void computeSecondaryVariable(std::vector<double> const& local_x);
+ virtual void computeSecondaryVariable(std::size_t const mesh_item_id,
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ GlobalVector const& x);
virtual void preTimestep(std::size_t const mesh_item_id,
NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -72,6 +74,9 @@ private:
}
virtual void postTimestepConcrete(std::vector<double> const& /*local_x*/) {}
+
+ virtual void computeSecondaryVariableConcrete
+ (std::vector<double> const& /*local_x*/) {}
};
} // namespace ProcessLib
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 27e20af297e..620ae625af4 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -242,6 +242,11 @@ void Process::postTimestep(GlobalVector const& x)
postTimestepConcreteProcess(x);
}
+void Process::computeSecondaryVariable(GlobalVector const& x)
+{
+ computeSecondaryVariableConcrete(x);
+}
+
void Process::preIteration(const unsigned iter, const GlobalVector &x)
{
// In every new iteration cached values of secondary variables are expired.
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 7ee405f8d4a..42d6b50a509 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -59,6 +59,9 @@ public:
void preIteration(const unsigned iter,
GlobalVector const& x) override final;
+ /// computeSecondaryVariable for the coupled equations or for output.
+ void computeSecondaryVariable(GlobalVector const& x);
+
NumLib::IterationResult postIteration(GlobalVector const& x) override final;
void initialize();
@@ -148,6 +151,8 @@ private:
virtual void preIterationConcreteProcess(const unsigned /*iter*/,
GlobalVector const& /*x*/){}
+ virtual void computeSecondaryVariableConcrete(GlobalVector const& /*x*/) {};
+
virtual NumLib::IterationResult postIterationConcreteProcess(
GlobalVector const& /*x*/)
{
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index aa5548d1a6c..fce35e90c7b 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -505,6 +505,8 @@ bool UncoupledProcessesTimeLoop::loop()
nonlinear_solver_succeeded = solveOneTimeStepOneProcess(
x, timestep, t, delta_t, *spd, *_output);
+ pcs.computeSecondaryVariable(x);
+
INFO("[time] Solving process #%u took %g s in timestep #%u.",
timestep, time_timestep.elapsed(), pcs_idx);
--
GitLab