peterrum
commented
8 months ago Personally, I would go one step further by merging the steady and transient operator. The diff is minor:
index 195e210..42d1ac2 100644
--- a/a.cc
+++ b/b.cc
@@ -1,18 +1,24 @@
+
+
+template <int dim, typename number>
+NavierStokesTransientOperator<dim, number>::NavierStokesTransientOperator() =
+ default;
+
/**
* The expressions calculated in this cell integral are:
- * (q,∇δu) + (v,(u·∇)δu) + (v,(δu·∇)u) - (∇·v,δp) + ν(∇v,∇δu) (Weak form
- * Jacobian),
+ * (q,∇δu) + (v,∂t δu) + (v,(u·∇)δu) + (v,(δu·∇)u) - (∇·v,δp) + ν(∇v,∇δu) (Weak
+ * form Jacobian),
* plus three additional terms in the case of SUPG-PSPG stabilization:
- * \+ ((u·∇)δu + (δu·∇)u + ∇δp - ν∆δu)τ·∇q (PSPG Jacobian)
- * \+ ((u·∇)δu + (δu·∇)u + ∇δp - ν∆δu)τu·∇v (SUPG Jacobian Part 1)
- * \+ ((u·∇)u + ∇p - ν∆u - f )τδu·∇v (SUPG Jacobian Part 2),
+ * \+ (∂t δu +(u·∇)δu + (δu·∇)u + ∇δp - ν∆δu)τ·∇q (PSPG Jacobian)
+ * \+ (∂t δu +(u·∇)δu + (δu·∇)u + ∇δp - ν∆δu)τu·∇v (SUPG Jacobian Part 1)
+ * \+ (∂t u +(u·∇)u + ∇p - ν∆u - f )τδu·∇v (SUPG Jacobian Part 2),
* plus two additional terms in the case of full gls stabilization:
- * \+ ((u·∇)δu + (δu·∇)u + ∇δp - ν∆δu)τ(−ν∆v) (GLS Jacobian)
+ * \+ (∂t δu +(u·∇)δu + (δu·∇)u + ∇δp - ν∆δu)τ(−ν∆v) (GLS Jacobian)
* \+ (∇·δu)τ'(∇·v) (LSIC Jacobian).
*/
template <int dim, typename number>
void
-NavierStokesSteadyOperator<dim, number>::do_cell_integral_local(
+NavierStokesTransientOperator<dim, number>::do_cell_integral_local(
FECellIntegrator &integrator) const
{
integrator.evaluate(EvaluationFlags::values | EvaluationFlags::gradients |
@@ -22,6 +28,15 @@ NavierStokesSteadyOperator<dim, number>::do_cell_integral_local(
const auto h = integrator.read_cell_data(this->get_element_size());
+ // Vector for the BDF coefficients
+ const Vector<double> &bdf_coefs =
+ this->simulation_control->get_bdf_coefficients();
+ const auto time_steps_vector =
+ this->simulation_control->get_time_steps_vector();
+ const double dt = time_steps_vector[0];
+ const double sdt = 1. / dt;
+
+
for (const auto q : integrator.quadrature_point_indices())
{
// Evaluate source term function
@@ -50,15 +65,20 @@ NavierStokesSteadyOperator<dim, number>::do_cell_integral_local(
auto previous_hessian_diagonal =
this->nonlinear_previous_hessian_diagonal(cell, q);
+ // Time derivatives of previous solutions
+ auto previous_time_derivatives =
+ this->time_derivatives_previous_solutions(cell, q);
+
// Calculate tau
VectorizedArray<number> u_mag_squared = 1e-12;
for (unsigned int k = 0; k < dim; ++k)
u_mag_squared += Utilities::fixed_power<2>(previous_values[k]);
const auto tau =
- 1. / std::sqrt(4. * u_mag_squared / h / h +
- 9. * Utilities::fixed_power<2>(
- 4. * this->kinematic_viscosity / (h * h)));
+ 1. /
+ std::sqrt(Utilities::fixed_power<2>(sdt) + 4. * u_mag_squared / h / h +
+ 9. * Utilities::fixed_power<2>(
+ 4. * this->kinematic_viscosity / (h * h)));
const auto tau_lsic = std::sqrt(u_mag_squared) * h * 0.5;
@@ -78,6 +98,8 @@ NavierStokesSteadyOperator<dim, number>::do_cell_integral_local(
value_result[i] += gradient[i][k] * previous_values[k] +
previous_gradient[i][k] * value[k];
}
+ // +(v,∂t δu)
+ value_result[i] += bdf_coefs[0] * value[i];
}
// PSPG Jacobian
@@ -91,6 +113,8 @@ NavierStokesSteadyOperator<dim, number>::do_cell_integral_local(
gradient[i][k] * previous_values[k] +
previous_gradient[i][k] * value[k]);
}
+ // +(∂t δu)·τ∇q
+ gradient_result[dim][i] += tau * bdf_coefs[0] * value[i];
}
// (∇δp)τ·∇q
gradient_result[dim] += tau * gradient[dim];
@@ -110,24 +134,27 @@ NavierStokesSteadyOperator<dim, number>::do_cell_integral_local(
previous_gradient[i][l] * value[l] -
this->kinematic_viscosity * hessian_diagonal[i][l]);
}
-
- // +(∇δp)τ(u·∇)v
+ // +(∇δp + ∂t δu)τ(u·∇)v
gradient_result[i][k] +=
- tau * previous_values[k] * gradient[dim][i];
+ tau * previous_values[k] *
+ (gradient[dim][i] + bdf_coefs[0] * value[i]);
// Part 2
for (unsigned int l = 0; l < dim; ++l)
{
- // +((u·∇)u -ν∆u)τ(δu·∇)v
+ // +((u·∇)u - ν∆u)τ(δu·∇)v
gradient_result[i][k] +=
tau * value[k] *
(previous_gradient[i][l] * previous_values[l] -
this->kinematic_viscosity *
previous_hessian_diagonal[i][l]);
}
- // +(∇p - f)τ(δu·∇)v
+ // +(∇p - f + ∂t u)τ(δu·∇)v
gradient_result[i][k] +=
- tau * value[k] * (previous_gradient[dim][i] - source_value[i]);
+ tau * value[k] *
+ (previous_gradient[dim][i] - source_value[i] +
+ bdf_coefs[0] * previous_values[i] +
+ previous_time_derivatives[i]);
}
}
@@ -149,9 +176,10 @@ NavierStokesSteadyOperator<dim, number>::do_cell_integral_local(
this->kinematic_viscosity * hessian_diagonal[i][l]);
}
- // +(∇δp)τ(−ν∆v)
+ // +(∇δp + ∂t δu)τ(−ν∆v)
hessian_result[i][k][k] +=
- tau * -this->kinematic_viscosity * gradient[dim][i];
+ tau * -this->kinematic_viscosity *
+ (gradient[dim][i] + bdf_coefs[0] * value[i]);
// LSIC term
// (∇·δu)τ'(∇·v)
@@ -160,7 +188,6 @@ NavierStokesSteadyOperator<dim, number>::do_cell_integral_local(
}
}
-
integrator.submit_gradient(gradient_result, q);
integrator.submit_value(value_result, q);
integrator.submit_hessian(hessian_result, q);
@@ -172,17 +199,17 @@ NavierStokesSteadyOperator<dim, number>::do_cell_integral_local(
/**
* The expressions calculated in this cell integral are:
- * (q, ∇·u) + (v,(u·∇)u) - (∇·v,p) + ν(∇v,∇u) - (v,f) (Weak form),
+ * (q, ∇·u) + (v,∂t u) + (v,(u·∇)u) - (∇·v,p) + ν(∇v,∇u) - (v,f) (Weak form),
* plus two additional terms in the case of SUPG-PSPG stabilization:
- * \+ ((u·∇)u + ∇p - ν∆u - f)τ∇·q (PSPG term)
- * \+ ((u·∇)u + ∇p - ν∆u - f)τu·∇v (SUPG term),
+ * \+ (∂t u +(u·∇)u + ∇p - ν∆u - f)τ∇·q (PSPG term)
+ * \+ (∂t u +(u·∇)u + ∇p - ν∆u - f)τu·∇v (SUPG term),
* plus two additional terms in the case of full gls stabilization:
- * \+ ((u·∇)u + ∇p - ν∆u - f)τ(−ν∆v) (GLS term)
+ * \+ (∂t u +(u·∇)u + ∇p - ν∆u - f)τ(−ν∆v) (GLS term)
* \+ (∇·u)τ'(∇·v) (LSIC term).
*/
template <int dim, typename number>
void
-NavierStokesSteadyOperator<dim, number>::local_evaluate_residual(
+NavierStokesTransientOperator<dim, number>::local_evaluate_residual(
const MatrixFree<dim, number> &matrix_free,
VectorType &dst,
const VectorType &src,
@@ -199,6 +226,14 @@ NavierStokesSteadyOperator<dim, number>::local_evaluate_residual(
const auto h = integrator.read_cell_data(this->get_element_size());
+ // Vector for the BDF coefficients
+ const Vector<double> &bdf_coefs =
+ this->simulation_control->get_bdf_coefficients();
+ const auto time_steps_vector =
+ this->simulation_control->get_time_steps_vector();
+ const double dt = time_steps_vector[0];
+ const double sdt = 1. / dt;
+
for (const auto q : integrator.quadrature_point_indices())
{
// Evaluate source term function
@@ -216,13 +251,18 @@ NavierStokesSteadyOperator<dim, number>::local_evaluate_residual(
typename FECellIntegrator::gradient_type hessian_diagonal =
integrator.get_hessian_diagonal(q);
+ // Time derivatives of previoussolutions
+ auto previous_time_derivatives =
+ this->time_derivatives_previous_solutions(cell, q);
+
// Calculate tau
VectorizedArray<number> u_mag_squared = 1e-12;
for (unsigned int k = 0; k < dim; ++k)
u_mag_squared += Utilities::fixed_power<2>(value[k]);
const auto tau =
- 1. / std::sqrt(4. * u_mag_squared / h / h +
+ 1. / std::sqrt(Utilities::fixed_power<2>(sdt) +
+ 4. * u_mag_squared / h / h +
9. * Utilities::fixed_power<2>(
4. * this->kinematic_viscosity / (h * h)));
@@ -240,8 +280,9 @@ NavierStokesSteadyOperator<dim, number>::local_evaluate_residual(
gradient_result[i] = this->kinematic_viscosity * gradient[i];
// -(∇·v,p)
gradient_result[i][i] += -value[dim];
- // -(v,f)
- value_result[i] = -source_value[i];
+ // +(v,-f + ∂t u)
+ value_result[i] = -source_value[i] + bdf_coefs[0] * value[i] +
+ previous_time_derivatives[i];
// +(q,∇·u)
value_result[dim] += gradient[i][i];
@@ -255,16 +296,17 @@ NavierStokesSteadyOperator<dim, number>::local_evaluate_residual(
// PSPG term
for (unsigned int i = 0; i < dim; ++i)
{
- // (-f)·τ∇q
- gradient_result[dim][i] += -tau * source_value[i];
-
for (unsigned int k = 0; k < dim; ++k)
{
- //(-ν∆u + (u·∇)u)·τ∇q
+ // (-ν∆u + (u·∇)u)·τ∇q
gradient_result[dim][i] +=
tau * (-this->kinematic_viscosity * hessian_diagonal[i][k] +
gradient[i][k] * value[k]);
}
+ // +(-f + ∂t u)·τ∇q
+ gradient_result[dim][i] +=
+ tau * (-source_value[i] + bdf_coefs[0] * value[i] +
+ previous_time_derivatives[i]);
}
// +(∇p)τ∇·q
gradient_result[dim] += tau * gradient[dim];
@@ -276,15 +318,20 @@ NavierStokesSteadyOperator<dim, number>::local_evaluate_residual(
{
for (unsigned int l = 0; l < dim; ++l)
{
- // (-ν∆u + (u·∇)u)τ(u·∇)v
+ // (-ν∆u)τ(u·∇)v
+ gradient_result[i][k] +=
+ -tau * this->kinematic_viscosity * value[k] *
+ hessian_diagonal[i][l];
+
+ // + ((u·∇)u)τ(u·∇)v
gradient_result[i][k] +=
- tau * value[k] *
- (-this->kinematic_viscosity * hessian_diagonal[i][l] +
- gradient[i][l] * value[l]);
+ tau * value[k] * gradient[i][l] * value[l];
}
- // + (∇p - f)τ(u·∇)v
+ // + (∇p - f + ∂t u)τ(u·∇)v
gradient_result[i][k] +=
- tau * value[k] * (gradient[dim][i] - source_value[i]);
+ tau * value[k] *
+ (gradient[dim][i] - source_value[i] +
+ bdf_coefs[0] * value[i] + previous_time_derivatives[i]);
}
}
@@ -305,10 +352,12 @@ NavierStokesSteadyOperator<dim, number>::local_evaluate_residual(
hessian_diagonal[i][l] +
gradient[i][l] * value[l]);
}
- // + (∇p - f)τ(−ν∆v)
+ // + (∇p - f + ∂t u)τ(−ν∆v)
hessian_result[i][k][k] +=
tau * -this->kinematic_viscosity *
- (gradient[dim][i] - source_value[i]);
+ (gradient[dim][i] - source_value[i] +
+ +bdf_coefs[0] * value[i] +
+ previous_time_derivatives[i]);
// LSIC term
// (∇·u)τ'(∇·v)
blaisb
lpsaavedra
Description of the problem
There was a lot of code duplication introduced in #1049 since the GLS operator only adds two additional terms to the SUPG/PSPG operator.
Description of the solution
This PR unifies the SUPG/PSPG and GLS steady operators, and the SUPG/PSPG and GLS transient operators, in only one operator that uses the stabilization type and the is_bdf function to identify what terms to take into account.
How Has This Been Tested?
All the matrix-free application tests pass without any issue.