Constraints

examples/crm/crm_frequency.cpp

@@ -1,6 +1,7 @@
 #include "JacobiDavidson.h"
 #include "TACSBuckling.h"
 #include "TACSIsoShellConstitutive.h"
+#include "TACSKSFailure.h"
 #include "TACSMeshLoader.h"
 #include "TACSShellElementDefs.h"
 
@@ -9,6 +10,11 @@ int main(int argc, char **argv) {
   MPI_Init(&argc, &argv);
   MPI_Comm comm = MPI_COMM_WORLD;
 
+  int rank;
+  MPI_Comm_rank(comm, &rank);
+
+  int use_cg = 0;
+  int num_vecs = 4;  // Number of vectors to include in the eigenvalue function
   int use_lanczos = 1;
   int use_tacs_freq = 0;
   for (int i = 0; i < argc; i++) {
@@ -19,11 +25,21 @@ int main(int argc, char **argv) {
       use_lanczos = 0;
       use_tacs_freq = 1;
     }
+    if (strcmp(argv[i], "use_cg") == 0) {
+      use_cg = 1;
+    }
+    if (sscanf(argv[i], "num_vecs=%d", &num_vecs) == 1) {
+      if (rank == 0) {
+        if (num_vecs < 1) {
+          num_vecs = 1;
+        } else if (num_vecs > 20) {
+          num_vecs = 20;
+        }
+        printf("num_vecs = %d\n", num_vecs);
+      }
+    }
   }
 
-  int rank;
-  MPI_Comm_rank(comm, &rank);
-
   // Write name of BDF file to be load to char array
   const char *filename = "CRM_box_2nd.bdf";
 
@@ -105,8 +121,8 @@ int main(int argc, char **argv) {
 
     // Perform the Frequency Analysis
     int max_lanczos = 60;
-    int num_eigvals = 20;
-    double eig_tol = 1e-8;
+    int num_eigvals = 20;  // 20;
+    double eig_tol = 1e-12;
     TacsScalar sigma = 10.0;
     int output_freq = 1;
 
@@ -118,7 +134,9 @@ int main(int argc, char **argv) {
     double t1 = MPI_Wtime();
     freq_analysis->solve(ksm_print);
     t1 = MPI_Wtime() - t1;
-    printf("Lanczos time: %15.5e\n", t1);
+    if (rank == 0) {
+      printf("Lanczos time: %15.5e\n", t1);
+    }
 
     TacsScalar freq;
     for (int k = 0; k < num_eigvals; k++) {
@@ -127,10 +145,94 @@ int main(int argc, char **argv) {
       eigvalue = sqrt(eigvalue);
       freq = TacsRealPart(eigvalue) / (2.0 * 3.14159);
 
-      printf("TACS frequency[%2d]: %15.6f %15.6f\n", k, TacsRealPart(eigvalue),
-             TacsRealPart(freq));
+      if (rank == 0) {
+        printf("TACS frequency[%2d]: %15.6f %15.6f\n", k,
+               TacsRealPart(eigvalue), TacsRealPart(freq));
+      }
+    }
+
+    TACSBVec *x = assembler->createDesignVec();
+    x->incref();
+    assembler->getDesignVars(x);
+
+    // Set the vector
+    TACSBVec *px = assembler->createDesignVec();
+    px->incref();
+
+    TacsScalar *px_array;
+    int size = px->getArray(&px_array);
+    for (int k = 0; k < size; k++) {
+      px_array[k] = 1.0 - 2.0 * (k % 3);
+    }
+
+    // The function values and the derivative
+    TacsScalar f0 = 0.0, f1 = 0.0;
+    TACSBVec *dfdx = assembler->createDesignVec();
+    dfdx->incref();
+
+    TACSFunction *func = new TACSKSFailure(assembler, 100.0);
+    func->incref();
+
+    // Evaluate the function and the derivative for the first
+    TACSBVec **dfdq = new TACSBVec *[num_vecs];
+    TacsScalar *dfdlam = new TacsScalar[num_vecs];
+    for (int i = 0; i < num_vecs; i++) {
+      dfdlam[i] = 0.0;
+      dfdq[i] = assembler->createVec();
+      dfdq[i]->incref();
+
+      TacsScalar error;
+      freq_analysis->extractEigenvector(i, vec, &error);
+
+      TacsScalar fval;
+      assembler->setVariables(vec);
+      assembler->evalFunctions(1, &func, &fval);
+      assembler->addDVSens(1.0, 1, &func, &dfdx);
+      assembler->addSVSens(1.0, 0.0, 0.0, 1, &func, &dfdq[i]);
+      f0 += fval;
+    }
+
+    // Compute the derivative
+    freq_analysis->addEigenSens(num_vecs, dfdlam, dfdq, dfdx, NULL, use_cg);
+
+    // Perturb the design variables
+    double dh = 1e-6;
+    x->axpy(dh, px);
+    assembler->setDesignVars(x);
+
+    // Solve the frequency analysis again
+    freq_analysis->solve(ksm_print);
+
+    for (int i = 0; i < num_vecs; i++) {
+      TacsScalar error;
+      freq_analysis->extractEigenvector(i, vec, &error);
+
+      TacsScalar fval;
+      assembler->setVariables(vec);
+      assembler->evalFunctions(1, &func, &fval);
+      f1 += fval;
+    }
+
+    // Compute the exact solution
+    TacsScalar ans = px->dot(dfdx);
+
+    // Compute the finite difference
+    TacsScalar fd = (f1 - f0) / dh;
+
+    if (rank == 0) {
+      printf("Ans: %25.10e  FD: %25.10e  Rel. Err: %25.10e \n",
+             TacsRealPart(ans), TacsRealPart(fd),
+             TacsRealPart((ans - fd) / fd));
     }
 
+    func->decref();
+    for (int i = 0; i < num_vecs; i++) {
+      dfdq[i]->decref();
+    }
+    x->decref();
+    dfdx->decref();
+    px->decref();
+
     pc->decref();
     vec->decref();
     ksm->decref();
@@ -158,14 +260,18 @@ int main(int argc, char **argv) {
       freq_analysis->solve(ksm_print);
       t1 = MPI_Wtime() - t1;
 
-      printf("Jacobi-Davidson time: %15.5e\n", t1);
+      if (rank == 0) {
+        printf("Jacobi-Davidson time: %15.5e\n", t1);
+      }
       for (int k = 0; k < num_eigvals; k++) {
         TacsScalar eigvalue = freq_analysis->extractEigenvalue(k, NULL);
         eigvalue = sqrt(eigvalue);
         TacsScalar freq = TacsRealPart(eigvalue) / (2.0 * 3.14159);
 
-        printf("TACS frequency[%2d]: %15.6f %15.6f\n", k,
-               TacsRealPart(eigvalue), TacsRealPart(freq));
+        if (rank == 0) {
+          printf("TACS frequency[%2d]: %15.6f %15.6f\n", k,
+                 TacsRealPart(eigvalue), TacsRealPart(freq));
+        }
       }
       freq_analysis->decref();
     } else {
@@ -185,15 +291,19 @@ int main(int argc, char **argv) {
       double t1 = MPI_Wtime();
       jd->solve(ksm_print);
       t1 = MPI_Wtime() - t1;
-      printf("Jacobi-Davidson time: %15.5e\n", t1);
+      if (rank == 0) {
+        printf("Jacobi-Davidson time: %15.5e\n", t1);
+      }
 
       for (int k = 0; k < num_eigvals; k++) {
         TacsScalar eigvalue = jd->extractEigenvalue(k, NULL);
         eigvalue = sqrt(eigvalue);
         TacsScalar freq = TacsRealPart(eigvalue) / (2.0 * 3.14159);
 
-        printf("TACS frequency[%2d]: %15.6f %15.6f\n", k,
-               TacsRealPart(eigvalue), TacsRealPart(freq));
+        if (rank == 0) {
+          printf("TACS frequency[%2d]: %15.6f %15.6f\n", k,
+                 TacsRealPart(eigvalue), TacsRealPart(freq));
+        }
       }
 
       jd->decref();

examples/cylinder/cylinder_verify.cpp

@@ -1105,9 +1105,11 @@ int main(int argc, char *argv[]) {
   if (grid_study_flag) {
     char file_name[128];
     if (orthotropic_flag) {
-      sprintf(file_name, "ortho_grid_study_order=%d_P=%d.dat", order, P);
+      snprintf(file_name, sizeof(file_name),
+               "ortho_grid_study_order=%d_P=%d.dat", order, P);
     } else {
-      sprintf(file_name, "iso_grid_study_order=%d_P=%d.dat", order, P);
+      snprintf(file_name, sizeof(file_name), "iso_grid_study_order=%d_P=%d.dat",
+               order, P);
     }
     meshStudy(file_name, P, transform, stiffness, load, R, L, alpha, beta,
               order, 8 - order);
@@ -1166,9 +1168,11 @@ int main(int argc, char *argv[]) {
 
   char file_name[128];
   if (orthotropic_flag) {
-    sprintf(file_name, "ortho_cylinder_func_order=%d_nx=%d.dat", order, nx);
+    snprintf(file_name, sizeof(file_name),
+             "ortho_cylinder_func_order=%d_nx=%d.dat", order, nx);
   } else {
-    sprintf(file_name, "iso_cylinder_func_order=%d_nx=%d.dat", order, nx);
+    snprintf(file_name, sizeof(file_name),
+             "iso_cylinder_func_order=%d_nx=%d.dat", order, nx);
   }
 
   FILE *fp = fopen(file_name, "w");

examples/four_bar_mechanism/four_bar_mechanism.cpp

@@ -400,7 +400,7 @@ int main(int argc, char *argv[]) {
   // Extra the data to a file
   for (int pt = 0; pt < 3; pt++) {
     char filename[128];
-    sprintf(filename, "mid_beam_%d.dat", pt + 1);
+    snprintf(filename, sizeof(filename), "mid_beam_%d.dat", pt + 1);
     FILE *fp = fopen(filename, "w");
 
     fprintf(fp, "Variables = t, u0, v0, w0, quantity\n");

examples/mg/mg.cpp

@@ -172,6 +172,9 @@ int main(int argc, char *argv[]) {
   TACSAssembler *assembler[max_nlevels];
   TACSCreator *creator[max_nlevels];
 
+  // Flag to indicate the type of solution method
+  int use_gmres = 0;
+
   // Set the dimension of the largest meshes
   int nx = 128;
   int ny = 128;
@@ -197,6 +200,9 @@ int main(int argc, char *argv[]) {
         nlevels = max_nlevels;
       }
     }
+    if (strcmp(argv[k], "use_gmres") == 0) {
+      use_gmres = 1;
+    }
   }
 
   // Create the multigrid object
@@ -315,16 +321,32 @@ int main(int argc, char *argv[]) {
   TACSBVec *ans = assembler[0]->createVec();
   ans->incref();
 
-  // Allocate the GMRES solution method
-  int gmres_iters = 25;
-  int nrestart = 8;
-  int is_flexible = 0;
-  GMRES *ksm = new GMRES(mg->getMat(0), mg, gmres_iters, nrestart, is_flexible);
-  ksm->incref();
+  TACSKsm *ksm = NULL;
+
+  if (use_gmres) {
+    // Allocate the GMRES solution method
+    int gmres_iters = 25;
+    int nrestart = 8;
+    int is_flexible = 1;
+    ksm = new GMRES(mg->getMat(0), mg, gmres_iters, nrestart, is_flexible);
+
+    // Set a monitor to check on solution progress
+    int freq = 1;
+    ksm->setMonitor(new KSMPrintStdout("GMRES", rank, freq));
 
-  // Set a monitor to check on solution progress
-  int freq = 1;
-  ksm->setMonitor(new KSMPrintStdout("GMRES", rank, freq));
+  } else {
+    // Create the PCG object
+    int reset_iter = 100;
+    int max_reset = 100;
+    ksm = new PCG(mg->getMat(0), mg, reset_iter, max_reset);
+
+    // Set a monitor to check on solution progress
+    int freq = 1;
+    ksm->setMonitor(new KSMPrintStdout("PCG", rank, freq));
+  }
+
+  ksm->incref();
+  ksm->setTolerances(1e-12, 1e-30);
 
   // The initial time
   double t0 = MPI_Wtime();
@@ -343,7 +365,7 @@ int main(int argc, char *argv[]) {
   // Factor the preconditioner
   mg->factor();
 
-  // Compute the solution using GMRES
+  // Compute the solution
   ksm->solve(force, ans);
 
   t0 = MPI_Wtime() - t0;

examples/stiffened_panel/skewed_plate.cpp

@@ -303,7 +303,8 @@ int main(int argc, char *argv[]) {
     // Set the local variables
     tacs->setVariables(vec);
     char file_name[256];
-    sprintf(file_name, "results/tacs_buckling_mode%02d.f5", k);
+    snprintf(file_name, sizeof(file_name), "results/tacs_buckling_mode%02d.f5",
+             k);
     f5->writeToFile(file_name);
   }
 
@@ -322,7 +323,7 @@ int main(int argc, char *argv[]) {
     // Set the local variables
     tacs->setVariables(vec);
     char file_name[256];
-    sprintf(file_name, "results/tacs_mode%02d.f5", k);
+    snprintf(file_name, sizeof(file_name), "results/tacs_mode%02d.f5", k);
     f5->writeToFile(file_name);
 
     if (rank == 0) {

examples/stiffened_panel/stiffened_panel.cpp

@@ -439,8 +439,8 @@ int main(int argc, char *argv[]) {
   //     // Set the local variables
   //     tacs->setVariables(vec);
   //     char file_name[256];
-  //     sprintf(file_name, "results/tacs_buckling_mode%02d.f5", k);
-  //     f5->writeToFile(file_name);
+  //     snprintf(file_name, sizeof(file_name),
+  //     "results/tacs_buckling_mode%02d.f5", k); f5->writeToFile(file_name);
   //   }
 
   //   linear_buckling->decref();

src/TACSAssembler.cpp

@@ -5258,11 +5258,6 @@ void TACSAssembler::addMatXptSensInnerProduct(TacsScalar scale,
   getDataPointers(elementData, &elemVars, &elemPsi, &elemPhi, NULL, &elemXpts,
                   &xptSens, NULL, NULL);
 
-  // Get the design variables from the elements on this process
-  const int maxDVs = maxElementDesignVars;
-  TacsScalar *fdvSens = elementSensData;
-  int *dvNums = elementSensIData;
-
   // Go through each element in the domain and compute the derivative
   // of the residuals with respect to each design variable and multiply by
   // the adjoint variables