backtrace_start

include/rf_solver.h

@@ -120,6 +120,7 @@ extern double var_damp[MAX_VARIABLE_TYPES]; /* variable specific damp factors */
 extern int Newt_Jacobian_Reformation_stride; /*Stride for reformation of jacobian for
                                    modified newton scheme               */
 extern int Time_Jacobian_Reformation_stride;
+extern int Newton_Line_Search_Type;
 extern int modified_newton;               /*boolean flag for modified Newton */
 extern int save_old_A;                    /*boolean flag for saving old A matrix
                                     for resolve reasons with AZTEC.   There

include/rf_solver_const.h

@@ -37,6 +37,9 @@
 #include "dpi.h"        /* defn of Dpi */
 #include "exo_struct.h" /* defn of Exo_DB */
 
+
+#define NLS_FULL_STEP 0
+#define NLS_BACKTRACK 1
 /*
  * Kinds of solvers available...
  */

src/dp_vif.c

@@ -706,6 +706,7 @@ void noahs_ark(void) {
   ddd_add_member(n, &custom_tol3, 1, MPI_DOUBLE);
   ddd_add_member(n, &Newt_Jacobian_Reformation_stride, 1, MPI_INT);
   ddd_add_member(n, &Time_Jacobian_Reformation_stride, 1, MPI_INT);
+  ddd_add_member(n, &Newton_Line_Search_Type, 1, MPI_INT);
   ddd_add_member(n, &modified_newton, 1, MPI_INT);
   ddd_add_member(n, &convergence_rate_tolerance, 1, MPI_DOUBLE);
   ddd_add_member(n, &modified_newt_norm_tol, 1, MPI_DOUBLE);

src/globals.c

@@ -239,6 +239,7 @@ double var_damp[MAX_VARIABLE_TYPES]; /* variable specific damp factors */
 int Newt_Jacobian_Reformation_stride; /*Stride for reformation of jacobian for
                                    modified newton scheme               */
 int Time_Jacobian_Reformation_stride;
+int Newton_Line_Search_Type;
 int modified_newton;               /*boolean flag for modified Newton */
 int save_old_A;                    /*boolean flag for saving old A matrix
                                     for resolve reasons with AZTEC.   There

src/mm_input.c

@@ -6075,6 +6075,10 @@ void rd_solver_specs(FILE *ifp, char *input) {
     strcpy(Stratimikos_File[0], "stratimikos.xml");
     ECHO(echo_string, echo_file);
   }
+  // copy Stratimikos_File[0] to Stratimikos_File[*]
+  for (i = 1; i < MAX_NUM_MATRICES; i++) {
+    strcpy(Stratimikos_File[i], Stratimikos_File[0]);
+  }
 
   strcpy(search_string, "Preconditioner");
 
@@ -6475,6 +6479,35 @@ void rd_solver_specs(FILE *ifp, char *input) {
     Time_Jacobian_Reformation_stride = 0;
   }
 
+  char ls_type[MAX_CHAR_IN_INPUT] = "FULL_STEP";;
+  Newton_Line_Search_Type = NLS_FULL_STEP;
+  int lsread = look_for_optional_string(ifp, "Newton line search type", ls_type, MAX_CHAR_IN_INPUT);
+  snprintf(echo_string, MAX_CHAR_ECHO_INPUT, "%s = %s", "Newton line search type", ls_type);
+  if (lsread) {
+    if (strcmp("FULL_STEP", ls_type) == 0) {
+      Newton_Line_Search_Type = NLS_FULL_STEP;
+    } else if (strcmp("BACKTRACK", ls_type) == 0) {
+      Newton_Line_Search_Type = NLS_BACKTRACK;
+    } else {
+      GOMA_EH(GOMA_ERROR, "Unknown Newton line search type: %s", ls_type);
+    }
+  }
+
+  if (fscanf(ifp, "%le %le %le %le %le", &custom_tol1, &damp_factor2, &custom_tol2, &damp_factor3,
+             &custom_tol3) != 5) {
+    damp_factor2 = damp_factor3 = -1.;
+    custom_tol1 = custom_tol2 = custom_tol3 = -1;
+    rewind(ifp); /* Added to make ibm happy when single relaxation value input. dal/1-6-99 */
+  } else {
+    if ((damp_factor1 <= 1. && damp_factor1 >= 0.) && (damp_factor2 <= 1. && damp_factor2 >= 0.) &&
+        (damp_factor3 <= 1. && damp_factor3 >= 0.)) {
+      SPF(endofstring(echo_string), " %.4g %.4g %.4g %.4g %.4g", custom_tol1, damp_factor2,
+          custom_tol2, damp_factor3, custom_tol3);
+    } else {
+      GOMA_EH(GOMA_ERROR, "All damping factors must be in the range 0 <= fact <=1");
+    }
+  }
+
   look_for(ifp, "Newton correction factor", input, '=');
   if (fscanf(ifp, "%le", &damp_factor1) != 1) {
     GOMA_EH(GOMA_ERROR, "error reading Newton correction factor, expected at least one float");
@@ -8086,10 +8119,7 @@ void rd_eq_specs(FILE *ifp, char *input, const int mn) {
       snprintf(echo_string, MAX_CHAR_ECHO_INPUT, "Stratimikos file = %s for matrix %d", input,
                mtrx_index1);
       ECHO(echo_string, echo_file);
-    } else {
-      // Set stratimikos.xml as defualt stratimikos file
-      strcpy(Stratimikos_File[imtrx], "stratimikos.xml");
-    }
+    } 
 
     iread = look_forward_optional_until(ifp, "Normalized Residual Tolerance", "MATRIX", input, '=');
     if (iread == 1) {

src/mm_sol_nonlinear.c

@@ -19,6 +19,7 @@
 #include "bc_contact.h"
 #include "sl_epetra_interface.h"
 #include "sl_util_structs.h"
+#include <az_aztec_defs.h>
 
 #define GOMA_MM_SOL_NONLINEAR_C
 /* Needed to declare POSIX function drand48 */
@@ -2009,17 +2010,135 @@ int solve_nonlinear_problem(struct GomaLinearSolverData *ams,
      *   UPDATE GOMA UNKNOWNS
      *
      *******************************************************************/
-    for (i = 0; i < NumUnknowns[pg->imtrx]; i++) {
-      x[i] -= damp_factor * var_damp[idv[pg->imtrx][i][0]] * delta_x[i];
-    }
-    exchange_dof(cx, dpi, x, pg->imtrx);
-    if (pd->TimeIntegration != STEADY) {
+    if (Newton_Line_Search_Type == NLS_BACKTRACK) {
+      dbl damp = 1.0;
+      dbl reduction_factor = 0.5;
+      dbl min_damp = 1e-6;
+      dbl *w = alloc_dbl_1(numProcUnknowns, 0.0);
+      dbl *R = alloc_dbl_1(numProcUnknowns, 0.0);
+      dbl *x_save = alloc_dbl_1(numProcUnknowns, 0.0);
+      dcopy1(numProcUnknowns, x, x_save);
+      dbl *xdot_save = alloc_dbl_1(numProcUnknowns, 0.0);
+      dcopy1(numProcUnknowns, xdot, xdot_save);
+
+      int save_jacobian = af->Assemble_Jacobian;
+      int save_residual = af->Assemble_Residual;
+      af->Assemble_Jacobian = FALSE;
+      af->Assemble_Residual = TRUE;
+      err = matrix_fill_full(ams, x, R, x_old, x_older, xdot, xdot_old, x_update, &delta_t, &theta,
+                             First_Elem_Side_BC_Array[pg->imtrx], &time_value, exo, dpi,
+                             &num_total_nodes, &h_elem_avg, &U_norm, NULL);
       for (i = 0; i < NumUnknowns[pg->imtrx]; i++) {
-        xdot[i] -=
-            damp_factor * var_damp[idv[pg->imtrx][i][0]] * delta_x[i] * (1.0 + 2 * theta) / delta_t;
+        x[i] -= damp * delta_x[i];
+      }
+      exchange_dof(cx, dpi, x, pg->imtrx);
+      if (pd->TimeIntegration != STEADY) {
+        for (i = 0; i < NumUnknowns[pg->imtrx]; i++) {
+          xdot[i] -= damp * delta_x[i] * (1.0 + 2 * theta) / delta_t;
+        }
+        exchange_dof(cx, dpi, xdot, pg->imtrx);
       }
-      exchange_dof(cx, dpi, xdot, pg->imtrx);
 
+      exchange_dof(cx, dpi, R, pg->imtrx);
+
+      dbl r_check = L2_norm(R, NumUnknowns[pg->imtrx]);
+
+      dbl best_damp = damp;
+      dbl best_norm = r_check;
+      r_check = 0.5 * sqrt(r_check * r_check);
+      dbl slope = 0;
+      if (strcmp(Matrix_Format, "msr") == 0) {
+        AZ_MATRIX *Amat =
+            AZ_matrix_create(ams->data_org[AZ_N_internal] + ams->data_org[AZ_N_border]);
+        AZ_set_MSR(Amat, ams->bindx, ams->val, ams->data_org, 0, NULL, AZ_LOCAL);
+        AZ_MSR_matvec_mult(delta_x, w, Amat, ams->proc_config);
+        for (int i = 0; i < numProcUnknowns; i++) {
+          slope += w[i] * R[i];
+        }
+        MPI_Allreduce(MPI_IN_PLACE, &slope, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+      }
+      if (slope > 0) {
+        slope = -slope;
+      } else if (slope == 0) {
+        slope = -1;
+      }
+      int skip = FALSE;
+      // Skip if converged
+      if (best_norm < Epsilon[pg->imtrx][0]) {
+        skip = TRUE;
+      }
+
+      while (!skip && (damp > min_damp)) {
+        damp *= reduction_factor;
+        init_vec_value(R, 0.0, numProcUnknowns);
+        for (i = 0; i < NumUnknowns[pg->imtrx]; i++) {
+          x[i] = x_save[i] - damp * delta_x[i];
+        }
+        exchange_dof(cx, dpi, x, pg->imtrx);
+        if (pd->TimeIntegration != STEADY) {
+          for (i = 0; i < NumUnknowns[pg->imtrx]; i++) {
+            xdot[i] = xdot_save[i] - damp * delta_x[i] * (1.0 + 2 * theta) / delta_t;
+          }
+          exchange_dof(cx, dpi, xdot, pg->imtrx);
+        }
+        err = matrix_fill_full(ams, x, R, x_old, x_older, xdot, xdot_old, x_update, &delta_t,
+                               &theta, First_Elem_Side_BC_Array[pg->imtrx], &time_value, exo, dpi,
+                               &num_total_nodes, &h_elem_avg, &U_norm, NULL);
+        exchange_dof(cx, dpi, R, pg->imtrx);
+
+        dbl g_check = L2_norm(R, NumUnknowns[pg->imtrx]);
+        if (isnan(g_check)) {
+          break;
+        }
+
+        if (g_check < best_norm) {
+          best_damp = damp;
+          best_norm = g_check;
+        }
+        if (best_norm < Epsilon[pg->imtrx][0]) {
+          break;
+        }
+        g_check = 0.5 * sqrt(g_check * g_check);
+
+        if (g_check <= r_check + 0.5 * slope * damp) {
+          break;
+        }
+      }
+
+      P0PRINTF("Newton Line Search: best damping factor: %f\n", best_damp);
+      for (i = 0; i < NumUnknowns[pg->imtrx]; i++) {
+        x[i] = x_save[i] - best_damp * delta_x[i];
+      }
+      exchange_dof(cx, dpi, x, pg->imtrx);
+      if (pd->TimeIntegration != STEADY) {
+        for (i = 0; i < NumUnknowns[pg->imtrx]; i++) {
+          xdot[i] = xdot_save[i] - (best_damp * delta_x[i] * (1.0 + 2 * theta) / delta_t);
+        }
+        exchange_dof(cx, dpi, xdot, pg->imtrx);
+      }
+
+      af->Assemble_Jacobian = save_jacobian;
+      af->Assemble_Residual = save_residual;
+      free(w);
+      free(R);
+      free(x_save);
+      free(xdot_save);
+
+    } else if (Newton_Line_Search_Type == NLS_FULL_STEP) {
+      for (i = 0; i < NumUnknowns[pg->imtrx]; i++) {
+        x[i] -= damp_factor * var_damp[idv[pg->imtrx][i][0]] * delta_x[i];
+      }
+      exchange_dof(cx, dpi, x, pg->imtrx);
+      if (pd->TimeIntegration != STEADY) {
+        for (i = 0; i < NumUnknowns[pg->imtrx]; i++) {
+          xdot[i] -= damp_factor * var_damp[idv[pg->imtrx][i][0]] * delta_x[i] * (1.0 + 2 * theta) /
+                     delta_t;
+        }
+        exchange_dof(cx, dpi, xdot, pg->imtrx);
+      }
+    }
+
+    if (pd->TimeIntegration != STEADY) {
       /* Now go back and correct all those dofs in solid regions undergoing newmark-beta
        * transient scheme */
       if (tran->solid_inertia) {

src/sl_stratimikos_interface.cpp

@@ -9,6 +9,7 @@
 #include "Teuchos_ENull.hpp"
 #include "Teuchos_FancyOStream.hpp"
 #include "Teuchos_ParameterList.hpp"
+#include "Teuchos_YamlParameterListCoreHelpers.hpp"
 #include "Teuchos_ParameterListExceptions.hpp"
 #include "Teuchos_Ptr.hpp"
 #include "Teuchos_RCP.hpp"
@@ -77,7 +78,9 @@ int stratimikos_solve(struct GomaLinearSolverData *ams,
     // Get parameters from file
     if (!param_set[imtrx]) {
       param_set[imtrx] = true;
-      solverParams_static[imtrx] = Teuchos::getParametersFromXmlFile(stratimikos_file[imtrx]);
+//      solverParams_static[imtrx] = Teuchos::getParametersFromXmlFile(stratimikos_file[imtrx]);
+      printf("stratimikos_file[imtrx] = %s\n", stratimikos_file[imtrx]);
+      solverParams_static[imtrx] = Teuchos::getParametersFromYamlFile(stratimikos_file[imtrx]);
     }
 
     RCP<Teuchos::ParameterList> solverParams = solverParams_static[imtrx];
@@ -91,6 +94,9 @@ int stratimikos_solve(struct GomaLinearSolverData *ams,
 
     linearSolverBuilder.setParameterList(solverParams);
 
+    auto valid_params = linearSolverBuilder.getValidParameters();
+    Teuchos::writeParameterListToYamlFile(*valid_params, "valid_params.yaml");
+
     // set up solver factory using base/params
     RCP<Thyra::LinearOpWithSolveFactoryBase<double>> solverFactory =
         linearSolverBuilder.createLinearSolveStrategy("");