Added DCDD from previous gradient (#1293)

Description
Previously, the tracer physics used the current time tracer gradient to compute the DCDD stabilization. Using the previous time gradient instead allows to reduce the number of non-linear iterations while still converging to proper solutions.

Testing
All tests still pass.

Documentation
Not applicable.

Former-commit-id: c2c8cfd

Author: oguevremont

CHANGELOG.md

@@ -11,6 +11,12 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/).
 
 ### Changed
 
+- MINOR The tracer gradient used for DCDD stabilization calculations is now from the previous solution, which reduces the number of non-linear iterations. [#1293](https://github.com/chaos-polymtl/lethe/pull/1293)
+
+## [Master] - 2024-09-25
+
+### Changed
+
 - MINOR Remove the GIFs from the repository to reduce its size. [#1289](https://github.com/chaos-polymtl/lethe/pull/1289)
 
 ## [Master] - 2024-09-24

applications_tests/lethe-fluid-sharp/tracer_around_sphere.output

@@ -31,4 +31,4 @@ Tracer statistics :
 	     Min : -0.0071
 	     Max : 1.007
 	 Average : 0.441
-	 Std-Dev : 0.3645
+	 Std-Dev : 0.3645

applications_tests/lethe-fluid/gls_tracer_advection_tanh.output

@@ -8,80 +8,80 @@ Running on 1 MPI rank(s)...
 Transient iteration: 1        Time: 0.02     Time step: 0.02     CFL: 0.567185
 *******************************************************************************
 L2 error velocity : 5.55112e-17
-L2 error tracer : 0.0101988
+L2 error tracer : 0.0118489
 
 *******************************************************************************
 Transient iteration: 2        Time: 0.04     Time step: 0.02     CFL: 0.567185
 *******************************************************************************
 L2 error velocity : 5.55112e-17
-L2 error tracer : 0.0186561
+L2 error tracer : 0.0203913
 
 *******************************************************************************
 Transient iteration: 3        Time: 0.06     Time step: 0.02     CFL: 0.567185
 *******************************************************************************
 L2 error velocity : 5.55112e-17
-L2 error tracer : 0.0253479
+L2 error tracer : 0.0265709
 
 *******************************************************************************
 Transient iteration: 4        Time: 0.08     Time step: 0.02     CFL: 0.567185
 *******************************************************************************
 L2 error velocity : 5.55112e-17
-L2 error tracer : 0.0320367
+L2 error tracer : 0.0329647
 
 *******************************************************************************
 Transient iteration: 5        Time: 0.1      Time step: 0.02     CFL: 0.567185
 *******************************************************************************
 L2 error velocity : 5.55112e-17
-L2 error tracer : 0.0369237
+L2 error tracer : 0.0376841
 
 *******************************************************************************
 Transient iteration: 6        Time: 0.12     Time step: 0.02     CFL: 0.567185
 *******************************************************************************
 L2 error velocity : 5.55112e-17
-L2 error tracer : 0.0416248
+L2 error tracer : 0.0422934
 
 *******************************************************************************
 Transient iteration: 7        Time: 0.14     Time step: 0.02     CFL: 0.567185
 *******************************************************************************
 L2 error velocity : 5.55112e-17
-L2 error tracer : 0.046606
+L2 error tracer : 0.0471852
 
 *******************************************************************************
 Transient iteration: 8        Time: 0.16     Time step: 0.02     CFL: 0.567185
 *******************************************************************************
 L2 error velocity : 5.55112e-17
-L2 error tracer : 0.0500243
+L2 error tracer : 0.050532
 
 *******************************************************************************
 Transient iteration: 9        Time: 0.18     Time step: 0.02     CFL: 0.567185
 *******************************************************************************
 L2 error velocity : 5.55112e-17
-L2 error tracer : 0.0540616
+L2 error tracer : 0.0545255
 
 *******************************************************************************
 Transient iteration: 10       Time: 0.2      Time step: 0.02     CFL: 0.567185
 *******************************************************************************
 L2 error velocity : 5.55112e-17
-L2 error tracer : 0.0579947
- time  error_velocity 
-0.0200   5.551115e-17 
-0.0400   5.551115e-17 
-0.0600   5.551115e-17 
-0.0800   5.551115e-17 
-0.1000   5.551115e-17 
-0.1200   5.551115e-17 
-0.1400   5.551115e-17 
-0.1600   5.551115e-17 
-0.1800   5.551115e-17 
-0.2000   5.551115e-17 
-cells error_tracer 
- 1024 1.019878e-02 
- 1024 1.865606e-02 
- 1024 2.534790e-02 
- 1024 3.203667e-02 
- 1024 3.692370e-02 
- 1024 4.162480e-02 
- 1024 4.660604e-02 
- 1024 5.002435e-02 
- 1024 5.406162e-02 
- 1024 5.799467e-02 
+L2 error tracer : 0.0584056
+ time  error_velocity
+0.0200   5.551115e-17
+0.0400   5.551115e-17
+0.0600   5.551115e-17
+0.0800   5.551115e-17
+0.1000   5.551115e-17
+0.1200   5.551115e-17
+0.1400   5.551115e-17
+0.1600   5.551115e-17
+0.1800   5.551115e-17
+0.2000   5.551115e-17
+cells error_tracer
+ 1024 1.184893e-02
+ 1024 2.039130e-02
+ 1024 2.657085e-02
+ 1024 3.296468e-02
+ 1024 3.768408e-02
+ 1024 4.229345e-02
+ 1024 4.718519e-02
+ 1024 5.053197e-02
+ 1024 5.452552e-02
+ 1024 5.840562e-02

applications_tests/lethe-fluid/time_dependent_boundaries_tracer.output

@@ -13,16 +13,16 @@ Tracer statistics :
 Transient iteration: 1        Time: 0.5      Time step: 0.5      CFL: 2.79578 
 *******************************************************************************
 Tracer statistics : 
-	     Min : 0.436485
-	     Max : 0.499339
-	 Average : 0.470525
-	 Std-Dev : 0.0191026
+	     Min : 0.436653
+	     Max : 0.499349
+	 Average : 0.470664
+	 Std-Dev : 0.019052
 
 *******************************************************************************
-Transient iteration: 2        Time: 1        Time step: 0.5      CFL: 3.0594  
+Transient iteration: 2        Time: 1        Time step: 0.5      CFL: 3.0594
 *******************************************************************************
-Tracer statistics : 
-	     Min : 0.0441899
-	     Max : 0.956945
-	 Average : 0.516389
-	 Std-Dev : 0.25533
+Tracer statistics :
+	     Min : 0.0444703
+	     Max : 0.956999
+	 Average : 0.516871
+	 Std-Dev : 0.255247

include/solvers/tracer_scratch_data.h

@@ -185,6 +185,8 @@ class TracerScratchData
                                                   this->tracer_gradients);
     this->fe_values_tracer.get_function_laplacians(current_solution,
                                                    this->tracer_laplacians);
+    this->fe_values_tracer.get_function_gradients(
+      previous_solutions[0], this->previous_tracer_gradients);
 
     // Gather previous tracer values
     for (unsigned int p = 0; p < previous_solutions.size(); ++p)
@@ -286,8 +288,6 @@ class TracerScratchData
   std::vector<double>                  tracer_diffusivity_0;
   std::vector<double>                  tracer_diffusivity_1;
 
-
-
   // FEValues for the Tracer problem
   FEValues<dim> fe_values_tracer;
   unsigned int  n_dofs;
@@ -303,6 +303,7 @@ class TracerScratchData
   std::vector<Tensor<1, dim>>      tracer_gradients;
   std::vector<double>              tracer_laplacians;
   std::vector<std::vector<double>> previous_tracer_values;
+  std::vector<Tensor<1, dim>>      previous_tracer_gradients;
 
   // Solid signed distance function
   std::vector<double> sdf_values;

include/solvers/vof_scratch_data.h

@@ -173,7 +173,6 @@ class VOFScratchData
                                          this->phase_gradients);
     fe_values_vof.get_function_laplacians(current_solution,
                                           this->phase_laplacians);
