fix the Newton failsafe for deal.ii versions before 9.0.0, and add th…

…e aplha to the unsymetrisized version of the equation.
geodynamics · Nov 17, 2017 · 52348d3 · 52348d3
1 parent 0c1fe4b
commit 52348d3
Show file tree

Hide file tree

Showing 2 changed files with 4 additions and 2 deletions.
diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc
@@ -143,7 +143,7 @@ namespace aspect
                       if (scratch.dof_component_indices[i] ==
                           scratch.dof_component_indices[j])
                         {
-                          data.local_matrix(i, j) += ((2.0 * eta * (scratch.grads_phi_u[i] * scratch.grads_phi_u[j]))
+                          data.local_matrix(i, j) += ((2.0 * eta * alpha * (scratch.grads_phi_u[i] * scratch.grads_phi_u[j]))
                                                       + derivative_scaling_factor * alpha * (scratch.grads_phi_u[i] * (viscosity_derivative_wrt_strain_rate * scratch.grads_phi_u[j]) * strain_rate)
                                                       + one_over_eta * pressure_scaling
                                                       * pressure_scaling
@@ -325,7 +325,7 @@ namespace aspect
                       for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
                         for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
                           {
-                            data.local_matrix(i,j) += ( derivative_scaling_factor * (scratch.grads_phi_u[i] * (viscosity_derivative_wrt_strain_rate * scratch.grads_phi_u[j]) * strain_rate)
+                            data.local_matrix(i,j) += ( derivative_scaling_factor * alpha * (scratch.grads_phi_u[i] * (viscosity_derivative_wrt_strain_rate * scratch.grads_phi_u[j]) * strain_rate)
                                                         + derivative_scaling_factor * pressure_scaling * scratch.grads_phi_u[i] * 2.0 * viscosity_derivative_wrt_pressure * scratch.phi_p[j] * strain_rate )
                                                       * JxW;
 

diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc
@@ -1091,6 +1091,8 @@ namespace aspect
         use_Newton_stabilisation_preconditioner = prm.get("Use Newton stabilization preconditioner");
         use_Newton_stabilisation_A_block = prm.get("Use Newton stabilization A block");
         use_Newton_failsafe = prm.get_bool("Use Newton failsafe");
+        AssertThrow(DEAL_II_VERSION_GTE(9,0,0) || !use_Newton_failsafe, ExcMessage("The failsafe option can't be used with a deal.ii "
+        		                                                                   " les then 9.0.0."));
       }
       prm.leave_subsection ();
       prm.enter_subsection ("AMG parameters");