update after rebase and fix for deal.ii 8.5.

geodynamics · Nov 23, 2017 · aad87cf · aad87cf
1 parent f60b507
commit aad87cf
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 28 deletions.
diff --git a/source/simulator/assemblers/newton_stokes.cc b/source/simulator/assemblers/newton_stokes.cc
@@ -138,8 +138,8 @@ namespace aspect
                           data.local_matrix(i, j) += ((2.0 * eta * (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
                                                                                              + scratch.grads_phi_u[j] * (viscosity_derivative_wrt_strain_rate * scratch.grads_phi_u[i]) * strain_rate)
-                                                      + one_over_eta * pressure_scaling
-                                                      * pressure_scaling
+                                                      + one_over_eta * this->get_pressure_scaling()
+                                                      * this->get_pressure_scaling()
                                                       * (scratch.phi_p[i] * scratch
                                                          .phi_p[j]))
                                                      * JxW;
@@ -154,8 +154,8 @@ namespace aspect
                         {
                           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
+                                                      + one_over_eta * this->get_pressure_scaling()
+                                                      * this->get_pressure_scaling()
                                                       * (scratch.phi_p[i] * scratch
                                                          .phi_p[j]))
                                                      * JxW;
@@ -314,32 +314,14 @@ namespace aspect
                                         :
                                         1;
 
-<<<<<<< 9c070c9492be916a7195c2cbb62f971e435fee95
-                  if(use_spd_factor)
-                  {
-                  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 * alpha * (scratch.grads_phi_u[i] * (viscosity_derivative_wrt_strain_rate * scratch.grads_phi_u[j]) * strain_rate
-                                                                                         + scratch.grads_phi_u[j] * (viscosity_derivative_wrt_strain_rate * scratch.grads_phi_u[i]) * strain_rate)
-                                                    + derivative_scaling_factor * this->get_pressure_scaling() * scratch.grads_phi_u[i] * 2.0 * viscosity_derivative_wrt_pressure * scratch.phi_p[j] * strain_rate )
-                                                  * JxW;
-
-                        Assert(dealii::numbers::is_finite(data.local_matrix(i,j)),
-                               ExcMessage ("Error: Assembly matrix is not finite." +
-                                           Utilities::to_string(data.local_matrix(i,j)) +
-                                           " = " + Utilities::to_string(eta)));
-                      }
-                  }
-=======
                   if (Newton_stabilisation.second == "SPD" || Newton_stabilisation.second == "symmetric")
                     {
                       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 * alpha * (scratch.grads_phi_u[i] * (viscosity_derivative_wrt_strain_rate * scratch.grads_phi_u[j]) * strain_rate
                                                                                              + scratch.grads_phi_u[j] * (viscosity_derivative_wrt_strain_rate * scratch.grads_phi_u[i]) * strain_rate)
-                                                        + derivative_scaling_factor * pressure_scaling * scratch.grads_phi_u[i] * 2.0 * viscosity_derivative_wrt_pressure * scratch.phi_p[j] * strain_rate )
+                                                        + derivative_scaling_factor * this->get_pressure_scaling() * scratch.grads_phi_u[i] * 2.0 * viscosity_derivative_wrt_pressure * scratch.phi_p[j] * strain_rate )
                                                       * JxW;
 
                             Assert(dealii::numbers::is_finite(data.local_matrix(i,j)),
@@ -348,14 +330,13 @@ namespace aspect
                                                " = " + Utilities::to_string(eta)));
                           }
                     }
->>>>>>> implement the failsafe in a more flexible way togheter with giving more control over the way both the preconditioner and the A block are stabelized.
                   else
                     {
                       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 * 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 )
+                                                        + derivative_scaling_factor * this->get_pressure_scaling() * scratch.grads_phi_u[i] * 2.0 * viscosity_derivative_wrt_pressure * scratch.phi_p[j] * strain_rate )
                                                       * JxW;
 
                             Assert(dealii::numbers::is_finite(data.local_matrix(i,j)),

diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc
@@ -364,7 +364,7 @@ namespace aspect
         prm.declare_entry ("Use Newton stabilization A block", "SPD",
                            Patterns::Selection ("SPD|PD|symmetric|none"),
                            "TODO");
-        prm.declare_entry ("Use Newton failsafe", "true",
+        prm.declare_entry ("Use Newton failsafe", "false",
                            Patterns::Bool (),
                            "Switches on SPD stabilization when solver fails.");
 
@@ -1091,8 +1091,9 @@ 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."));
+
+        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");