[cig-commits] r16812 - short/3D/PyLith/trunk/libsrc/faults
brad at geodynamics.org
brad at geodynamics.org
Thu May 27 13:45:43 PDT 2010
Author: brad
Date: 2010-05-27 13:45:43 -0700 (Thu, 27 May 2010)
New Revision: 16812
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
Log:
Fixed bug in full cell matrix fault preconditioner implementation (but left use of full precondigioner disabled). Compute -(CAiC^T)^(-1).
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-05-27 20:22:46 UTC (rev 16811)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-05-27 20:45:43 UTC (rev 16812)
@@ -706,7 +706,7 @@
jacobianInvVertexP[iDim] = 1.0/jacobianVertexP[iDim*spaceDim+iDim];
} // for
- // Compute [C] [Adiag]^(-1) [C]^T
+ // Compute -[C] [Adiag]^(-1) [C]^T
// C_{ij} = orientationVertex[i*spaceDim+j]
// C^T_{ij} = orientationVertex[j*spaceDim+i]
// Adiag^{-1}_{ii} = jacobianInvVertexN[i] + jacobianInvVertexP[i] (BRAD: Are you sure its not a minus sign here?)
@@ -715,12 +715,12 @@
for (int kDim=0; kDim < spaceDim; ++kDim) {
double value = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim)
- value +=
+ value -=
orientationVertex[kDim*spaceDim+iDim] *
orientationVertex[kDim*spaceDim+iDim] *
(jacobianInvVertexN[iDim] + jacobianInvVertexP[iDim]);
if (fabs(value) > 1.0e-8)
- precondVertexL[kDim] = -1.0 / value;
+ precondVertexL[kDim] = 1.0 / value;
} // for
@@ -758,7 +758,7 @@
#else // FULL PRECONDITIONER
- // Compute [C] [A]^(-1) [C]^T for cell.
+ // Compute -( [C] [A]^(-1) [C]^T )^-1 for cell.
//
// numBasis = number of corners in fault cell
// spaceDim = spatial dimension
@@ -800,6 +800,7 @@
// Allocate vectors for vertex values
double_array preconditionerCell(matrixSizeF);
+ double_array preconditionerInvCell(matrixSizeF);
int_array indicesN(nrowsF);
int_array indicesP(nrowsF);
int_array indicesLagrange(nrowsF);
@@ -900,6 +901,7 @@
jacobianCellP = 0.0;
jacobianCellN = 0.0;
+ preconditionerInvCell = 0.0;
preconditionerCell = 0.0;
#if defined(DETAILED_EVENT_LOGGING)
@@ -1051,7 +1053,7 @@
for (int lLagrange=0; lLagrange < numBasis; ++lLagrange) {
// Exploit structure of C in matrix multiplication.
- // C_ij Ai_jk C_lk - Structure of C means k == l;
+ // -C_ij Ai_jk C_lk - Structure of C means k == l;
const int kBasis = lLagrange;
for (int iDim=0; iDim < spaceDim; ++iDim) {
@@ -1066,7 +1068,7 @@
for (int kDim=0; kDim < spaceDim; ++kDim) {
const int kB = kBasis*spaceDim + kDim;
- preconditionerCell[iL*nrowsF+lL] +=
+ preconditionerInvCell[iL*nrowsF+lL] -=
orientationCell[iLagrange*orientationSize+iDim*spaceDim+jDim] *
jacobianInvCellN[jB*nrowsF+kB] *
orientationCell[lLagrange*orientationSize+kDim*spaceDim+lDim];
@@ -1078,30 +1080,88 @@
} // for
} // for
+#if 1
+ std::cout << "1/P_cell " << *c_iter << std::endl;
+ for(int i = 0; i < nrowsF; ++i) {
+ for(int j = 0; j < nrowsF; ++j) {
+ std::cout << " " << preconditionerInvCell[i*nrowsF+j];
+ }
+ std::cout << std::endl;
+ }
+#endif
+
+ // Invert (C Ai C^T)
+
+ // Transpose matrices so we can call LAPACK
+ for(int i = 0; i < nrowsF; ++i) {
+ for(int j = 0; j < i; ++j) {
+ PetscInt k = i*nrowsF+j;
+ PetscInt kp = j*nrowsF+i;
+ PetscScalar tmp;
+ tmp = preconditionerInvCell[k];
+ preconditionerInvCell[k] = preconditionerInvCell[kp];
+ preconditionerInvCell[kp] = tmp;
+ }
+ }
+ LAPACKgesvd_("A", "A", &elemRows, &elemRows, &preconditionerInvCell[0], &elemRows, &singularValuesN[0], &UN[0], &elemRows, &VNt[0], &elemRows, &work[0], &workSize, &berr);
+ CHECK_PETSC_ERROR_MSG(berr, "Inversion of preconditioner element matrix failed.");
+
+#if 1
+ for(int i = 0; i < nrowsF; ++i) {
+ std::cout << "sigmaN["<<i<<"]: " << singularValuesN[i] << " sigmaP["<<i<<"]: " << singularValuesP[i] << std::endl;
+ }
+ std::cout << "UN_cell " << *c_iter << std::endl;
+ for(int i = 0; i < nrowsF; ++i) {
+ for(int j = 0; j < nrowsF; ++j) {
+ std::cout << " " << UN[j*nrowsF+i];
+ }
+ std::cout << std::endl;
+ }
+ std::cout << "VNt_cell " << *c_iter << std::endl;
+ for(int i = 0; i < nrowsF; ++i) {
+ for(int j = 0; j < nrowsF; ++j) {
+ std::cout << " " << VNt[j*nrowsF+i];
+ }
+ std::cout << std::endl;
+ }
+#endif
+
+ // Row scale Vt by the inverse of the singular values
+ for(int i = 0; i < nrowsF; ++i) {
+ const PetscReal invN = singularValuesN[i] > 1.0e-10 ? 1.0/singularValuesN[i] : 0.0;
+ const PetscReal invP = singularValuesP[i] > 1.0e-10 ? 1.0/singularValuesP[i] : 0.0;
+
+ for(int j = 0; j < nrowsF; ++j) {
+ VNt[j*nrowsF+i] *= invN;
+ VPt[j*nrowsF+i] *= invP;
+ }
+ }
+ BLASgemm_("N", "N", &elemRows, &elemRows, &elemRows, &one, &UN[0], &elemRows, &VNt[0], &elemRows, &one, &preconditionerCell[0], &elemRows);
+
+ // Transpose matrices from LAPACK
+ for(int i = 0; i < nrowsF; ++i) {
+ for(int j = 0; j < i; ++j) {
+ PetscInt k = i*nrowsF+j;
+ PetscInt kp = j*nrowsF+i;
+ PetscScalar tmp;
+ tmp = preconditionerCell[k];
+ preconditionerCell[k] = preconditionerCell[kp];
+ preconditionerCell[kp] = tmp;
+ }
+ }
+
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
_logger->eventBegin(updateEvent);
#endif
-#if 0
err = MatSetValues(*precondMatrix,
indicesLagrange.size(), &indicesLagrange[0],
indicesLagrange.size(), &indicesLagrange[0],
&preconditionerCell[0],
ADD_VALUES);
CHECK_PETSC_ERROR_MSG(err, "Setting values in fault preconditioner failed.");
-#else
- for (int iLagrange=0; iLagrange < numBasis; ++iLagrange)
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- const int iL = iLagrange*spaceDim + iDim;
- const int indexL = indicesLagrange[iL];
- MatSetValue(*precondMatrix,
- indexL, indexL,
- preconditionerCell[iL*nrowsF+iL],
- ADD_VALUES);
- } // for
-
-#endif
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(updateEvent);
More information about the CIG-COMMITS
mailing list