-
     fe_values_vof.get_function_gradients(previous_solutions[0],
                                          this->previous_phase_gradients);
 

source/solvers/tracer_assemblers.cc

@@ -43,6 +43,14 @@ TracerAssemblerCore<dim>::assemble_matrix(TracerScratchData<dim> &scratch_data,
       // Store JxW in local variable for faster access;
       const double JxW = JxW_vec[q];
 
+      // We use the previous concentration for the shock capture if the
+      // simulation is transient
+      const Tensor<1, dim> tracer_gradient_for_dcdd =
+        is_steady(method) ? tracer_gradient :
+                            scratch_data.previous_tracer_gradients[q];
+      const double tracer_gradient_for_dcdd_norm =
+        tracer_gradient_for_dcdd.norm();
+
       // Shock capturing viscosity term
       const double order = scratch_data.fe_values_tracer.get_fe().degree;
 
@@ -51,7 +59,8 @@ TracerAssemblerCore<dim>::assemble_matrix(TracerScratchData<dim> &scratch_data,
 
       const double   tolerance = 1e-12;
       Tensor<1, dim> s         = velocity / (velocity.norm() + tolerance);
-      Tensor<1, dim> r = tracer_gradient / (tracer_gradient.norm() + tolerance);
+      Tensor<1, dim> r =
+        tracer_gradient_for_dcdd / (tracer_gradient_for_dcdd_norm + tolerance);
 
       const Tensor<2, dim> k_corr      = (r * s) * outer_product(s, s);
       const Tensor<2, dim> rr          = outer_product(r, r);
@@ -60,7 +69,7 @@ TracerAssemblerCore<dim>::assemble_matrix(TracerScratchData<dim> &scratch_data,
 
       const double d_vdcdd = order * (0.5 * h * h) *
                              (velocity.norm() * velocity.norm()) *
-                             pow(tracer_gradient.norm() * h, order - 1);
+                             pow(tracer_gradient_for_dcdd_norm * h, order - 1);
 
 
 
@@ -149,24 +158,34 @@ TracerAssemblerCore<dim>::assemble_rhs(TracerScratchData<dim>    &scratch_data,
   for (unsigned int q = 0; q < n_q_points; ++q)
     {
       // Gather into local variables the relevant fields
-      const double         diffusivity      = diffusivity_vector[q];
-      const Tensor<1, dim> tracer_gradient  = scratch_data.tracer_gradients[q];
+      const double         diffusivity     = diffusivity_vector[q];
+      const Tensor<1, dim> tracer_gradient = scratch_data.tracer_gradients[q];
+      const Tensor<1, dim> previous_tracer_gradient =
+        scratch_data.previous_tracer_gradients[q];
       const double         tracer_laplacian = scratch_data.tracer_laplacians[q];
       const Tensor<1, dim> velocity         = scratch_data.velocity_values[q];
 
       // Store JxW in local variable for faster access;
       const double JxW = JxW_vec[q];
 
+      // We use the previous concentration for the shock capture if the
+      // simulation is transient
+      const Tensor<1, dim> tracer_gradient_for_dcdd =
+        is_steady(method) ? tracer_gradient : previous_tracer_gradient;
+      const double tracer_gradient_for_dcdd_norm =
+        tracer_gradient_for_dcdd.norm();
+
       // Shock capturing viscosity term
       const double order = scratch_data.fe_values_tracer.get_fe().degree;
 
 
       const double vdcdd = (0.5 * h) * (velocity.norm() * velocity.norm()) *
-                           pow(tracer_gradient.norm() * h, order);
+                           pow(tracer_gradient_for_dcdd_norm * h, order);
 
       const double   tolerance = 1e-12;
       Tensor<1, dim> s         = velocity / (velocity.norm() + tolerance);
-      Tensor<1, dim> r = tracer_gradient / (tracer_gradient.norm() + tolerance);
+      Tensor<1, dim> r =
+        tracer_gradient_for_dcdd / (tracer_gradient_for_dcdd_norm + tolerance);
 
       const Tensor<2, dim> k_corr      = (r * s) * outer_product(s, s);
       const Tensor<2, dim> rr          = outer_product(r, r);

source/solvers/tracer_scratch_data.cc

@@ -21,8 +21,10 @@ TracerScratchData<dim>::allocate()
   // Velocity
   this->velocity_values = std::vector<Tensor<1, dim>>(n_q_points);
   // Tracer
-  this->tracer_values        = std::vector<double>(n_q_points);
-  this->tracer_gradients     = std::vector<Tensor<1, dim>>(n_q_points);
+  this->tracer_values    = std::vector<double>(n_q_points);
+  this->tracer_gradients = std::vector<Tensor<1, dim>>(n_q_points);
+  this->previous_tracer_gradients =
+    std::vector<Tensor<1, dim>>(this->n_q_points);
   this->tracer_laplacians    = std::vector<double>(n_q_points);
   this->tracer_diffusivity   = std::vector<double>(n_q_points);
   this->tracer_diffusivity_0 = std::vector<double>(n_q_points);