diff --git a/examples/crm/crm_frequency.cpp b/examples/crm/crm_frequency.cpp index 527890643..be30c7dfa 100644 --- a/examples/crm/crm_frequency.cpp +++ b/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(); diff --git a/examples/cylinder/cylinder_verify.cpp b/examples/cylinder/cylinder_verify.cpp index e3a3d2b8c..2416691d4 100644 --- a/examples/cylinder/cylinder_verify.cpp +++ b/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"); diff --git a/examples/four_bar_mechanism/four_bar_mechanism.cpp b/examples/four_bar_mechanism/four_bar_mechanism.cpp index 647d3b2e0..f81a660b0 100644 --- a/examples/four_bar_mechanism/four_bar_mechanism.cpp +++ b/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"); diff --git a/examples/mg/mg.cpp b/examples/mg/mg.cpp index 71fe079c5..322ab224f 100644 --- a/examples/mg/mg.cpp +++ b/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; diff --git a/examples/stiffened_panel/skewed_plate.cpp b/examples/stiffened_panel/skewed_plate.cpp index 478b8bd16..02641291c 100644 --- a/examples/stiffened_panel/skewed_plate.cpp +++ b/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) { diff --git a/examples/stiffened_panel/stiffened_panel.cpp b/examples/stiffened_panel/stiffened_panel.cpp index b63ebfe82..abd3536b7 100644 --- a/examples/stiffened_panel/stiffened_panel.cpp +++ b/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(); diff --git a/src/TACSAssembler.cpp b/src/TACSAssembler.cpp index 082686485..696785924 100644 --- a/src/TACSAssembler.cpp +++ b/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 diff --git a/src/TACSBuckling.cpp b/src/TACSBuckling.cpp index 12f5ed55b..765ba7a27 100644 --- a/src/TACSBuckling.cpp +++ b/src/TACSBuckling.cpp @@ -20,6 +20,56 @@ #include "tacslapack.h" +class TACSMatCombo : public TACSMat { + public: + TACSMatCombo(TacsScalar _Acoef, TacsScalar _Bcoef, TACSMat *_A, TACSMat *_B, + TACSAssembler *_assembler) { + Acoef = _Acoef; + Bcoef = _Bcoef; + + A = _A; + B = _B; + A->incref(); + B->incref(); + + assembler = _assembler; + assembler->incref(); + + temp = assembler->createVec(); + temp->incref(); + } + ~TACSMatCombo() { + A->decref(); + B->decref(); + assembler->decref(); + temp->decref(); + } + + TACSVec *createVec() { return A->createVec(); } + + void setCoefficients(TacsScalar _Acoef, TacsScalar _Bcoef) { + Acoef = _Acoef; + Bcoef = _Bcoef; + } + + void mult(TACSVec *x, TACSVec *y) { + A->mult(x, y); + if (!(TacsRealPart(Acoef) == 1.0)) { + y->scale(Acoef); + } + + B->mult(x, temp); + y->axpy(Bcoef, temp); + assembler->applyBCs(y); + } + + private: + TacsScalar Acoef, Bcoef; + TACSMat *A, *B; + TACSAssembler *assembler; + TACSBVec *temp; +}; + /* Implementation of the buckling/frequency analysis */ @@ -827,7 +877,7 @@ void TACSFrequencyAnalysis::solve(KSMPrint *ksm_print, int print_level) { t0 = MPI_Wtime() - t0; char line[256]; - sprintf(line, "JD computational time: %15.6f\n", t0); + snprintf(line, sizeof(line), "JD computational time: %15.6f\n", t0); ksm_print->print(line); } } else { @@ -873,7 +923,7 @@ void TACSFrequencyAnalysis::solve(KSMPrint *ksm_print, int print_level) { t0 = MPI_Wtime() - t0; char line[256]; - sprintf(line, "Lanczos computational time: %15.6f\n", t0); + snprintf(line, sizeof(line), "Lanczos computational time: %15.6f\n", t0); ksm_print->print(line); } } @@ -954,6 +1004,7 @@ void TACSFrequencyAnalysis::evalEigenDVSens(int n, TACSBVec *dfdx) { eigvec, eigvec, dfdx); // Finish computing the derivative + // Note - this shouldn't be needed since the eigenvalues are M-orthogonal if (mmat) { mmat->mult(eigvec, res); } @@ -1011,6 +1062,168 @@ void TACSFrequencyAnalysis::evalEigenXptSens(int n, TACSBVec *dfdXpt) { dfdXpt->scale(scale); } +/* + Compute the derivative a function of interest with respect to design variables + given the derivative of the eigenvalues and eigenvectors with respect to a + function of interest. +*/ +void TACSFrequencyAnalysis::addEigenSens(int neigs, TacsScalar dfdlam[], + TACSBVec *dfdq[], TACSBVec *dfdx, + TACSBVec *dfdXpt, int use_cg, + double rtol, double atol) { + TacsScalar *eigs = new TacsScalar[neigs]; + TACSVec **Cvecs = new TACSVec *[neigs]; + TACSBVec **eigvecs = new TACSBVec *[neigs]; + + // Assemble the mass matrix + assembler->assembleMatType(TACS_MASS_MATRIX, mmat); + + for (int i = 0; i < neigs; i++) { + // Create space for the eigenvector + eigvecs[i] = assembler->createVec(); + eigvecs[i]->incref(); + + // Extract the eigenvalue and eigenvector + TacsScalar error; + eigs[i] = extractEigenvector(i, eigvecs[i], &error); + + // Now, create the constraint space vector + TACSBVec *vec = assembler->createVec(); + vec->incref(); + + // Form the matrix-vector product C = M * eigvec + mmat->mult(eigvecs[i], vec); + Cvecs[i] = vec; + } + + // Allocate the constraint + TACSKSMConstraint *con = new TACSKSMConstraint(neigs, Cvecs); + con->incref(); + + if (mg) { + assembler->assembleMatType(TACS_STIFFNESS_MATRIX, kmat); + + // Factor the matrix - associated with the stiffness matrix + pc->factor(); + } else { + // Assemble a linear combination of the mass and stiffness matrices + // kmat = K - fact * eigs[0] * M + double fact = 0.99; // This should be a parameter + ElementMatrixType kmat_types[2] = {TACS_STIFFNESS_MATRIX, TACS_MASS_MATRIX}; + TacsScalar kscale[2]; + kscale[0] = 1.0; + kscale[1] = 0.0; // -fact * eigs[0]; + assembler->assembleMatCombo(kmat_types, kscale, 2, kmat); + + // Factor the matrix + pc->factor(); + + // Now assemble the stiffness and mass matrices again + assembler->assembleMatType(TACS_STIFFNESS_MATRIX, kmat); + } + + // Create the combo matrix operator + TACSMatCombo *combo = new TACSMatCombo(1.0, 0.0, kmat, mmat, assembler); + combo->incref(); + + TACSKsm *ksm = NULL; + + if (use_cg) { + ksm = new PCG(combo, pc, 100, 10, con); + ksm->incref(); + + int mpi_rank; + MPI_Comm_rank(assembler->getMPIComm(), &mpi_rank); + ksm->setMonitor(new KSMPrintStdout("PCG", mpi_rank, 1)); + } else { + int gmres_iters = 25; + int nrestart = 5; + int is_flexible = 0; + ksm = new GMRES(combo, pc, gmres_iters, nrestart, is_flexible, con); + ksm->incref(); + + int mpi_rank; + MPI_Comm_rank(assembler->getMPIComm(), &mpi_rank); + ksm->setMonitor(new KSMPrintStdout("GMRES", mpi_rank, 1)); + } + + ksm->setTolerances(rtol, atol); + + // Set the adjoint vector + TACSBVec *psi = res; + + for (int i = 0; i < neigs; i++) { + // Assemble the linear system mmat = K - eigs[i] * M + combo->setCoefficients(1.0, -eigs[i]); + + // Solve the constrained problem + ksm->solve(dfdq[i], psi); + psi->scale(-1.0); + + for (int j = 0; j < neigs; j++) { + if (i != j) { + TacsScalar theta = dfdq[i]->dot(eigvecs[j]) / (eigs[i] - eigs[j]); + psi->axpy(theta, eigvecs[j]); + } + } + + // Add the contribution from the derivative of the eigenvalues + psi->axpy(dfdlam[i], eigvecs[i]); + + // Add the material derivatives + if (dfdx) { + assembler->addMatDVSensInnerProduct(1.0, TACS_STIFFNESS_MATRIX, psi, + eigvecs[i], dfdx); + assembler->addMatDVSensInnerProduct(-eigs[i], TACS_MASS_MATRIX, psi, + eigvecs[i], dfdx); + } + + // Add the coordinate derivatives + if (dfdXpt) { + assembler->addMatXptSensInnerProduct(1.0, TACS_STIFFNESS_MATRIX, psi, + eigvecs[i], dfdXpt); + assembler->addMatXptSensInnerProduct(-eigs[i], TACS_MASS_MATRIX, psi, + eigvecs[i], dfdXpt); + } + + // Add the contribution from alpha + TacsScalar alpha = -0.5 * dfdq[i]->dot(eigvecs[i]); + if (dfdx) { + assembler->addMatDVSensInnerProduct(alpha, TACS_MASS_MATRIX, eigvecs[i], + eigvecs[i], dfdx); + } + + // Add the coordinate derivatives + if (dfdXpt) { + assembler->addMatXptSensInnerProduct(alpha, TACS_MASS_MATRIX, eigvecs[i], + eigvecs[i], dfdXpt); + } + } + + if (dfdx) { + dfdx->beginSetValues(TACS_ADD_VALUES); + dfdx->endSetValues(TACS_ADD_VALUES); + } + + if (dfdXpt) { + dfdXpt->beginSetValues(TACS_ADD_VALUES); + dfdXpt->endSetValues(TACS_ADD_VALUES); + } + + // Free memory + for (int i = 0; i < neigs; i++) { + eigvecs[i]->decref(); + Cvecs[i]->decref(); + } + delete[] eigvecs; + delete[] Cvecs; + delete[] eigs; + + combo->decref(); + con->decref(); + ksm->decref(); +} + /*! Check the actual residual for the given eigenvalue */ diff --git a/src/TACSBuckling.h b/src/TACSBuckling.h index 83de19bda..3a6ab8b56 100644 --- a/src/TACSBuckling.h +++ b/src/TACSBuckling.h @@ -149,8 +149,6 @@ class TACSFrequencyAnalysis : public TACSObject { TacsScalar getSigma(); void setSigma(TacsScalar _sigma); void solve(KSMPrint *ksm_print = NULL, int print_level = 0); - void evalEigenDVSens(int n, TACSBVec *dfdx); - void evalEigenXptSens(int n, TACSBVec *dfdX); // Extract and check the solution // ------------------------------ @@ -159,6 +157,18 @@ class TACSFrequencyAnalysis : public TACSObject { void checkEigenvector(int n); TacsScalar checkOrthogonality(); + // Compute derivatives of the eigenvalues wrt design variables and coordinates + // --------------------------------------------------------------------------- + void evalEigenDVSens(int n, TACSBVec *dfdx); + void evalEigenXptSens(int n, TACSBVec *dfdX); + + // Given the derivative of a function wrt eigenvalues and eigenvectors, + // compute the derivative wrt the design variabels and coordinates + // ---------------------------------------------------------------------------- + void addEigenSens(int neigs, TacsScalar dfdlam[], TACSBVec *dfdq[], + TACSBVec *dfdx = NULL, TACSBVec *dfdXpts = NULL, + int use_cg = 0, double rtol = 1e-8, double atol = 1e-30); + private: // The TACS assembler object TACSAssembler *assembler; diff --git a/src/TACSContinuation.cpp b/src/TACSContinuation.cpp index 2314ce2db..ea4206aca 100644 --- a/src/TACSContinuation.cpp +++ b/src/TACSContinuation.cpp @@ -310,7 +310,7 @@ void TACSContinuation::solve_tangent(TACSMat *mat, TACSPc *pc, TACSKsm *ksm, if (ksm_print) { char line[256]; ksm_print->print("Performing initial Newton iterations\n"); - sprintf(line, "%5s %9s %10s\n", "Iter", "t", "|R|"); + snprintf(line, sizeof(line), "%5s %9s %10s\n", "Iter", "t", "|R|"); ksm_print->print(line); } @@ -343,8 +343,8 @@ void TACSContinuation::solve_tangent(TACSMat *mat, TACSPc *pc, TACSKsm *ksm, TacsScalar res_norm = res->norm(); if (ksm_print) { char line[256]; - sprintf(line, "%5d %9.4f %10.4e\n", k + 1, MPI_Wtime() - t0, - TacsRealPart(res_norm)); + snprintf(line, sizeof(line), "%5d %9.4f %10.4e\n", k + 1, + MPI_Wtime() - t0, TacsRealPart(res_norm)); ksm_print->print(line); } if (k == 0) { @@ -469,11 +469,12 @@ void TACSContinuation::solve_tangent(TACSMat *mat, TACSPc *pc, TACSKsm *ksm, if (ksm_print) { char line[256]; - sprintf(line, "Outer iteration %3d: t: %9.4f dp_ds: %10.4e\n", - iteration_count, MPI_Wtime() - t0, TacsRealPart(dlambda_ds)); + snprintf(line, sizeof(line), + "Outer iteration %3d: t: %9.4f dp_ds: %10.4e\n", iteration_count, + MPI_Wtime() - t0, TacsRealPart(dlambda_ds)); ksm_print->print(line); - sprintf(line, "%5s %9s %10s %10s %10s\n", "Iter", "t", "|R|", "lambda", - "|u|"); + snprintf(line, sizeof(line), "%5s %9s %10s %10s %10s\n", "Iter", "t", + "|R|", "lambda", "|u|"); ksm_print->print(line); } @@ -507,9 +508,9 @@ void TACSContinuation::solve_tangent(TACSMat *mat, TACSPc *pc, TACSKsm *ksm, TacsScalar res_norm = res->norm(); if (ksm_print) { char line[256]; - sprintf(line, "%5d %9.4f %10.3e %10.3e %10.3e\n", j, MPI_Wtime() - t0, - TacsRealPart(res_norm), TacsRealPart(lambda), - TacsRealPart(vars->norm())); + snprintf(line, sizeof(line), "%5d %9.4f %10.3e %10.3e %10.3e\n", j, + MPI_Wtime() - t0, TacsRealPart(res_norm), + TacsRealPart(lambda), TacsRealPart(vars->norm())); ksm_print->print(line); } @@ -542,9 +543,9 @@ void TACSContinuation::solve_tangent(TACSMat *mat, TACSPc *pc, TACSKsm *ksm, if (ksm_print) { char line[256]; - sprintf(line, - "Failed to converge, retrying with step size = %10.3e\n", - TacsRealPart(delta_s)); + snprintf(line, sizeof(line), + "Failed to converge, retrying with step size = %10.3e\n", + TacsRealPart(delta_s)); ksm_print->print(line); } } diff --git a/src/TACSIntegrator.cpp b/src/TACSIntegrator.cpp index f2e6fbbd7..1a773517b 100644 --- a/src/TACSIntegrator.cpp +++ b/src/TACSIntegrator.cpp @@ -35,7 +35,7 @@ TACSIntegrator::TACSIntegrator(TACSAssembler *_assembler, double tinit, assembler->incref(); // Set the default prefix = results - sprintf(prefix, "./"); + snprintf(prefix, sizeof(prefix), "./"); // Allocate the total number of time steps num_time_steps = int(num_steps); @@ -352,7 +352,7 @@ void TACSIntegrator::setFunctions(int _num_funcs, TACSFunction **_funcs, Set the output directory prefix */ void TACSIntegrator::setOutputPrefix(const char *_prefix) { - strncpy(prefix, _prefix, sizeof(prefix)); + strncpy(prefix, _prefix, sizeof(prefix) - 1); } /* @@ -483,7 +483,7 @@ void TACSIntegrator::writeStepToF5(int step_num) { // Write the f5 file if (f5) { - char fname[256]; + char fname[512]; snprintf(fname, sizeof(fname), "%s/output_%06d.f5", prefix, step_num); f5->writeToFile(fname); } @@ -926,7 +926,7 @@ int TACSIntegrator::lapackNaturalFrequencies(int use_gyroscopic, TACSBVec *q, // Write the mode to disk as f5 assembler->setVariables(mode, mode, mode); if (f5) { - char fname[256]; + char fname[512]; snprintf(fname, sizeof(fname), "%s/mode_freq_%d.f5", prefix, index); f5->writeToFile(fname); } @@ -991,7 +991,7 @@ int TACSIntegrator::lapackNaturalFrequencies(int use_gyroscopic, TACSBVec *q, // Write the mode to disk as f5 assembler->setVariables(mode, mode, mode); if (f5) { - char fname[256]; + char fname[512]; snprintf(fname, sizeof(fname), "%s/mode_freq_%d.f5", prefix, index); f5->writeToFile(fname); } diff --git a/src/TACSMg.cpp b/src/TACSMg.cpp index ace0c636c..c831bf373 100644 --- a/src/TACSMg.cpp +++ b/src/TACSMg.cpp @@ -692,8 +692,9 @@ void TACSMg::applyFactor(TACSVec *bvec, TACSVec *xvec) { if (monitor) { for (int k = 0; k < nlevels; k++) { char descript[128]; - sprintf(descript, "TACSMg cumulative level %2d time %15.8e\n", k, - cumulative_level_time[k]); + snprintf(descript, sizeof(descript), + "TACSMg cumulative level %2d time %15.8e\n", k, + cumulative_level_time[k]); monitor->print(descript); } } diff --git a/src/bpmat/GSEP.cpp b/src/bpmat/GSEP.cpp index b26dd3eae..e22b749dd 100644 --- a/src/bpmat/GSEP.cpp +++ b/src/bpmat/GSEP.cpp @@ -572,25 +572,25 @@ void SEP::solve(KSMPrint *ksm_print, KSMPrint *ksm_file) { // Print out a summary of the eigenvalues and errors if (ksm_print) { char line[256]; - sprintf(line, "%3s %18s %18s %10s\n", " ", "eigenvalue", "shift-invert eig", - "error"); + snprintf(line, sizeof(line), "%3s %18s %18s %10s\n", " ", "eigenvalue", + "shift-invert eig", "error"); ksm_print->print(line); for (int i = 0; i < niters; i++) { int index = perm[i]; char line[256]; - sprintf(line, "%3d %18.10e %18.10e %10.3e\n", i, - TacsRealPart(Op->convertEigenvalue(eigs[index])), - TacsRealPart(eigs[index]), - fabs(TacsRealPart(Beta[niters - 1] * - eigvecs[index * niters + (niters - 1)] * er))); + snprintf(line, sizeof(line), "%3d %18.10e %18.10e %10.3e\n", i, + TacsRealPart(Op->convertEigenvalue(eigs[index])), + TacsRealPart(eigs[index]), + fabs(TacsRealPart(Beta[niters - 1] * + eigvecs[index * niters + (niters - 1)] * er))); ksm_print->print(line); } } // Print the iteration count to file if (ksm_file) { char line[256]; - sprintf(line, "%2d\n", niters); + snprintf(line, sizeof(line), "%2d\n", niters); ksm_file->print(line); } } diff --git a/src/bpmat/JacobiDavidson.cpp b/src/bpmat/JacobiDavidson.cpp index d90e08570..b1c3b974d 100644 --- a/src/bpmat/JacobiDavidson.cpp +++ b/src/bpmat/JacobiDavidson.cpp @@ -555,8 +555,8 @@ void TACSJacobiDavidson::solve(KSMPrint *ksm_print, int print_level) { if (ksm_print && print_level > 0) { char line[256]; - sprintf(line, "%4s %15s %15s %10s\n", "Iter", "JD Residual", "Ritz value", - "toler"); + snprintf(line, sizeof(line), "%4s %15s %15s %10s\n", "Iter", "JD Residual", + "Ritz value", "toler"); ksm_print->print(line); } @@ -653,8 +653,8 @@ void TACSJacobiDavidson::solve(KSMPrint *ksm_print, int print_level) { if (ksm_print && print_level > 0) { char line[256]; - sprintf(line, "%4d %15.5e %15.5e %10.2e\n", iteration, - TacsRealPart(w_norm), theta, toler); + snprintf(line, sizeof(line), "%4d %15.5e %15.5e %10.2e\n", iteration, + TacsRealPart(w_norm), theta, toler); ksm_print->print(line); } @@ -897,15 +897,17 @@ void TACSJacobiDavidson::solve(KSMPrint *ksm_print, int print_level) { char line[256]; for (int i = 0; i < nconverged; i++) { int index = eigindex[i]; - sprintf(line, "Eigenvalue[%2d]: %25.10e Eig. error[%2d]: %25.10e\n", i, - TacsRealPart(eigvals[index]), i, TacsRealPart(eigerror[index])); + snprintf(line, sizeof(line), + "Eigenvalue[%2d]: %25.10e Eig. error[%2d]: %25.10e\n", i, + TacsRealPart(eigvals[index]), i, TacsRealPart(eigerror[index])); ksm_print->print(line); } - sprintf(line, "JD number of outer iterations: %2d\n", iteration); + snprintf(line, sizeof(line), "JD number of outer iterations: %2d\n", + iteration); ksm_print->print(line); - sprintf(line, "JD number of inner GMRES iterations: %2d\n", - gmres_iteration); + snprintf(line, sizeof(line), "JD number of inner GMRES iterations: %2d\n", + gmres_iteration); ksm_print->print(line); } diff --git a/src/bpmat/KSM.cpp b/src/bpmat/KSM.cpp index fc10c832a..40fd40fce 100644 --- a/src/bpmat/KSM.cpp +++ b/src/bpmat/KSM.cpp @@ -21,6 +21,35 @@ #include #include +/* + Classical Gram-Schmidt orthogonalization. + + Given an input vector q, and the set of orthogonal vectors w. Make q + orthogonal to the set of vectors w. + + Return an array h = w^{T} q. + q' = q - w h + w^{T} q' = w^{T} q - w^{T} w w^{T} q +*/ +static void ClassicalGramSchmidt(TacsScalar *h, TACSVec *q, TACSVec **w, + int nvecs) { + q->mdot(w, h, nvecs); + for (int j = 0; j < nvecs; j++) { + q->axpy(-h[j], w[j]); + } +} + +/* + Modified Gram-Schmidt orthogonalization +*/ +static void ModifiedGramSchmidt(TacsScalar *h, TACSVec *q, TACSVec **w, + int nvecs) { + for (int j = 0; j < nvecs; j++) { + h[j] = w[j]->dot(q); + q->axpy(-h[j], w[j]); + } +} + /* Implementation of various Krylov-subspace methods */ @@ -332,6 +361,45 @@ void KSMPrintFile::print(const char *cstr) { } } +TACSKSMConstraint::TACSKSMConstraint(int _nvecs, TACSVec *_vecs[], + OrthogonalizationType otype) { + nvecs = _nvecs; + alpha = new TacsScalar[nvecs]; + vecs = new TACSVec *[nvecs]; + for (int i = 0; i < nvecs; i++) { + vecs[i] = _vecs[i]; + vecs[i]->incref(); + } + + if (otype == CLASSICAL_GRAM_SCHMIDT) { + ortho = ClassicalGramSchmidt; + } else { + ortho = ModifiedGramSchmidt; + } + + orthogonalize(); +} + +TACSKSMConstraint::~TACSKSMConstraint() { + for (int i = 0; i < nvecs; i++) { + vecs[i]->decref(); + } + delete[] vecs; + delete[] alpha; +} + +void TACSKSMConstraint::orthogonalize() { + for (int i = 0; i < nvecs; i++) { + // Orthogonalize this vector against all previous vectors + ortho(alpha, vecs[i], vecs, i); + + TacsScalar nrm = vecs[i]->norm(); + vecs[i]->scale(1.0 / nrm); + } +} + +void TACSKSMConstraint::apply(TACSVec *x) { ortho(alpha, x, vecs, nvecs); } + /* The preconditioned conjugate gradient method @@ -344,14 +412,19 @@ void KSMPrintFile::print(const char *cstr) { reset: reset the CG iterations every 'reset' iterations nouter: the number of resets to try before giving up */ -PCG::PCG(TACSMat *_mat, TACSPc *_pc, int _reset, int _nouter) { +PCG::PCG(TACSMat *_mat, TACSPc *_pc, int _reset, int _nouter, + TACSKSMConstraint *_con) { monitor = NULL; mat = _mat; pc = _pc; + con = _con; mat->incref(); pc->incref(); + if (con) { + con->incref(); + } reset = _reset; nouter = _nouter; @@ -375,6 +448,9 @@ PCG::PCG(TACSMat *_mat, TACSPc *_pc, int _reset, int _nouter) { PCG::~PCG() { mat->decref(); pc->decref(); + if (con) { + con->decref(); + } work->decref(); R->decref(); @@ -428,6 +504,9 @@ void PCG::setMonitor(KSMPrint *_monitor) { Solve the linear system with the preconditioned conjugate gradient method + This uses the variant of PCG from the paper "Inexact Preconditioned + Conjugate Gradient Method with Inner-Outer Iteration" by Golub and Ye. + input: b: the right-hand-side x: the solution vector @@ -454,6 +533,11 @@ int PCG::solve(TACSVec *b, TACSVec *x, int zero_guess) { R->axpby(1.0, -1.0, b); // R = b - A*x } + if (con) { + con->apply(R); + con->apply(x); + } + if (count == 0) { rhs_norm = R->norm(); resNorm = rhs_norm; @@ -464,21 +548,40 @@ int PCG::solve(TACSVec *b, TACSVec *x, int zero_guess) { } if (TacsRealPart(rhs_norm) > atol) { - // Apply the preconditioner + // Apply the preconditioner Z = M^{-1} * R pc->applyFactor(R, Z); + if (con) { + con->apply(Z); + } // Copy Z to P, ie. P = Z P->copyValues(Z); + // Compute rz = (R, Z) + TacsScalar rz = R->dot(Z); + for (int i = 0; i < reset; i++) { - mat->mult(P, work); // work = A*P - TacsScalar temp = R->dot(Z); // (R,Z) - TacsScalar alpha = temp / (work->dot(P)); // alpha = (R,Z)/(A*P,P) - x->axpy(alpha, P); // x = x + alpha*P - R->axpy(-alpha, work); // R' = R - alpha*A*P - pc->applyFactor(R, Z); // Z' = M^{-1} R - TacsScalar beta = R->dot(Z) / temp; // beta = (R',Z')/(R,Z) - P->axpby(1.0, beta, Z); // P' = Z' + beta*P + mat->mult(P, work); // work = A * P + if (con) { // work = (I - Q * Q^{T}) * work + con->apply(work); + } + + TacsScalar alpha = rz / (work->dot(P)); // alpha = (R, Z)/(A * P, P) + x->axpy(alpha, P); // x = x + alpha * P + R->axpy(-alpha, work); // R' = R - alpha * A * P + + pc->applyFactor(R, work); // Z' = M^{-1} R + if (con) { // Z' = (I - Q * Q^{T}) * Z' + con->apply(work); + } + + TacsScalar rz_new = R->dot(work); // rz_new = (R', Z') + TacsScalar rz_old = R->dot(Z); // rz_old = (R', Z) + TacsScalar beta = + (rz_new - rz_old) / rz; // beta = ((R',Z') - rz0) / (R,Z) + P->axpby(1.0, beta, work); // P' = Z' + beta*P + Z->copyValues(work); // Z <- work + rz = rz_new; // rz <- (R', Z') iterCount++; resNorm = R->norm(); @@ -502,35 +605,6 @@ int PCG::solve(TACSVec *b, TACSVec *x, int zero_guess) { return solve_flag; } -/* - Classical Gram-Schmidt orthogonalization. - - Given an input vector q, and the set of orthogonal vectors w. Make q - orthogonal to the set of vectors w. - - Return an array h = w^{T} q. - q' = q - w h - w^{T} q' = w^{T} q - w^{T} w w^{T} q -*/ -static void ClassicalGramSchmidt(TacsScalar *h, TACSVec *q, TACSVec **w, - int nvecs) { - q->mdot(w, h, nvecs); - for (int j = 0; j < nvecs; j++) { - q->axpy(-h[j], w[j]); - } -} - -/* - Modified Gram-Schmidt orthogonalization -*/ -static void ModifiedGramSchmidt(TacsScalar *h, TACSVec *q, TACSVec **w, - int nvecs) { - for (int j = 0; j < nvecs; j++) { - h[j] = w[j]->dot(q); - q->axpy(-h[j], w[j]); - } -} - /* Create a GMRES object for solving a linear system. @@ -544,9 +618,9 @@ static void ModifiedGramSchmidt(TacsScalar *h, TACSVec *q, TACSVec **w, nrestart: the number of restarts before we give up isFlexible: is the preconditioner actually flexible? If so use FGMRES */ -GMRES::GMRES(TACSMat *_mat, TACSPc *_pc, int _m, int _nrestart, - int _isFlexible) { - init(_mat, _pc, _m, _nrestart, _isFlexible); +GMRES::GMRES(TACSMat *_mat, TACSPc *_pc, int _m, int _nrestart, int _isFlexible, + TACSKSMConstraint *_con) { + init(_mat, _pc, _m, _nrestart, _isFlexible, _con); } /* @@ -558,8 +632,8 @@ GMRES::GMRES(TACSMat *_mat, TACSPc *_pc, int _m, int _nrestart, m: the size of the Krylov subspace nrestart: try this many times before giving up */ -GMRES::GMRES(TACSMat *_mat, int _m, int _nrestart) { - init(_mat, NULL, _m, _nrestart, 0); +GMRES::GMRES(TACSMat *_mat, int _m, int _nrestart, TACSKSMConstraint *_con) { + init(_mat, NULL, _m, _nrestart, 0, _con); } /* @@ -568,7 +642,7 @@ GMRES::GMRES(TACSMat *_mat, int _m, int _nrestart) { This is called by both of the two constructors above. */ void GMRES::init(TACSMat *_mat, TACSPc *_pc, int _m, int _nrestart, - int _isFlexible) { + int _isFlexible, TACSKSMConstraint *_con) { orthogonalize = ModifiedGramSchmidt; monitor = NULL; monitor_time = 0; @@ -578,8 +652,12 @@ void GMRES::init(TACSMat *_mat, TACSPc *_pc, int _m, int _nrestart, mat = _mat; pc = _pc; - mat->incref(); + con = _con; + mat->incref(); + if (con) { + con->incref(); + } if (pc) { pc->incref(); } else { @@ -648,6 +726,9 @@ GMRES::~GMRES() { if (pc) { pc->decref(); } + if (con) { + con->decref(); + } if (work) { work->decref(); } @@ -750,7 +831,7 @@ void GMRES::setMonitor(KSMPrint *_monitor) { Unless you have a good reason, you should use modified Gram-Schmidt. */ -void GMRES::setOrthoType(enum OrthoType otype) { +void GMRES::setOrthoType(OrthogonalizationType otype) { if (otype == CLASSICAL_GRAM_SCHMIDT) { orthogonalize = ClassicalGramSchmidt; } else { @@ -800,13 +881,22 @@ int GMRES::solve(TACSVec *b, TACSVec *x, int zero_guess) { // If the initial guess is zero x->zeroEntries(); // Set x = 0 W[0]->copyValues(b); // W[0] = b + if (con) { // Compute W[0] = (I - Q * Q^{T}) * W[0] + con->apply(W[0]); + } res[0] = W[0]->norm(); W[0]->scale(1.0 / res[0]); // W[0] = b/|| b || } else { // If the initial guess is non-zero or restarting + if (con) { + con->apply(x); + } mat->mult(x, W[0]); W[0]->axpy(-1.0, b); // W[0] = A*x - b + if (con) { // Compute W[0] = (I - Q * Q^{T}) * W[0] + con->apply(W[0]); + } res[0] = W[0]->norm(); W[0]->scale(-1.0 / res[0]); // W[0] = (b - A*x)/|| b - A*x || @@ -835,14 +925,29 @@ int GMRES::solve(TACSVec *b, TACSVec *x, int zero_guess) { if (isFlexible) { // Apply the preconditioner, Z[i] = M^{-1} W[i] pc->applyFactor(W[i], Z[i]); + if (con) { + con->apply(Z[i]); + } mat->mult(Z[i], W[i + 1]); // W[i+1] = A*Z[i] = A*M^{-1}*W[i] + if (con) { + con->apply(W[i + 1]); + } } else { if (pc) { // Apply the preconditioner, work = M^{-1} W[i] pc->applyFactor(W[i], work); + if (con) { + con->apply(work); + } mat->mult(work, W[i + 1]); // W[i+1] = A*work = A*M^{-1}*W[i] + if (con) { + con->apply(W[i + 1]); + } } else { mat->mult(W[i], W[i + 1]); // Compute W[i+1] = A*W[i] + if (con) { + con->apply(W[i + 1]); + } } } if (monitor_time) { @@ -934,6 +1039,9 @@ int GMRES::solve(TACSVec *b, TACSVec *x, int zero_guess) { // Apply M^{-1} to the linear combination pc->applyFactor(work, W[0]); x->axpy(1.0, W[0]); + if (con) { + con->apply(x); + } } if (solve_flag) { @@ -944,9 +1052,9 @@ int GMRES::solve(TACSVec *b, TACSVec *x, int zero_guess) { if (monitor_time && monitor) { t_total = MPI_Wtime() - t_total; char str_mat[80], str_ort[80], str_tot[80]; - sprintf(str_mat, "pc-mat time %10.6f\n", t_pc); - sprintf(str_ort, "ortho time %10.6f\n", t_ortho); - sprintf(str_tot, "total time %10.6f\n", t_total); + snprintf(str_mat, sizeof(str_mat), "pc-mat time %10.6f\n", t_pc); + snprintf(str_ort, sizeof(str_ort), "ortho time %10.6f\n", t_ortho); + snprintf(str_tot, sizeof(str_tot), "total time %10.6f\n", t_total); monitor->print(str_mat); monitor->print(str_ort); monitor->print(str_tot); diff --git a/src/bpmat/KSM.h b/src/bpmat/KSM.h index c304ecdb2..96f24a643 100644 --- a/src/bpmat/KSM.h +++ b/src/bpmat/KSM.h @@ -33,6 +33,11 @@ */ enum MatrixOrientation { TACS_MAT_NORMAL, TACS_MAT_TRANSPOSE }; +/* + The type of orthogonalization to use +*/ +enum OrthogonalizationType { CLASSICAL_GRAM_SCHMIDT, MODIFIED_GRAM_SCHMIDT }; + /* Define boundary conditions that are applied after all the matrix/vector values have been set. @@ -237,6 +242,46 @@ class KSMPrint : public TACSObject { static const char *printName; }; +/* + Constraint class for KSM methods + + Construct a constraint for Krylov methods that takes the form + + V^{T} * x = 0 + + This is enforced by computing the action of the operator + + x = (I - V * (V^{T} * V)^{-1} * V^{T}) * x + + The vectors are orthogonalized using Gram Schmidt such that + + V = Q * R + + As a result, the operator can be computed as + + x = (I - Q * Q^{T}) * x +*/ +class TACSKSMConstraint : public TACSObject { + public: + TACSKSMConstraint(int _nvecs, TACSVec *_vecs[], + OrthogonalizationType otype = MODIFIED_GRAM_SCHMIDT); + ~TACSKSMConstraint(); + + void orthogonalize(); + void apply(TACSVec *x); + + private: + // Orthogonalize a vector against a set of vectors + void (*ortho)(TacsScalar *, TACSVec *, TACSVec **, int); + + // The size of the constraint space + int nvecs; + TACSVec **vecs; + + // Temporary scalar values + TacsScalar *alpha; +}; + /*! The abstract Krylov-subspace method class @@ -327,7 +372,8 @@ class KSMPrintFile : public KSMPrint { */ class PCG : public TACSKsm { public: - PCG(TACSMat *_mat, TACSPc *_pc, int reset, int _nouter); + PCG(TACSMat *_mat, TACSPc *_pc, int reset, int _nouter, + TACSKSMConstraint *con = NULL); ~PCG(); TACSVec *createVec() { return mat->createVec(); } @@ -342,6 +388,9 @@ class PCG : public TACSKsm { TACSMat *mat; TACSPc *pc; + // Constraints (if any) + TACSKSMConstraint *con; + // The relative/absolute tolerances double rtol, atol; @@ -390,9 +439,9 @@ class PCG : public TACSKsm { */ class GMRES : public TACSKsm { public: - enum OrthoType { CLASSICAL_GRAM_SCHMIDT, MODIFIED_GRAM_SCHMIDT }; - GMRES(TACSMat *_mat, TACSPc *_pc, int _m, int _nrestart, int _isFlexible); - GMRES(TACSMat *_mat, int _m, int _nrestart); + GMRES(TACSMat *_mat, TACSPc *_pc, int _m, int _nrestart, int _isFlexible, + TACSKSMConstraint *con = NULL); + GMRES(TACSMat *_mat, int _m, int _nrestart, TACSKSMConstraint *con = NULL); ~GMRES(); TACSVec *createVec() { return mat->createVec(); } @@ -401,19 +450,21 @@ class GMRES : public TACSKsm { void getOperators(TACSMat **_mat, TACSPc **_pc); void setTolerances(double _rtol, double _atol); void setMonitor(KSMPrint *_monitor); - void setOrthoType(enum OrthoType otype); + void setOrthoType(OrthogonalizationType otype); void setTimeMonitor(); const char *getObjectName(); private: // Initialize the class - void init(TACSMat *_mat, TACSPc *_pc, int _m, int _nrestart, int _isFlexible); + void init(TACSMat *_mat, TACSPc *_pc, int _m, int _nrestart, int _isFlexible, + TACSKSMConstraint *_con); // Orthogonalize a vector against a set of vectors void (*orthogonalize)(TacsScalar *, TACSVec *, TACSVec **, int); TACSMat *mat; TACSPc *pc; + TACSKSMConstraint *con; int msub; int nrestart; int isFlexible; diff --git a/src/constitutive/TACSPanelAnalysis.cpp b/src/constitutive/TACSPanelAnalysis.cpp index 6ec2f3020..024692506 100644 --- a/src/constitutive/TACSPanelAnalysis.cpp +++ b/src/constitutive/TACSPanelAnalysis.cpp @@ -1717,7 +1717,8 @@ int TACSPanelAnalysis::computeBucklingLoads(TacsScalar Nx, TacsScalar Nxy, char *file_name = new char[file_len]; for (int i = 0; i < nloads; i++) { - sprintf(file_name, "%sbuckling_mode%02d.dat", prefix, i); + snprintf(file_name, sizeof(file_name), "%sbuckling_mode%02d.dat", prefix, + i); printPanelMode(file_name, &eigvecs[nvars * i], 125); } @@ -2176,7 +2177,7 @@ int TACSPanelAnalysis::computeFrequencies(TacsScalar freq[], int nfreq, char *file_name = new char[file_len]; for (int i = 0; i < nfreq; i++) { - sprintf(file_name, "%spanel_mode%02d.dat", prefix, i); + snprintf(file_name, sizeof(file_name), "%spanel_mode%02d.dat", prefix, i); printPanelMode(file_name, &eigvecs[nvars * i], 125); } diff --git a/src/elements/dynamics/MITC3.cpp b/src/elements/dynamics/MITC3.cpp index e5b21c756..658c77550 100644 --- a/src/elements/dynamics/MITC3.cpp +++ b/src/elements/dynamics/MITC3.cpp @@ -3068,7 +3068,7 @@ void MITC3::testStrain(const TacsScalar X[]) { // Write out the error components char descript[64]; - sprintf(descript, "strain after rigid rotation"); + snprintf(descript, sizeof(descript), "strain after rigid rotation"); writeErrorComponents(stdout, descript, e, fd, 6); // Compute the bmatrix @@ -3111,7 +3111,7 @@ void MITC3::testStrain(const TacsScalar X[]) { // Write out the error components char descript[64]; - sprintf(descript, "B%d", k); + snprintf(descript, sizeof(descript), "B%d", k); writeErrorComponents(stdout, descript, &B[6 * k], fd, 6); } @@ -3148,6 +3148,6 @@ void MITC3::testStrain(const TacsScalar X[]) { getStrain(&u, X, vars, fd); // Write out the error components of the strain - sprintf(descript, "strain before/after rigid rotation"); + snprintf(descript, sizeof(descript), "strain before/after rigid rotation"); writeErrorComponents(stdout, descript, e, fd, 6); } diff --git a/src/elements/dynamics/MITC9.cpp b/src/elements/dynamics/MITC9.cpp index cfc6525f2..7a758e7b1 100644 --- a/src/elements/dynamics/MITC9.cpp +++ b/src/elements/dynamics/MITC9.cpp @@ -6142,7 +6142,7 @@ void MITC9::testStrain(const TacsScalar X[]) { // Write out the error components char descript[64]; - sprintf(descript, "B%d", k); + snprintf(descript, sizeof(descript), "B%d", k); writeErrorComponents(stdout, descript, &B[8 * k], fd, 8); } } diff --git a/src/elements/dynamics/TACSRigidBody.cpp b/src/elements/dynamics/TACSRigidBody.cpp index fb9eb9f41..7b0b7dac5 100644 --- a/src/elements/dynamics/TACSRigidBody.cpp +++ b/src/elements/dynamics/TACSRigidBody.cpp @@ -1517,7 +1517,7 @@ void TACSRigidBody::testJacobian(double dh, TacsScalar alpha, TacsScalar beta, // Print out the results to stdout char outname[128]; - sprintf(outname, "Jacobian col %d", ii); + snprintf(outname, sizeof(outname), "Jacobian col %d", ii); writeErrorComponents(stdout, outname, res, fd, 8); } } diff --git a/src/elements/shell/TACSDirector.h b/src/elements/shell/TACSDirector.h index ed1bd51ed..37484c174 100644 --- a/src/elements/shell/TACSDirector.h +++ b/src/elements/shell/TACSDirector.h @@ -1060,7 +1060,6 @@ class TACSQuadraticRotation { TacsScalar res[]) { TacsScalar *r = &res[offset]; const TacsScalar *q = &vars[offset]; - const TacsScalar *qdot = &dvars[offset]; for (int i = 0; i < num_nodes; i++) { // D = (t^{x} + 0.5*(q^{x}t)^{x} - 0.5*t^{x}*q^{x}) @@ -1085,7 +1084,6 @@ class TACSQuadraticRotation { r += vars_per_node; q += vars_per_node; - qdot += vars_per_node; dd += 3; t += 3; @@ -1154,8 +1152,6 @@ class TACSQuadraticRotation { Diddot[8] = 0.5 * (qiddot[0] * ti[0] + qiddot[1] * ti[1]); } - const TacsScalar *q = &vars[offset]; - const TacsScalar *qdot = &dvars[offset]; TacsScalar *r = NULL; if (res) { r = &res[offset]; @@ -1250,9 +1246,6 @@ class TACSQuadraticRotation { } } - q += vars_per_node; - qdot += vars_per_node; - dd += 3; t += 3; m += vars_per_node * size; @@ -2146,7 +2139,6 @@ class TACSQuaternionRotation { TacsScalar res[]) { TacsScalar *r = &res[offset]; const TacsScalar *q = &vars[offset]; - const TacsScalar *qdot = &dvars[offset]; for (int i = 0; i < num_nodes; i++) { // D = d(Qdot*t)/d(qdot) @@ -2173,7 +2165,6 @@ class TACSQuaternionRotation { r += vars_per_node; q += vars_per_node; - qdot += vars_per_node; dd += 3; t += 3; @@ -2264,8 +2255,6 @@ class TACSQuaternionRotation { Diddot[11] = 2.0 * (qiddot[1] * ti[0] + qiddot[2] * ti[1]); } - const TacsScalar *q = &vars[offset]; - const TacsScalar *qdot = &dvars[offset]; TacsScalar *r = NULL; if (res) { r = &res[offset]; @@ -2375,9 +2364,6 @@ class TACSQuaternionRotation { } } - q += vars_per_node; - qdot += vars_per_node; - dd += 3; t += 3; m += vars_per_node * size; diff --git a/src/elements/shell/TACSRBE3.cpp b/src/elements/shell/TACSRBE3.cpp index fba5ed0fe..afaa4e015 100644 --- a/src/elements/shell/TACSRBE3.cpp +++ b/src/elements/shell/TACSRBE3.cpp @@ -700,7 +700,7 @@ void TACSRBE3::addAdjResXptProduct( const TacsScalar Xpts[], const TacsScalar vars[], const TacsScalar dvars[], const TacsScalar ddvars[], TacsScalar fXptSens[]) { TacsScalar Xcg[3], sXcg[3], Jcg[3][3], sJcg[3][3], W[3], Lc, sLc; - const TacsScalar *X0, *Xn, *un, *F0, *M0, *u0, *t0, *tn; + const TacsScalar *X0, *Xn, *un, *F0, *M0, *tn; TacsScalar *sXpts, *sXn, *sX0; TacsScalar *maskedVars; @@ -918,8 +918,6 @@ void TACSRBE3::addAdjResXptProduct( transpose of the load transfer matrix above (using the Lagrange multiplier method) to figure out what the coefficients need to be. */ // start with dependent node (again, trivial) - u0 = &vars[0]; // dep node displacements - t0 = &vars[3]; // dep nodes rotations // displacements res[0] = 0.0; // u @@ -1329,7 +1327,7 @@ TacsScalar TACSRBE3::getMomentsOfInertiaSens( TacsScalar Icg[9], r[3], sr[3], temp[9], sIcg[9], stemp[9], c[6], Lc; Lc = 0.0; *sLc = 0.0; - TacsScalar detIcg, sdetIcg; + TacsScalar sdetIcg; for (int j = 0; j < 9; j++) { Icg[j] = sIcg[j] = 0.0; } @@ -1413,7 +1411,7 @@ TacsScalar TACSRBE3::getMomentsOfInertiaSens( } // Find sensitivity of inverse matrix - detIcg = inv3x3Sens(Icg, sIcg, temp, stemp, &sdetIcg); + inv3x3Sens(Icg, sIcg, temp, stemp, &sdetIcg); // Recast into 2D array for (int row = 0; row < 3; row++) { diff --git a/src/io/TACSMeshLoader.cpp b/src/io/TACSMeshLoader.cpp index 94ebca70c..7408a78cf 100644 --- a/src/io/TACSMeshLoader.cpp +++ b/src/io/TACSMeshLoader.cpp @@ -1078,7 +1078,7 @@ TACSToFH5 *TACSMeshLoader::createTACSToFH5(TACSAssembler *tacs, for (int k = 0; k < num_components; k++) { if (strlen(&component_descript[33 * k]) == 0) { char name[64]; - sprintf(name, "Component %d", k); + snprintf(name, sizeof(name), "Component %d", k); f5->setComponentName(k, name); } else { f5->setComponentName(k, &component_descript[33 * k]); diff --git a/src/io/TACSToFH5.cpp b/src/io/TACSToFH5.cpp index d350ba460..ddf1a8292 100644 --- a/src/io/TACSToFH5.cpp +++ b/src/io/TACSToFH5.cpp @@ -66,7 +66,7 @@ TACSToFH5::TACSToFH5(TACSAssembler *_assembler, ElementType _elem_type, for (int k = 0; k < num_components; k++) { char comp_name[128]; - sprintf(comp_name, "Component %d", k); + snprintf(comp_name, sizeof(comp_name), "Component %d", k); setComponentName(k, comp_name); } } @@ -154,7 +154,7 @@ int TACSToFH5::writeToFile(const char *filename) { int vars_per_node = assembler->getVarsPerNode(); // Find the maximum string length - int str_len = strlen("X,Y,Z") + 1; + size_t str_len = strlen("X,Y,Z") + 1; int nd = TacsGetOutputComponentCount(elem_type, TACS_OUTPUT_DISPLACEMENTS); int k = 0; for (; (k < nd && k < vars_per_node); k++) { @@ -164,7 +164,7 @@ int TACSToFH5::writeToFile(const char *filename) { } for (; k < vars_per_node; k++) { char stemp[64]; - sprintf(stemp, "v%d", k); + snprintf(stemp, sizeof(stemp), "v%d", k); str_len += strlen(stemp) + 1; } int nl = TacsGetOutputComponentCount(elem_type, TACS_OUTPUT_LOADS); @@ -176,17 +176,16 @@ int TACSToFH5::writeToFile(const char *filename) { } for (; k < vars_per_node; k++) { char stemp[64]; - sprintf(stemp, "f%d", k); + snprintf(stemp, sizeof(stemp), "f%d", k); str_len += strlen(stemp) + 1; } char *var_names = new char[str_len]; var_names[0] = '\0'; if (write_flag & TACS_OUTPUT_NODES) { - sprintf(var_names, "X,Y,Z"); + snprintf(var_names, str_len, "X,Y,Z"); } if (write_flag & TACS_OUTPUT_DISPLACEMENTS) { - str_len = strlen(var_names); nd = TacsGetOutputComponentCount(elem_type, TACS_OUTPUT_DISPLACEMENTS); k = 0; for (; (k < nd && k < vars_per_node); k++) { @@ -194,18 +193,17 @@ int TACSToFH5::writeToFile(const char *filename) { TacsGetOutputComponentName(elem_type, TACS_OUTPUT_DISPLACEMENTS, k); size_t len = strlen(var_names); if (k == 0 && !(write_flag & TACS_OUTPUT_NODES)) { - sprintf(&(var_names[len]), "%s", stemp); + snprintf(&(var_names[len]), str_len - len, "%s", stemp); } else { - sprintf(&(var_names[len]), ",%s", stemp); + snprintf(&(var_names[len]), str_len - len, ",%s", stemp); } } for (; k < vars_per_node; k++) { size_t len = strlen(var_names); - sprintf(&(var_names[len]), ",v%d", k); + snprintf(&(var_names[len]), str_len - len, ",v%d", k); } } if (write_flag & TACS_OUTPUT_LOADS) { - str_len = strlen(var_names); nl = TacsGetOutputComponentCount(elem_type, TACS_OUTPUT_LOADS); k = 0; for (; (k < nl && k < vars_per_node); k++) { @@ -214,14 +212,14 @@ int TACSToFH5::writeToFile(const char *filename) { size_t len = strlen(var_names); if (k == 0 && !(write_flag & TACS_OUTPUT_NODES || write_flag & TACS_OUTPUT_DISPLACEMENTS)) { - sprintf(&(var_names[len]), "%s", stemp); + snprintf(&(var_names[len]), str_len - len, "%s", stemp); } else { - sprintf(&(var_names[len]), ",%s", stemp); + snprintf(&(var_names[len]), str_len - len, ",%s", stemp); } } for (; k < vars_per_node; k++) { size_t len = strlen(var_names); - sprintf(&(var_names[len]), ",f%d", k); + snprintf(&(var_names[len]), str_len - len, ",f%d", k); } } @@ -324,7 +322,7 @@ int TACSToFH5::writeToFile(const char *filename) { // Write the data with a time stamp from the simulation in TACS char data_name[128]; double t = assembler->getSimulationTime(); - sprintf(data_name, "continuous data t=%.10e", t); + snprintf(data_name, sizeof(data_name), "continuous data t=%.10e", t); file->writeZoneData(data_name, var_names, TACSFH5File::FH5_FLOAT, dim1, dim2, float_data); delete[] float_data; @@ -351,7 +349,7 @@ int TACSToFH5::writeToFile(const char *filename) { // Write the data with a time stamp from the simulation in TACS char data_name[128]; double t = assembler->getSimulationTime(); - sprintf(data_name, "element data t=%.10e", t); + snprintf(data_name, sizeof(data_name), "element data t=%.10e", t); file->writeZoneData(data_name, variable_names, TACSFH5File::FH5_FLOAT, dim1, dim2, float_data); delete[] float_data; @@ -518,7 +516,7 @@ char *TACSToFH5::getElementVarNames(int flag) { for (int i = 1; i < nd; i++) { stemp = TacsGetOutputComponentName(elem_type, out_types[k], i); size_t len = strlen(temp); - sprintf(&(temp[len]), ",%s", stemp); + snprintf(&(temp[len]), str_len - len, ",%s", stemp); } } output_names[k] = temp; @@ -550,7 +548,7 @@ char *TACSToFH5::getElementVarNames(int flag) { for (; k < 3; k++) { if (output_names[k]) { int len = strlen(elem_vars); - sprintf(&elem_vars[len], ",%s", output_names[k]); + snprintf(&elem_vars[len], elem_size - len, ",%s", output_names[k]); } } diff --git a/tacs/TACS.pxd b/tacs/TACS.pxd index f3b25ffbb..2c291263a 100644 --- a/tacs/TACS.pxd +++ b/tacs/TACS.pxd @@ -12,18 +12,10 @@ # distutils: language=c++ -# For MPI capabilities -from mpi4py.libmpi cimport * -cimport mpi4py.MPI as MPI - # Import numpy from libc.string cimport const_char from libcpp cimport bool -# Import numpy -cimport numpy as np -import numpy as np - # Import TACS c++ headers from tacs.cpp_headers.TACS cimport * diff --git a/tacs/TACS.pyx b/tacs/TACS.pyx index 358b289ec..85543578f 100644 --- a/tacs/TACS.pyx +++ b/tacs/TACS.pyx @@ -28,7 +28,6 @@ from libc.string cimport const_char from libc.stdlib cimport malloc, free from libcpp cimport bool - # Import C methods for python from cpython cimport PyObject, Py_INCREF @@ -3262,6 +3261,37 @@ cdef class FrequencyAnalysis: """ self.ptr.evalEigenXptSens(index, xptsens.getBVecPtr()) + def addEigenSens(self, int neigs, dfdlam, dfdq, dfdx=None, dfdXpts=None, + int use_cg=0, rtol=1e-8, atol=1e-30): + """ + Add the derivative of the function wrt the eigenvalues and eigenvectors + """ + + cdef TacsScalar *dfdlam_ = NULL + cdef TACSBVec **dfdq_ = NULL + cdef TACSBVec *dfdx_ = NULL + cdef TACSBVec *dfdXpts_ = NULL + + dfdlam_ = malloc(neigs * sizeof(TacsScalar)) + for i in range(neigs): + dfdlam_[i] = dfdlam[i] + + dfdq_ = malloc(neigs * sizeof(TACSBVec*)) + for i in range(neigs): + dfdq_[i] = (dfdq[i]).getBVecPtr() + + if dfdx is not None: + dfdx_ = (dfdx).getBVecPtr() + if dfdXpts is not None: + dfdXpts_ = (dfdXpts).getBVecPtr() + + self.ptr.addEigenSens(neigs, dfdlam_, dfdq_, dfdx_, dfdXpts_, use_cg, rtol, atol) + + free(dfdlam_) + free(dfdq_) + + return + cdef class BucklingAnalysis: cdef TACSLinearBuckling *ptr def __cinit__(self, Assembler assembler, TacsScalar sigma, diff --git a/tacs/cpp_headers/TACS.pxd b/tacs/cpp_headers/TACS.pxd index 220a819c8..d3f2bbba5 100644 --- a/tacs/cpp_headers/TACS.pxd +++ b/tacs/cpp_headers/TACS.pxd @@ -495,6 +495,9 @@ cdef extern from "TACSBuckling.h": TacsScalar extractEigenvector(int, TACSBVec*, TacsScalar*) void evalEigenDVSens(int, TACSBVec*) void evalEigenXptSens(int, TACSBVec*) + void addEigenSens(int neigs, TacsScalar dfdlam[], + TACSBVec *dfdq[], TACSBVec *dfdx, + TACSBVec *dfdXpt, int, double, double) cdef cppclass TACSLinearBuckling(TACSObject): TACSLinearBuckling( TACSAssembler *, diff --git a/tacs/elements.pxd b/tacs/elements.pxd index e33caea17..e40471167 100644 --- a/tacs/elements.pxd +++ b/tacs/elements.pxd @@ -12,10 +12,6 @@ #distutils: language=c++ -# Import numpy -cimport numpy as np -import numpy as np - # Import from constitutive for definitions from tacs.cpp_headers.constitutive cimport * from tacs.cpp_headers.TACS cimport *