diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
index 388d69d93b906e1af1d220d0761bc35e245347d6..6886303a067ec6116c31a67d5e6647e7f5c5d712 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 39a20ecf653cbb6877b664ddff1ee34c2ddb623b..b75a21fecd3f7a23d42cc2c9199650ecc6923bb0 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 e06ff0db8b1c4ac87d1de67aaf44c293d57b9d39..ef4ca9047bc305a057f8c88197730f8fa4dbba07 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 ff1a5aac2052550c6cfd6c6cf592925b79c09ffc..550d6038f16c1020541b6653236d6be799b4c993 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 6ce33e700734f74cd3c9d19f1f5131bc3d9962c1..e756ae1d92e10f44670deefa148e865bb6d9924b 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 83084e12bea1e5c62081635f8c1b6d8317966438..372b6101fc273b288098e7c533a52589daf657a0 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 27e20af297e6476cb3d73ea6d0237f84d03d32d7..620ae625af48a2b2bd3838c459aa890cd8689b00 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 7ee405f8d4ac0456d12ab1b3e3dd1e42e47bc539..42d6b50a5098df4a3d31f6bf4cb36eade3180500 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 aa5548d1a6c519be479558af526dd1f677b07da7..fce35e90c7b4959298206308ff3627dfe813359c 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);