#ifdef PETSC_RCS_HEADER static char vcid[] = "$Id: gmat.c,v 1.34 2000/07/16 23:20:03 knepley Exp $"; #endif /* This file provides routines for grid matrices */ #include "src/gvec/gvecimpl.h" /*I "gvec.h" I*/ #include "petscts.h" #include "gsolver.h" /* Logging support */ int GMAT_CreateRectangular, GMAT_EvaluateOperatorGalerkin, GMAT_EvaluateSystemMatrix, GMAT_SetBoundary; int GMAT_MatMultConstrained, GMAT_MatMultTransposeConstrained; extern int MatMult_MPIAIJ(Mat, Vec, Vec); /* For GMatMatMultConstrained */ #undef __FUNCT__ #define __FUNCT__ "GMatDestroy" /*@C GMatDestroy - Destroys a grid matrix. Input Parameter: . mat - the matrix Level: beginner .keywords: matrix, destroy @*/ int GMatDestroy(Mat mat) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); if (--mat->refct > 0) PetscFunctionReturn(0); ierr = MatDestroy(mat); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatView" /*@ GMatView - Views a grid matrix. Input Parameters: + mat - the grid matrix - viewer - an optional visualization context Notes: GMatView() supports the same viewers as MatView(). The only difference is that for all multiprocessor cases, the output vector employs the natural ordering of the grid, so it many cases this corresponds to the ordering that would have been used for the uniprocessor case. The available visualization contexts include $ VIEWER_STDOUT_SELF - standard output (default) $ VIEWER_STDOUT_WORLD - synchronized standard $ output where only the first processor opens $ the file. All other processors send their $ data to the first processor to print. The user can open alternative visualization contexts with $ PetscViewerFileOpenASCII() - output vector to a specified file $ PetscViewerFileOpenBinary() - output in binary to a $ specified file; corresponding input uses VecLoad() $ PetscViewerDrawOpenX() - output vector to an X window display $ DrawLGCreate() - output vector as a line graph to an X window display $ PetscViewerMatlabOpen() - output vector to Matlab viewer Level: beginner .keywords: view, visualize, output, print, write, draw .seealso: MatView() @*/ int GMatView(GMat mat, PetscViewer viewer) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE); if (!viewer) { viewer = PETSC_VIEWER_STDOUT_SELF; } else { PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE); } ierr = GMatGetGrid(mat, &grid); CHKERRQ(ierr); ierr = (*grid->ops->gmatview)(mat, viewer); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatSerialize" /*@ GMatSerialize - This function stores or recreates a grid matrix using a viewer for a binary file. Input Parameters: . viewer - The viewer context . store - This flag is PETSC_TRUE is data is being written, otherwise it will be read Output Parameter: . m - The grid matrix Level: beginner .keywords: grid vector, serialize .seealso: GridSerialize() @*/ int GMatSerialize(Grid grid, GMat *m, PetscViewer viewer, PetscTruth store) { int fd; GMat mat; MatInfo info; PetscScalar *vals, *vals2; int *cols, *cols2; int *diag; int *offdiag; int *firstCol; int type, rowVars, rowLocVars, colVars, colLocVars, numNonZeros; int rowStart, rowEnd, colStart, colEnd, offset, size; int numProcs, rank; int proc, row, col; PetscTruth match; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidHeaderSpecific(viewer, PETSC_VIEWER_COOKIE); PetscValidPointer(m); ierr = PetscTypeCompare((PetscObject) viewer, PETSC_VIEWER_BINARY, &match); CHKERRQ(ierr); if (match == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_WRONG, "Must be binary viewer"); ierr = PetscViewerBinaryGetDescriptor(viewer, &fd); CHKERRQ(ierr); if (store) { PetscValidHeaderSpecific(*m, MAT_COOKIE); ierr = MatGetSize(*m, &rowVars, &colVars); CHKERRQ(ierr); ierr = MatGetLocalSize(*m, &rowLocVars, &colLocVars); CHKERRQ(ierr); ierr = MatGetInfo(*m, MAT_LOCAL, &info); CHKERRQ(ierr); ierr = MatGetOwnershipRange(*m, &rowStart, &rowEnd); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &(*m)->cookie, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &rowVars, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &rowLocVars, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &colVars, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &colLocVars, 1, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, &info.nz_used, 1, PETSC_INT, 0); CHKERRQ(ierr); numNonZeros = (int) info.nz_used; ierr = PetscMalloc((rowLocVars*2 + numNonZeros) * sizeof(int) + numNonZeros * sizeof(PetscScalar), &diag); CHKERRQ(ierr); offdiag = diag + rowLocVars; cols = offdiag + rowLocVars; vals = (PetscScalar *) (cols + numNonZeros); ierr = PetscMemzero(diag, rowLocVars*2 * sizeof(int)); CHKERRQ(ierr); ierr = MPI_Comm_size(grid->comm, &numProcs); CHKERRQ(ierr); ierr = MPI_Comm_rank(grid->comm, &rank); CHKERRQ(ierr); ierr = PetscMalloc((numProcs+1) * sizeof(int), &firstCol); CHKERRQ(ierr); MPI_Allgather(&colLocVars, 1, MPI_INT, &firstCol[1], 1, MPI_INT, grid->comm); for(proc = 1, firstCol[0] = 0; proc <= numProcs; proc++) firstCol[proc] += firstCol[proc-1]; if (firstCol[numProcs] != colVars) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid column partition"); colStart = firstCol[rank]; colEnd = firstCol[rank+1]; for(row = rowStart, offset = 0; row < rowEnd; row++, offset += size) { ierr = MatGetRow(*m, row, &size, &cols2, &vals2); CHKERRQ(ierr); for(col = 0; col < size; col++) if ((col >= colStart) && (col < colEnd)) { diag[row]++; } else { offdiag[row]++; } ierr = PetscMemcpy(&cols[offset], cols2, size * sizeof(int)); CHKERRQ(ierr); ierr = PetscMemcpy(&vals[offset], vals2, size * sizeof(PetscScalar)); CHKERRQ(ierr); ierr = MatRestoreRow(*m, row, &size, &cols2, &vals2); CHKERRQ(ierr); } ierr = PetscBinaryWrite(fd, diag, rowLocVars, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, offdiag, rowLocVars, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, cols, numNonZeros, PETSC_INT, 0); CHKERRQ(ierr); ierr = PetscBinaryWrite(fd, vals, numNonZeros, PETSC_SCALAR, 0); CHKERRQ(ierr); ierr = PetscFree(diag); CHKERRQ(ierr); ierr = PetscFree(offdiag); CHKERRQ(ierr); ierr = PetscFree(cols); CHKERRQ(ierr); ierr = PetscFree(vals); CHKERRQ(ierr); } else { ierr = PetscBinaryRead(fd, &type, 1, PETSC_INT); CHKERRQ(ierr); if (type != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_ARG_WRONG, "Non-matrix object"); ierr = PetscBinaryRead(fd, &rowVars, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &rowLocVars, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &colVars, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &colLocVars, 1, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, &numNonZeros, 1, PETSC_INT); CHKERRQ(ierr); MPI_Reduce(&rowLocVars, &size, 1, MPI_INT, MPI_SUM, 0, grid->comm); if (size != rowVars) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid row partition"); MPI_Reduce(&colLocVars, &size, 1, MPI_INT, MPI_SUM, 0, grid->comm); if (size != colVars) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid column partition"); ierr = PetscMalloc((rowLocVars*2 + numNonZeros) * sizeof(int) + numNonZeros * sizeof(PetscScalar), &diag); CHKERRQ(ierr); offdiag = diag + rowLocVars; cols = offdiag + rowLocVars; vals = (PetscScalar *) (cols + numNonZeros); ierr = PetscBinaryRead(fd, diag, rowLocVars, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, offdiag, rowLocVars, PETSC_INT); CHKERRQ(ierr); ierr = MatCreateMPIAIJ(grid->comm, rowLocVars, colLocVars, rowVars, colVars, 0, diag, 0, offdiag, &mat); CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) mat, "Grid", (PetscObject) grid); CHKERRQ(ierr); ierr = MatGetOwnershipRange(mat, &rowStart, &rowEnd); CHKERRQ(ierr); if (rowEnd - rowStart + 1 != rowLocVars) SETERRQ(PETSC_ERR_ARG_CORRUPT, "Invalid row partition"); ierr = PetscBinaryRead(fd, cols, numNonZeros, PETSC_INT); CHKERRQ(ierr); ierr = PetscBinaryRead(fd, vals, numNonZeros, PETSC_SCALAR); CHKERRQ(ierr); for(row = rowStart, offset = 0; row < rowEnd; row++, offset += size) { size = diag[row] + offdiag[row]; ierr = MatSetValues(mat, 1, &row, size, &cols[offset], &vals[offset], INSERT_VALUES); CHKERRQ(ierr); } ierr = PetscFree(diag); CHKERRQ(ierr); ierr = PetscFree(offdiag); CHKERRQ(ierr); ierr = PetscFree(cols); CHKERRQ(ierr); ierr = PetscFree(vals); CHKERRQ(ierr); ierr = MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); ierr = MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); *m = mat; } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatDuplicate" /*@C GMatDuplicate - Duplicates a grid vector. Input Parameters: . mat - The matrix Level: beginner .keywords: grid matrix, destroy .seealso: MatDuplicate(), GMatDestroy() @*/ int GMatDuplicate(GMat mat, GMat *newmat) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); PetscValidPointer(newmat); ierr = MatConvert(mat, MATSAME, newmat); CHKERRQ(ierr); PetscFunctionReturn(ierr); } #undef __FUNCT__ #define __FUNCT__ "GMatGetSize" /*@ GMatGetSize - Returns the numbers of rows and columns in a matrix. Not Collective Input Parameter: . mat - the matrix Output Parameters: + M - the number of global rows - N - the number of global columns Level: intermediate .keywords: matrix, dimension, size, rows, columns, global, get .seealso: GMatGetLocalSize() @*/ int GMatGetSize(GMat mat, int *M, int* N) { Grid grid; PetscTruth isConstrained = PETSC_FALSE; PetscTruth explicitConstraints = PETSC_TRUE; int gM, gN; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); ierr = GMatGetGrid(mat, &grid); ierr = GridIsConstrained(grid, &isConstrained); CHKERRQ(ierr); ierr = GridGetExplicitConstraints(grid, &explicitConstraints); CHKERRQ(ierr); if ((isConstrained == PETSC_FALSE) || (explicitConstraints == PETSC_TRUE)) { if (M) *M = mat->M; if (N) *N = mat->N; } else { ierr = GridGetConstraints(grid, PETSC_NULL, PETSC_NULL, &gM, PETSC_NULL); CHKERRQ(ierr); gN = gM; /* KLUDGE - Must catch matrices which arise from nonstandard orderings if (mat->N == x->N) gN = x->N; if (mat->M == y->N) gM = y->N; */ if (M) *M = gM; if (N) *N = gN; } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatGetLocalSize" /*@ GMatGetLocalSize - Returns the number of rows and columns in a matrix stored locally. This information may be implementation dependent, so use with care. Not Collective Input Parameters: . mat - the matrix Output Parameters: + m - the number of local rows - n - the number of local columns Level: intermediate .keywords: matrix, dimension, size, local, rows, columns, get .seealso: GMatGetSize() @*/ int GMatGetLocalSize(GMat mat, int *m, int* n) { Grid grid; PetscTruth isConstrained = PETSC_FALSE; PetscTruth explicitConstraints = PETSC_TRUE; int gm, gn; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); ierr = GMatGetGrid(mat, &grid); ierr = GridIsConstrained(grid, &isConstrained); CHKERRQ(ierr); ierr = GridGetExplicitConstraints(grid, &explicitConstraints); CHKERRQ(ierr); if ((isConstrained == PETSC_FALSE) || (explicitConstraints == PETSC_TRUE)) { if (m) *m = mat->m; if (n) *n = mat->n; } else { ierr = GridGetConstraints(grid, PETSC_NULL, PETSC_NULL, PETSC_NULL, &gm); CHKERRQ(ierr); gn = gm; /* KLUDGE - Must catch matrices which arise from nonstandard orderings if (mat->n == x->n) gn = x->n; if (mat->m == y->n) gm = y->n; */ if (m) *m = gm; if (n) *n = gn; } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatGetGrid" /*@ GMatGetGrid - This function returns the grid from a grid matrix. Not collective Input Parameter: . m - The grid matrix Output Parameter: . grid - The grid Level: intermediate .keywords: grid matrix, grid, get .seealso: GridGetMesh(), GVecGetGrid() @*/ int GMatGetGrid(GMat m, Grid *grid) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(m, MAT_COOKIE); PetscValidPointer(grid); ierr = PetscObjectQuery((PetscObject) m, "Grid", (PetscObject *) grid); CHKERRQ(ierr); PetscValidHeaderSpecific(*grid, GRID_COOKIE); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatGetOrder" /*@ GMatGetOrder - This function returns the orderings from a grid matrix. Not collective Input Parameter: . m - The grid matrix Output Parameters: + rowOrder - The row (or test function) ordering - colOrder - The column (or shape function) ordering Level: intermediate .keywords: grid matrix, variable ordering, get .seealso: GMatGetGrid(), GVecGetOrder() @*/ int GMatGetOrder(GMat m, VarOrdering *rowOrder, VarOrdering *colOrder) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(m, MAT_COOKIE); if (rowOrder != PETSC_NULL) { PetscValidPointer(rowOrder); ierr = PetscObjectQuery((PetscObject) m, "RowOrder", (PetscObject *) rowOrder); CHKERRQ(ierr); PetscValidHeaderSpecific(*rowOrder, VAR_ORDER_COOKIE); } if (colOrder != PETSC_NULL) { PetscValidPointer(colOrder); ierr = PetscObjectQuery((PetscObject) m, "ColOrder", (PetscObject *) colOrder); CHKERRQ(ierr); PetscValidHeaderSpecific(*colOrder, VAR_ORDER_COOKIE); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatGetDiagonalConstrained" /*@ GMatGetDiagonalConstrained - This function returns the diagonal of a constrained matrix. Input Paramter: . mat - The grid matrix Output Paramter: . diag - A constrained grid vector containing the diagonal elements Level: advanced .keyowrds grid matrix, constraint .seealso MatGetDiagonal() @*/ int GMatGetDiagonalConstrained(GMat mat, GVec diag) { Grid grid; PetscTruth isConstrained, explicitConstraints; Vec diagLong; PetscScalar *array; PetscScalar *arrayC; int *ordering; int numInteriorVars; int var, size; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); PetscValidHeaderSpecific(diag, VEC_COOKIE); ierr = GMatGetGrid(mat, &grid); CHKERRQ(ierr); ierr = GridIsConstrained(grid, &isConstrained); CHKERRQ(ierr); ierr = GridGetExplicitConstraints(grid, &explicitConstraints); CHKERRQ(ierr); if ((isConstrained == PETSC_FALSE) || (explicitConstraints == PETSC_TRUE)) { ierr = MatGetDiagonal(mat, diag); CHKERRQ(ierr); PetscFunctionReturn(0); } ierr = VecGetLocalSize(diag, &size); CHKERRQ(ierr); if (size != grid->constraintOrder->numLocVars) SETERRQ(PETSC_ERR_ARG_INCOMP, "Vector wrong size for matrix"); ierr = ISGetSize(grid->constraintOrdering, &size); CHKERRQ(ierr); if (size != grid->order->numLocVars) SETERRQ(PETSC_ERR_ARG_INCOMP, "Constraint mapping wrong size for matrix"); /* First copy all the interior variables */ ierr = GVecCreate(grid, &diagLong); CHKERRQ(ierr); ierr = MatGetDiagonal(mat, diagLong); CHKERRQ(ierr); ierr = VecGetArray(diagLong, &array); CHKERRQ(ierr); ierr = VecGetArray(diag, &arrayC); CHKERRQ(ierr); ierr = ISGetIndices(grid->constraintOrdering, &ordering); CHKERRQ(ierr); numInteriorVars = grid->constraintOrder->numLocVars - grid->constraintOrder->numLocNewVars; for(var = 0; var < grid->order->numLocVars; var++) { if (ordering[var] < numInteriorVars) { if (PetscAbsScalar(array[var]) > PETSC_MACHINE_EPSILON) arrayC[ordering[var]] = array[var]; else arrayC[ordering[var]] = 1.0; } } ierr = ISRestoreIndices(grid->constraintOrdering, &ordering); CHKERRQ(ierr); ierr = VecRestoreArray(diagLong, &array); CHKERRQ(ierr); ierr = VecRestoreArray(diag, &arrayC); CHKERRQ(ierr); ierr = GVecDestroy(diagLong); CHKERRQ(ierr); /* Now ask grid for diagonal of extra variables */ if ((isConstrained == PETSC_TRUE) && (grid->constraintCtx->ops->jacgetdiag)) { ierr = (*grid->constraintCtx->ops->jacgetdiag)(mat, diag); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatGetDiagonalMF" /*@ GMatGetDiagonalMF - This function returns the diagonal of the matrix. Input Parameter: . mat - The grid matrix Output Parameter: . diag - A grid vector containing the diagonal elements Notes: See comments in GMatGetDiagonalConstrained() about dealing with implicit constraints Level: intermediate .keywords grid matrix, constraint .seealso GMatGetDiagonalConstrained @*/ int GMatGetDiagonalMF(GMat mat, GVec diag) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); PetscValidHeaderSpecific(diag, VEC_COOKIE); ierr = GMatGetGrid(mat, &grid); CHKERRQ(ierr); ierr = GVecEvaluateJacobianDiagonal(diag, grid->matrixFreeArg, grid->matrixFreeContext); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatDiagonalScaleConstrained" /*@ GMatDiagonalScaleConstrained - Scales a constrained matrix on the left and right by diagonal matrices that are stored as vectors. Either of the two scaling matrices can be PETSC_NULL. Input Parameters: . mat - The grid matrix to be scaled . l - The left scaling vector (or PETSC_NULL) . r - The right scaling vector (or PETSC_NULL) Notes: GMatDiagonalScaleConstrained() computes A <- LAR, where $ L = a diagonal matrix $ R = a diagonal matrix Level: advanced .keywords: matrix, diagonal, scale .seealso: MatDiagonalScale() @*/ int GMatDiagonalScaleConstrained(GMat mat, GVec l, GVec r) { Grid grid; GVec newL, newR; PetscScalar *arrayL, *arrayR; PetscScalar *arrayNewL, *arrayNewR; PetscScalar one = 1.0; int *ordering; int numInteriorVars; int var, size; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); ierr = GMatGetGrid(mat, &grid); CHKERRQ(ierr); PetscValidHeaderSpecific(grid, GRID_COOKIE); if (grid->isConstrained == PETSC_FALSE) { ierr = MatDiagonalScale(mat, l, r); CHKERRQ(ierr); PetscFunctionReturn(0); } ierr = ISGetSize(grid->constraintOrdering, &size); CHKERRQ(ierr); if (size != grid->order->numLocVars) SETERRQ(PETSC_ERR_ARG_INCOMP, "Constraint mapping wrong size for matrix"); newL = PETSC_NULL; newR = PETSC_NULL; arrayL = PETSC_NULL; arrayR = PETSC_NULL; numInteriorVars = grid->constraintOrder->numLocVars - grid->constraintOrder->numLocNewVars; ierr = ISGetIndices(grid->constraintOrdering, &ordering); CHKERRQ(ierr); if (l != PETSC_NULL) { PetscValidHeaderSpecific(l, VEC_COOKIE); ierr = VecGetLocalSize(l, &size); CHKERRQ(ierr); if (size != grid->constraintOrder->numLocVars) SETERRQ(PETSC_ERR_ARG_INCOMP, "Left vector wrong size for matrix"); ierr = GVecCreate(grid, &newL); CHKERRQ(ierr); ierr = VecSet(&one, newL); CHKERRQ(ierr); ierr = VecGetArray(l, &arrayL); CHKERRQ(ierr); ierr = VecGetArray(newL, &arrayNewL); CHKERRQ(ierr); for(var = 0; var < grid->order->numLocVars; var++) { if (ordering[var] < numInteriorVars) arrayNewL[var] = arrayL[ordering[var]]; } ierr = VecRestoreArray(l, &arrayL); CHKERRQ(ierr); ierr = VecRestoreArray(newL, &arrayNewL); CHKERRQ(ierr); } if (r != PETSC_NULL) { PetscValidHeaderSpecific(r, VEC_COOKIE); ierr = VecGetLocalSize(r, &size); CHKERRQ(ierr); if (size != grid->constraintOrder->numLocVars) SETERRQ(PETSC_ERR_ARG_INCOMP, "Right vector wrong size for matrix"); ierr = GVecCreate(grid, &newR); CHKERRQ(ierr); ierr = VecSet(&one, newR); CHKERRQ(ierr); ierr = VecGetArray(r, &arrayR); CHKERRQ(ierr); ierr = VecGetArray(newR, &arrayNewR); CHKERRQ(ierr); for(var = 0; var < grid->order->numLocVars; var++) { if (ordering[var] < numInteriorVars) arrayNewR[var] = arrayR[ordering[var]]; } ierr = VecRestoreArray(r, &arrayR); CHKERRQ(ierr); ierr = VecRestoreArray(newR, &arrayNewR); CHKERRQ(ierr); } ierr = ISRestoreIndices(grid->constraintOrdering, &ordering); CHKERRQ(ierr); ierr = MatDiagonalScale(mat, newL, newR); CHKERRQ(ierr); if (l != PETSC_NULL) { ierr = VecDestroy(newL); CHKERRQ(ierr); } if (r != PETSC_NULL) { ierr = VecDestroy(newR); CHKERRQ(ierr); } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatOrderConstrained" /*@ GMatOrderConstrained - This function creates an ordering which places all interior variables before variables eliminated by constraints. Input Paramter: . mat - The grid matrix . type - The reordering type Output Paramter: . rowIS - The row reordering . colIS - The column reordering Level: advanced .keyowrds grid matrix, constraint .seealso GMatGetDiagonalConstrained() @*/ int GMatOrderConstrained(GMat mat, MatOrderingType type, IS *rowIS, IS *colIS) { Grid grid; PetscTruth opt, match; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); PetscValidPointer(rowIS); PetscValidPointer(colIS); ierr = PetscStrcmp(type, MATORDERING_CONSTRAINED, &match); CHKERRQ(ierr); if (match == PETSC_FALSE) SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid ordering type %s", type); ierr = GMatGetGrid(mat, &grid); CHKERRQ(ierr); if (grid->constraintOrdering) { ierr = ISDuplicate(grid->constraintOrdering, rowIS); CHKERRQ(ierr); ierr = ISDuplicate(grid->constraintOrdering, colIS); CHKERRQ(ierr); ierr = PetscOptionsHasName(PETSC_NULL, "-gmat_order_nonzero_diagonal", &opt); CHKERRQ(ierr); if (opt == PETSC_TRUE) { /* Remove the zeros from the diagonal in both blocks */ ierr = GMatReorderForNonzeroDiagonalConstrained(mat, 1e-10, *rowIS, *colIS); CHKERRQ(ierr); } } else { *rowIS = PETSC_NULL; *colIS = PETSC_NULL; } PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatFindPreviousPivot_Private" /* GMatFindPreviousPivot_Private - Finds a pivot element in the lower triangle which does not introduce a zero on the above diagonal. Input Parameters: . mat - The grid matrix . prow - The row with a zero pivot . rowPerm - The row permutation . colPerm - The column permutation . atol - The smallest acceptable pivot Output Paramters: . replCol - The column to interchange with . replVal - The magnitude of the new pivot Level: developer */ int GMatFindPreviousPivot_Private(GMat mat, int prow, int *rowPerm, int *colPerm, double atol, int *replCol, double *replVal) { int nz; int *cols; PetscScalar *vals; int newNz; int *newCols; PetscScalar *newVals; int newRow; int repl; int colIndex, newColIndex; int ierr; PetscFunctionBegin; newRow = rowPerm[prow]; ierr = MatGetRow(mat, prow, &nz, &cols, &vals); CHKERRQ(ierr); for(colIndex = 0; colIndex < nz; colIndex++) { /* Find an acceptable pivot in the lower triangle */ if ((colPerm[cols[colIndex]] < newRow) && (PetscAbsScalar(vals[colIndex]) > atol)) { /* Check previous diagonal for zero pivot */ repl = cols[colIndex]; ierr = MatGetRow(mat, repl, &newNz, &newCols, &newVals); CHKERRQ(ierr); for(newColIndex = 0; newColIndex < newNz; newColIndex++) { if ((colPerm[newCols[newColIndex]] == newRow) && (PetscAbsScalar(newVals[newColIndex]) > atol)) { *replCol = repl; *replVal = PetscAbsScalar(vals[colIndex]); ierr = MatRestoreRow(mat, repl, &newNz, &newCols, &newVals); CHKERRQ(ierr); ierr = MatRestoreRow(mat, prow, &nz, &cols, &vals); CHKERRQ(ierr); PetscFunctionReturn(0); } } ierr = MatRestoreRow(mat, repl, &newNz, &newCols, &newVals); CHKERRQ(ierr); } } ierr = MatRestoreRow(mat, prow, &nz, &cols, &vals); CHKERRQ(ierr); /* Signal error since no acceptable pivot was found */ PetscFunctionReturn(1); } #undef __FUNCT__ #define __FUNCT__ "GMatReorderForNonzeroDiagonal_Private" /* GMatReorderForNonzeroDiagonal_Private - This is identical to MatReorderForNonzeroDiagonal(), but allows reordering of a block of the matrix. Input Parameters: . mat - The grid matrix to reorder . rowStart - The first row in the block . rowEnd - The last row in the block . rowIS - The row permutation, usually obtained from GMatOrderConstrained(). . colIS - The column permutation Level: developer */ int GMatReorderForNonzeroDiagonal_Private(GMat mat, int rowStart, int rowEnd, double atol, IS rowIS, IS colIS) { int *rowPerm; int *colPerm; int m, n, nz; int *cols; PetscScalar *vals; int replCol; /* Replacement column in the original matrix */ double replVal; /* Replacement pivot magnitude */ int prow; /* Pivot row in the original matrix */ int newRow; /* Pivot row in the permuted matrix */ int temp; int colIndex; int ierr; PetscFunctionBegin; ierr = ISGetIndices(rowIS, &rowPerm); CHKERRQ(ierr); ierr = ISGetIndices(colIS, &colPerm); CHKERRQ(ierr); ierr = MatGetSize(mat, &m, &n); CHKERRQ(ierr); for(prow = 0; prow < m; prow++) { newRow = rowPerm[prow]; /* Row in the permuted matrix */ /* Act only on a block in the reordered matrix */ if ((newRow < rowStart) || (newRow >= rowEnd)) continue; ierr = MatGetRow(mat, prow, &nz, &cols, &vals); CHKERRQ(ierr); /* Find diagonal element */ for(colIndex = 0; colIndex < nz; colIndex++) if (colPerm[cols[colIndex]] == newRow) break; /* Check for zero pivot */ if ((colIndex >= nz) || (PetscAbsScalar(vals[colIndex]) <= atol)) { /* Find the best candidate in the upper triangular part */ replCol = prow; if (colIndex >= nz) replVal = 0.0; else replVal = PetscAbsScalar(vals[colIndex]); for(colIndex = 0; colIndex < nz; colIndex++) { /* Stay within block and upper triangle */ if ((colPerm[cols[colIndex]] <= newRow) || (colPerm[cols[colIndex]] >= rowEnd)) continue; /* Check for acceptable pivot */ if (PetscAbsScalar(vals[colIndex]) > atol) { replCol = cols[colIndex]; replVal = PetscAbsScalar(vals[colIndex]); break; } } /* No candidate was found */ if (prow == replCol) { /* Now we need to look for an element that allows us to pivot with a previous column. To do this, we need to be sure that we don't introduce a zero in a previous diagonal */ if (GMatFindPreviousPivot_Private(mat, prow, rowPerm, colPerm, atol, &replCol, &replVal)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Can not reorder matrix for nonzero diagonal"); } } /* Interchange columns */ temp = colPerm[prow]; colPerm[prow] = colPerm[replCol]; colPerm[replCol] = temp; } ierr = MatRestoreRow(mat, prow, &nz, &cols, &vals); CHKERRQ(ierr); } ierr = ISRestoreIndices(rowIS, &rowPerm); CHKERRQ(ierr); ierr = ISRestoreIndices(colIS, &colPerm); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatReorderForNonzeroDiagonalConstrained" /*@ GMatReorderForNonzeroDiagonalConstrained - Changes matrix ordering to remove zeros from diagonal. This may help in the LU factorization to prevent a zero pivot. Input Parameters: . mat - The grid matrix to reorder . atol - The smallest acceptable pivot . rowIS - The row permutation, usually obtained from GMatOrderConstrained(). . colIS - The column permutation Notes: This is not intended as a replacement for pivoting for matrices that have ``bad'' structure. It is only a stop-gap measure. Should be called after a call to MatGetReordering(), this routine changes the column ordering defined in cis. Options Database Keys: (When using SLES) . -gmat_order_nonzero_diagonal Algorithm: Column pivoting is used. Choice of column is made by looking at the non-zero elements in the row. This algorithm is simple and fast but does NOT guarantee that a non-singular or well conditioned principle submatrix will be produced. Level: developer @*/ int GMatReorderForNonzeroDiagonalConstrained(GMat mat, double atol, IS rowIS, IS colIS) { Grid grid; int m; int mInt; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); PetscValidHeaderSpecific(rowIS, IS_COOKIE); PetscValidHeaderSpecific(colIS, IS_COOKIE); ierr = GMatGetGrid(mat, &grid); CHKERRQ(ierr); if (grid->isConstrained == PETSC_FALSE) { ierr = MatReorderForNonzeroDiagonal(mat, atol, rowIS, colIS); CHKERRQ(ierr); PetscFunctionReturn(0); } /* Get size of matrix blocks */ mInt = grid->constraintOrder->numLocVars - grid->constraintOrder->numLocNewVars; m = grid->order->numVars; /* Reorder interior block */ ierr = GMatReorderForNonzeroDiagonal_Private(mat, 0, mInt, atol, rowIS, colIS); CHKERRQ(ierr); /* Reorder constrained block */ ierr = GMatReorderForNonzeroDiagonal_Private(mat, mInt, m, atol, rowIS, colIS); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatReorder" /*@ GMatReorder - Reorders the matrix based on the ordering type provided. Collective on GMat Input Parameters: . mat - The grid matrix to reorder . rowIS - The row permutation, usually obtained from MatGetOrdering() . colIS - The column permutation . sparse - The flag for sparsification . bw - [Optional] The sparsified bandwidth, PETSC_DECIDE gives the default . frac - [Optional] The sparsified fractional bandwidth, 0.0 gives the default . tol - [Optional] The sparsification drop tolerance, 0.0 is the default Output Parameter: . newmat - The reordered, and possibly sparsified, matrix Options Database Keys: . -gmat_order_nonzero_diagonal Level: advanced .keywords: grid matrix, ordering .seealso: MatGetOrdering(), MatPermute(), MatPermuteSparsify() @*/ int GMatReorder(GMat mat, IS rowIS, IS colIS, PetscTruth sparse, int bw, double frac, double tol, GMat *newmat) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); PetscValidHeaderSpecific(rowIS, IS_COOKIE); PetscValidHeaderSpecific(colIS, IS_COOKIE); PetscValidPointer(newmat); if (sparse == PETSC_TRUE) { ierr = MatPermuteSparsify(mat, bw, frac, tol, rowIS, colIS, newmat); CHKERRQ(ierr); } else { ierr = MatPermute(mat, rowIS, colIS, newmat); CHKERRQ(ierr); } ierr = GMatGetGrid(mat, &grid); CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *newmat, "Grid", (PetscObject) grid); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatEvaluateALEOperatorGalerkin" /*@ GMatEvaluateALEOperatorGalerkin - Evaluates the weak form of an operator over a field on the locations defined by the underlying grid and its discretization. Collective on GMat Input Parameters: + M - The grid matrix . numFields - The number of fields in sFields and tFields . sFields - The shape function fields . sOrder - The global variable ordering for the shape functions . sLocOrder - The local variable ordering for the shape functions . tFields - The test function fields . tOrder - The global variable ordering for the test functions . tLocOrder - The local variable ordering for the test functions . op - The operator . alpha - A scalar multiple for the operator . type - The matrix assembly type - ctx - An optional user provided context for the function Level: intermediate .keywords: grid matrix, evaluate, ALE, operator .seealso: GMatEvaluateSystemMatrix() @*/ int GMatEvaluateALEOperatorGalerkin(GMat M, int numFields, int *sFields, VarOrdering sOrder, LocalVarOrdering sLocOrder, int *tFields, VarOrdering tOrder, LocalVarOrdering tLocOrder, int op, PetscScalar alpha, MatAssemblyType type, void *ctx) { Grid grid; int field; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(M, MAT_COOKIE); PetscValidHeaderSpecific(sOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(sLocOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(tOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(tLocOrder, VAR_ORDER_COOKIE); PetscValidIntPointer(sFields); PetscValidIntPointer(tFields); ierr = GMatGetGrid(M, &grid); CHKERRQ(ierr); PetscValidHeaderSpecific(grid, GRID_COOKIE); for(field = 0; field < numFields; field++) { /* Check for valid fields */ if ((sFields[field] < 0) || (sFields[field] > grid->numFields)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid field number"); } if ((tFields[field] < 0) || (tFields[field] > grid->numFields)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid field number"); } /* Check for valid operator */ if ((op < 0) || (op >= grid->fields[sFields[field]].disc->numOps) || (!grid->fields[sFields[field]].disc->operators[op])) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid operator"); } if ((grid->fields[sFields[field]].disc->operators[op]->precompInt == PETSC_NULL) && (grid->fields[sFields[field]].disc->operators[op]->opFunc == PETSC_NULL) && (grid->fields[sFields[field]].disc->operators[op]->ALEOpFunc == PETSC_NULL)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid operator"); } if (grid->fields[sFields[field]].disc->numQuadPoints != grid->fields[tFields[field]].disc->numQuadPoints) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible quadrature sets"); } } /* Check for compatible orderings */ if ((tOrder->numVars != M->M) || (tOrder->numLocVars != M->m)) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function ordering"); } if ((sOrder->numVars != M->N) || (sOrder->numLocVars != M->n)) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible shape function ordering"); } ierr = PetscLogEventBegin(GMAT_EvaluateOperatorGalerkin, M, grid, sOrder, tOrder); CHKERRQ(ierr); ierr = (*grid->ops->gmatevaluatealeoperatorgalerkin)(grid, M, numFields, sFields, sOrder, sLocOrder, tFields, tOrder, tLocOrder, op, alpha, type, ctx); CHKERRQ(ierr); ierr = PetscLogEventEnd(GMAT_EvaluateOperatorGalerkin, M, grid, sOrder, tOrder); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatEvaluateOperatorGalerkin" /*@ GMatEvaluateOperatorGalerkin - Evaluates the weak form of an operator over a field on the locations defined by the underlying grid and its discretization. Collective on GMat Input Parameter: + M - The grid matrix . x - The argument vector . numFields - The number of fields in sFields and tFields . sFields - The shape function fields . tFields - The test function fields . op - The operator . alpha - The scalar multiple for the operator . type - The matrix assembly type - ctx - [Optional] A user provided context for the function Level: intermediate .keywords: grid matrix, evaluate, operator, galerkin .seealso: GMatEvaluateOperatorGalerkin(), GMatEvaluateSystemMatrix() @*/ int GMatEvaluateOperatorGalerkin(GMat M, GVec x, int numFields, int *sFields, int *tFields, int op, PetscScalar alpha, MatAssemblyType type, void *ctx) { Grid grid; VarOrdering sOldOrder, tOldOrder; VarOrdering sOrder, tOrder; LocalVarOrdering sLocOrder, tLocOrder; VarOrdering xOrder; PetscTruth compat; int numTotalFields; int f; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(M, MAT_COOKIE); PetscValidIntPointer(sFields); PetscValidIntPointer(tFields); ierr = GMatGetGrid(M, &grid); CHKERRQ(ierr); ierr = GMatGetOrder(M, &tOldOrder, &sOldOrder); CHKERRQ(ierr); if (x != PETSC_NULL) { PetscValidHeaderSpecific(x, VEC_COOKIE); ierr = GVecGetOrder(x, &xOrder); CHKERRQ(ierr); ierr = VarOrderingCompatible(xOrder, tOldOrder, &compat); CHKERRQ(ierr); if (compat == PETSC_FALSE) SETERRQ(PETSC_ERR_ARG_INCOMP, "Matrix and vector must have compatible maps"); } ierr = GridGetNumFields(grid, &numTotalFields); CHKERRQ(ierr); for(f = 0; f < numFields; f++) { /* Check for valid fields */ if ((sFields[f] < 0) || (sFields[f] >= numTotalFields)) SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid field number %d", sFields[f]); if ((tFields[f] < 0) || (tFields[f] >= numTotalFields)) SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid field number %d", tFields[f]); /* Check for valid operator */ if ((op < 0) || (op >= grid->fields[sFields[f]].disc->numOps) || (!grid->fields[sFields[f]].disc->operators[op])) { SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid operator %d", op); } if ((grid->fields[sFields[f]].disc->operators[op]->precompInt == PETSC_NULL) && (grid->fields[sFields[f]].disc->operators[op]->opFunc == PETSC_NULL) && (grid->fields[sFields[f]].disc->operators[op]->ALEOpFunc == PETSC_NULL)) { SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid operator %d", op); } if (grid->fields[sFields[f]].disc->numQuadPoints != grid->fields[tFields[f]].disc->numQuadPoints) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible quadrature sets"); } } /* Create orderings */ ierr = VarOrderingCreateSubset(sOldOrder, numFields, sFields, PETSC_FALSE, &sOrder); CHKERRQ(ierr); ierr = LocalVarOrderingCreate(grid, numFields, sFields, &sLocOrder); CHKERRQ(ierr); ierr = VarOrderingCreateSubset(tOldOrder, numFields, tFields, PETSC_FALSE, &tOrder); CHKERRQ(ierr); ierr = LocalVarOrderingCreate(grid, numFields, tFields, &tLocOrder); CHKERRQ(ierr); /* Check for compatible orderings */ if ((tOrder->numVars != M->M) || (tOrder->numLocVars != M->m)) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function ordering"); } if ((sOrder->numVars != M->N) || (sOrder->numLocVars != M->n)) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible shape function ordering"); } /* Calculate operator */ ierr = PetscLogEventBegin(GMAT_EvaluateOperatorGalerkin, M, grid, sOrder, tOrder); CHKERRQ(ierr); ierr = (*grid->ops->gmatevaluateoperatorgalerkin)(grid, M, x, sOrder, sLocOrder, tOrder, tLocOrder, op, alpha, type, ctx); CHKERRQ(ierr); ierr = PetscLogEventEnd(GMAT_EvaluateOperatorGalerkin, M, grid, sOrder, tOrder); CHKERRQ(ierr); /* Destroy orderings */ ierr = VarOrderingDestroy(sOrder); CHKERRQ(ierr); ierr = LocalVarOrderingDestroy(sLocOrder); CHKERRQ(ierr); ierr = VarOrderingDestroy(tOrder); CHKERRQ(ierr); ierr = LocalVarOrderingDestroy(tLocOrder); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatEvaluateALEConstrainedOperatorGalerkin" /*@ GMatEvaluateALEConstrainedOperatorGalerkin - Evaluates the weak form of an operator over a field on the locations defined by the underlying grid and its discretization. The constrained variable space is used for the shape and test functions. Collective on GMat Input Parameters: + M - The grid matrix . numFields - The number of fields in sFields and tFields . sFields - The shape function fields . sOrder - The global variable ordering for the shape functions . sLocOrder - The local variable ordering for the shape functions . tFields - The test function fields . tOrder - The global variable ordering for the test functions . tLocOrder - The local variable ordering for the test functions . op - The operator . alpha - A scalar multiple for the operator . type - The matrix assembly type - ctx - An optional user provided context for the function Level: intermediate .keywords: grid matrix, evaluate, ALE, operator, constraint .seealso: GMatEvaluateALEOperatorGalerkin(), GMatEvaluateSystemMatrix() @*/ int GMatEvaluateALEConstrainedOperatorGalerkin(GMat M, int numFields, int *sFields, VarOrdering sOrder, LocalVarOrdering sLocOrder, int *tFields, VarOrdering tOrder, LocalVarOrdering tLocOrder, int op, PetscScalar alpha, MatAssemblyType type, void *ctx) { Grid grid; int field; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(M, MAT_COOKIE); PetscValidHeaderSpecific(sOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(sLocOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(tOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(tLocOrder, VAR_ORDER_COOKIE); PetscValidIntPointer(sFields); PetscValidIntPointer(tFields); ierr = GMatGetGrid(M, &grid); CHKERRQ(ierr); PetscValidHeaderSpecific(grid, GRID_COOKIE); if (grid->isConstrained == PETSC_FALSE) { ierr = GMatEvaluateALEOperatorGalerkin(M, numFields, sFields, sOrder, sLocOrder, tFields, tOrder, tLocOrder, op, alpha, type, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } for(field = 0; field < numFields; field++) { /* Check for valid fields */ if ((sFields[field] < 0) || (sFields[field] > grid->numFields)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid field number"); } if ((tFields[field] < 0) || (tFields[field] > grid->numFields)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid field number"); } /* Check for valid operator */ if ((op < 0) || (op >= grid->fields[sFields[field]].disc->numOps) || (!grid->fields[sFields[field]].disc->operators[op])) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid operator"); } if ((grid->fields[sFields[field]].disc->operators[op]->precompInt == PETSC_NULL) && (grid->fields[sFields[field]].disc->operators[op]->opFunc == PETSC_NULL) && (grid->fields[sFields[field]].disc->operators[op]->ALEOpFunc == PETSC_NULL)) { SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid operator"); } if (grid->fields[sFields[field]].disc->numQuadPoints != grid->fields[tFields[field]].disc->numQuadPoints) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible quadrature sets"); } } /* Check for compatible orderings */ if ((tOrder->numVars != M->M) || (tOrder->numLocVars != M->m)) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function ordering"); } if ((sOrder->numVars != M->N) || (sOrder->numLocVars != M->n)) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible shape function ordering"); } ierr = (*grid->ops->gmatevaluatealeconstrainedoperatorgalerkin)(grid, M, numFields, sFields, sOrder, sLocOrder, tFields, tOrder, tLocOrder, op, alpha, type, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatEvaluateBoundaryOperatorGalerkin" /*@ GMatEvaluateBoundaryOperatorGalerkin - Evaluates the weak form of an operator over a field on the locations defined by the boundary of the underlying grid and its discretization. The test function are defined on the entire grid, but the only nonzeros come from functions with support on the boundary. Input Parameter: . M - The grid matrix . numFields - The number of fields in sFields and tFields . sFields - The shape function fields . sOrder - The global variable ordering for the shape functions . sLocOrder - The local variable ordering for the shape functions . tFields - The test function fields . tOrder - The global variable ordering for the test functions . tLocOrder - The local variable ordering for the test functions . op - The operator . alpha - A scalar multiple for the operator . type - The matrix assembly type . ctx - An optional user provided context for the function Level: intermediate .seealso: GMatEvaluateOperatorGalerkin, GMatEvaluateSystemMatrix @*/ int GMatEvaluateBoundaryOperatorGalerkin(GMat M, int numFields, int *sFields, VarOrdering sOrder, LocalVarOrdering sLocOrder, int *tFields, VarOrdering tOrder, LocalVarOrdering tLocOrder, int op, PetscScalar alpha, MatAssemblyType type, void *ctx) { Grid grid; int numTotalFields; int f; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(M, MAT_COOKIE); PetscValidHeaderSpecific(sOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(sLocOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(tOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(tLocOrder, VAR_ORDER_COOKIE); PetscValidIntPointer(sFields); PetscValidIntPointer(tFields); ierr = GMatGetGrid(M, &grid); CHKERRQ(ierr); ierr = GridGetNumFields(grid, &numTotalFields); CHKERRQ(ierr); for(f = 0; f < numFields; f++) { /* Check for valid fields */ if ((sFields[f] < 0) || (sFields[f] >= numTotalFields)) SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid field number %d", sFields[f]); if ((tFields[f] < 0) || (tFields[f] >= numTotalFields)) SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid field number %d", tFields[f]); /* Check for valid operator */ if ((op < 0) || (op >= grid->fields[sFields[f]].disc->numOps) || (!grid->fields[sFields[f]].disc->operators[op])) { SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid operator %d", op); } if ((grid->fields[sFields[f]].disc->operators[op]->precompInt == PETSC_NULL) && (grid->fields[sFields[f]].disc->operators[op]->opFunc == PETSC_NULL) && (grid->fields[sFields[f]].disc->operators[op]->ALEOpFunc == PETSC_NULL)) { SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid operator %d", op); } if (grid->fields[sFields[f]].disc->numQuadPoints != grid->fields[tFields[f]].disc->numQuadPoints) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible quadrature sets"); } } /* Check for compatible orderings */ if ((tOrder->numVars != M->M) || (tOrder->numLocVars != M->m)) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible test function ordering"); } if ((sOrder->numVars != M->N) || (sOrder->numLocVars != M->n)) { SETERRQ(PETSC_ERR_ARG_INCOMP, "Incompatible shape function ordering"); } ierr = (*grid->ops->gmatevaluateboundaryoperatorgalerkin)(grid, M, PETSC_NULL, sOrder, sLocOrder, tOrder, tLocOrder, op, alpha, type, ctx); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatEvaluateNewFields" /*@ GMatEvaluateNewFields - Evaluates the weak form of the system operators over the new fields introduced by constraints. These fields are generally not discretized using the computational mesh. Collective on GMat Input Parameters: + M - The grid matrix - type - The matrix assembly type Level: intermediate .keywords grid matrix, evaluate, new field .seealso: GMatEvaluateOperatorGalerkin() @*/ int GMatEvaluateNewFields(GMat M, int numFields, int *sFields, VarOrdering sOrder, LocalVarOrdering sLocOrder, int *tFields, VarOrdering tOrder, LocalVarOrdering tLocOrder, PetscScalar alpha, MatAssemblyType type) { Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(M, MAT_COOKIE); ierr = GMatGetGrid(M, &grid); CHKERRQ(ierr); PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidHeaderSpecific(sOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(sLocOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(tOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(tLocOrder, VAR_ORDER_COOKIE); PetscValidIntPointer(sFields); PetscValidIntPointer(tFields); if (numFields != 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Only one constrained field allowed currently"); if (sFields[0] != tFields[0]) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Must form with constrained field currently"); if (grid->fields[sFields[0]].isConstrained == PETSC_FALSE) { SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Only valid with constrained field currently"); } #if 0 ierr = PetscLogEventBegin(GMAT_EvaluateSystemMatrix, M, grid, 0, 0); CHKERRQ(ierr); ierr = GMatEvaluateNewFields_Triangular_2D(grid, M, numFields, sFields, sOrder, sLocOrder, tFields, tOrder, tLocOrder, alpha, type, grid->constraintCtx); CHKERRQ(ierr); ierr = PetscLogEventEnd(GMAT_EvaluateSystemMatrix, M, grid, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); #else SETERRQ(PETSC_ERR_SUP, "Coming soon"); #endif } #undef __FUNCT__ #define __FUNCT__ "GMatMatMultConstrained" /*@ GMatMatMultConstrained - This function applies P^T A P to a vector, where P is the constraint matrix for the grid. Input Parameters: . A - The grid matrix . x - The input grid vector Output Parameter: . y - The output grid vector P^T A P x Level: intermediate .seealso: GMatEvaluateOperatorGalerkin @*/ int GMatMatMultConstrained(GMat mat, GVec x, GVec y) { Grid grid; Mat P; PetscConstraintObject ctx; GTSContext GTSctx; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); PetscValidHeaderSpecific(x, VEC_COOKIE); PetscValidHeaderSpecific(y, VEC_COOKIE); ierr = GMatGetGrid(mat, &grid); CHKERRQ(ierr); PetscValidHeaderSpecific(grid, GRID_COOKIE); ierr = GridGetConstraintMatrix(grid, &P); CHKERRQ(ierr); ierr = GridGetConstraintContext(grid, &ctx); CHKERRQ(ierr); ierr = PetscObjectQuery((PetscObject) ctx, "GTSContext", (PetscObject *) >Sctx); CHKERRQ(ierr); ierr = PetscLogEventBegin(GMAT_MatMultConstrained, mat, x, y, 0); CHKERRQ(ierr); /* End the current MatMult() timing */ ierr = PetscLogEventEnd(MAT_Mult, mat, x, y, 0); CHKERRQ(ierr); ierr = MatMult(P, x, GTSctx->work[0]); CHKERRQ(ierr); /* Time the multiply with the unconstrained matrix */ ierr = PetscLogEventBegin(MAT_Mult, mat, GTSctx->work[0], GTSctx->work[1], 0); CHKERRQ(ierr); ierr = MatMultConstrained(mat, GTSctx->work[0], GTSctx->work[1]); CHKERRQ(ierr); ierr = PetscLogEventEnd(MAT_Mult, mat, GTSctx->work[0], GTSctx->work[1], 0); CHKERRQ(ierr); ierr = MatMultTranspose(P, GTSctx->work[1], y); CHKERRQ(ierr); /* Add in new fields */ ierr = (*grid->constraintCtx->ops->applyjac)(grid, x, y); CHKERRQ(ierr); /* Restart the current MatMult() timing */ ierr = PetscLogEventBegin(MAT_Mult, mat, x, y, 0); CHKERRQ(ierr); ierr = PetscLogEventEnd(GMAT_MatMultConstrained, mat, x, y, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatMatMultMF" /*@ GMatMatMultMF - This function applies the operator matrix-free to the vector. Input Parameters: + A - The grid matrix - x - The input grid vector Output Parameter: . y - The output grid vector A x Level: intermediate .seealso: GMatMatMultMFConstrained @*/ int GMatMatMultMF(GMat mat, GVec x, GVec y) { PetscScalar zero = 0.0; Grid grid; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); PetscValidHeaderSpecific(x, VEC_COOKIE); PetscValidHeaderSpecific(y, VEC_COOKIE); ierr = GMatGetGrid(mat, &grid); CHKERRQ(ierr); PetscValidHeaderSpecific(grid, GRID_COOKIE); ierr = VecSet(&zero, y); CHKERRQ(ierr); ierr = GVecEvaluateJacobian(y, grid->matrixFreeArg, x, grid->matrixFreeContext); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatMatMultTransposeConstrained" /*@ GMatMatMultTransposeConstrained - This function applies P^T A^T P to a vector, where P is the constraint matrix for the grid. Input Parameters: . A - The grid matrix . x - The input grid vector Output Parameter: . y - The output grid vector P^T A^T P x Level: intermediate .seealso: GMatEvaluateOperatorGalerkin @*/ int GMatMatMultTransposeConstrained(GMat mat, GVec x, GVec y) { Grid grid; Mat P; PetscConstraintObject ctx; GTSContext GTSctx; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(mat, MAT_COOKIE); PetscValidHeaderSpecific(x, VEC_COOKIE); PetscValidHeaderSpecific(y, VEC_COOKIE); ierr = GMatGetGrid(mat, &grid); CHKERRQ(ierr); PetscValidHeaderSpecific(grid, GRID_COOKIE); ierr = GridGetConstraintMatrix(grid, &P); CHKERRQ(ierr); ierr = GridGetConstraintContext(grid, &ctx); CHKERRQ(ierr); ierr = PetscObjectQuery((PetscObject) ctx, "GTSContext", (PetscObject *) >Sctx); CHKERRQ(ierr); ierr = PetscLogEventBegin(GMAT_MatMultTransposeConstrained, mat, x, y, 0); CHKERRQ(ierr); /* End the current MatMult() timing */ ierr = PetscLogEventEnd(MAT_MultTranspose, mat, x, y, 0); CHKERRQ(ierr); ierr = MatMult(P, x, GTSctx->work[0]); CHKERRQ(ierr); /* Time the multiply with the unconstrained matrix */ ierr = PetscLogEventBegin(MAT_MultTranspose, mat, GTSctx->work[0], GTSctx->work[1], 0); CHKERRQ(ierr); ierr = MatMultTransposeConstrained(mat, GTSctx->work[0], GTSctx->work[1]); CHKERRQ(ierr); ierr = PetscLogEventEnd(MAT_MultTranspose, mat, GTSctx->work[0], GTSctx->work[1], 0); CHKERRQ(ierr); ierr = MatMultTranspose(P, GTSctx->work[1], y); CHKERRQ(ierr); /* Add in new fields */ ierr = (*grid->constraintCtx->ops->applyjac)(grid, x, y); CHKERRQ(ierr); /* Restart the current MatMult() timing */ ierr = PetscLogEventBegin(MAT_MultTranspose, mat, x, y, 0); CHKERRQ(ierr); ierr = PetscLogEventEnd(GMAT_MatMultTransposeConstrained, mat, x, y, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatSetBoundary" /*@ GMatSetBoundary - Applies the specified Dirichlet boundary conditions to this matrix. Input Parameter: . M - The grid matrix . diag - The scalar to place on the diagonal . ctx - An optional user provided context for the function Level: intermediate .seealso: GridSetBC(), GridAddBC() @*/ int GMatSetBoundary(GMat M, PetscScalar diag, void *ctx) { Grid grid; int bc; int num; int *bd; int *field; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(M, MAT_COOKIE); ierr = GMatGetGrid(M, &grid); CHKERRQ(ierr); PetscValidHeaderSpecific(grid, GRID_COOKIE); ierr = PetscLogEventBegin(GMAT_SetBoundary, M, grid, 0, 0); CHKERRQ(ierr); for(bc = 0, num = 0; bc < grid->numBC; bc++) { if (grid->bc[bc].reduce == PETSC_FALSE) num++; } if (num > 0) { ierr = PetscMalloc(num * sizeof(int), &bd); CHKERRQ(ierr); ierr = PetscMalloc(num * sizeof(int), &field); CHKERRQ(ierr); for(bc = 0, num = 0; bc < grid->numBC; bc++) { if (grid->bc[bc].reduce == PETSC_FALSE) { bd[num] = grid->bc[bc].boundary; field[num++] = grid->bc[bc].field; } } ierr = GridSetMatBoundaryRectangular(num, bd, field, diag, grid->order, M, ctx); CHKERRQ(ierr); ierr = PetscFree(bd); CHKERRQ(ierr); ierr = PetscFree(field); CHKERRQ(ierr); } for(bc = 0; bc < grid->numPointBC; bc++) { if (grid->pointBC[bc].reduce == PETSC_FALSE) { ierr = GridSetMatPointBoundary(grid->pointBC[bc].node, grid->pointBC[bc].field, diag, M, ctx); CHKERRQ(ierr); } } ierr = PetscLogEventEnd(GMAT_SetBoundary, M, grid, 0, 0); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatCreate" /*@ GMatCreate - Creates a sparse matrix with the appropriate preallocation of nonzeros for a discretization of an operator on the grid. Input Parameter: . grid - The grid Output Parameters: . gmat - The discrete operator Level: beginner .seealso: GVecCreate() @*/ int GMatCreate(Grid grid, GMat *gmat) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(gmat); ierr = GMatCreateRectangular(grid, grid->order, grid->order, gmat); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatCreateRectangular" /*@ GMatCreateRectangular - Creates a sparse matrix with the appropriate preallocation of nonzeros for a discretization of an operator on the grid. A different set of fields may be speicifed for the shape and test discretizations. Input Parameter: . grid - The grid . sOrder - The shape function ordering . tOrder - The test function ordering Output Parameters: . gmat - The discrete operator Level: beginner .seealso: GMatCreate() @*/ int GMatCreateRectangular(Grid grid, VarOrdering sOrder, VarOrdering tOrder, GMat *gmat) { VarOrdering shapeOrder = sOrder; VarOrdering testOrder = tOrder; int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(gmat); ierr = GridSetUp(grid); CHKERRQ(ierr); if (shapeOrder == PETSC_NULL) shapeOrder = grid->order; if (testOrder == PETSC_NULL) testOrder = grid->order; PetscValidHeaderSpecific(shapeOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(testOrder, VAR_ORDER_COOKIE); ierr = PetscLogEventBegin(GMAT_CreateRectangular, grid, shapeOrder, testOrder, 0); CHKERRQ(ierr); ierr = (*grid->ops->gridcreategmat)(grid, shapeOrder, testOrder, PETSC_FALSE, gmat); CHKERRQ(ierr); ierr = PetscLogEventEnd(GMAT_CreateRectangular, grid, shapeOrder, testOrder, 0); CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *gmat, "RowOrder", (PetscObject) tOrder); CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *gmat, "ColOrder", (PetscObject) sOrder); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatDummy0" int GMatDummy0(GMat mat) { PetscFunctionBegin; PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatCreateMF" /*@ GMatCreateMF - Creates a matrix-free operator. A different set of fields may be speicifed for the shape and test discretizations. Input Parameter: . grid - The grid . sOrder - The shape function ordering . tOrder - The test function ordering Output Parameters: . gmat - The discrete operator Level: beginner .seealso: GMatCreateRectangular() @*/ int GMatCreateMF(Grid grid, VarOrdering sOrder, VarOrdering tOrder, GMat *gmat) { VarOrdering shapeOrder = sOrder; VarOrdering testOrder = tOrder; int ierr; PetscFunctionBegin; if (shapeOrder == PETSC_NULL) shapeOrder = grid->order; if (testOrder == PETSC_NULL) testOrder = grid->order; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidHeaderSpecific(shapeOrder, VAR_ORDER_COOKIE); PetscValidHeaderSpecific(testOrder, VAR_ORDER_COOKIE); PetscValidPointer(gmat); ierr = GridSetUp(grid); CHKERRQ(ierr); ierr = PetscLogEventBegin(GMAT_CreateRectangular, grid, shapeOrder, testOrder, 0); CHKERRQ(ierr); ierr = MatCreateShell(grid->comm, tOrder->numLocVars, sOrder->numLocVars, tOrder->numVars, sOrder->numVars, grid, gmat); CHKERRQ(ierr); if (grid->isConstrained == PETSC_TRUE) { if (grid->explicitConstraints == PETSC_FALSE) { ierr = MatShellSetOperation(*gmat, MATOP_MULT, (void (*)(void)) GMatMatMultConstrained); CHKERRQ(ierr); ierr = MatShellSetOperation(*gmat, MATOP_MULT_CONSTRAINED, (void (*)(void)) GMatMatMultMF); CHKERRQ(ierr); } else { ierr = MatShellSetOperation(*gmat, MATOP_MULT, (void (*)(void)) GMatMatMultMF); CHKERRQ(ierr); } } else { ierr = MatShellSetOperation(*gmat, MATOP_MULT, (void (*)(void)) GMatMatMultMF); CHKERRQ(ierr); } ierr = MatShellSetOperation(*gmat, MATOP_ZERO_ENTRIES, (void (*)(void)) GMatDummy0); CHKERRQ(ierr); ierr = MatShellSetOperation(*gmat, MATOP_GET_DIAGONAL, (void (*)(void)) GMatGetDiagonalMF); CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *gmat, "Grid", (PetscObject) grid); CHKERRQ(ierr); ierr = MatSetOption(*gmat, MAT_NEW_NONZERO_ALLOCATION_ERR); CHKERRQ(ierr); ierr = PetscLogEventEnd(GMAT_CreateRectangular, grid, shapeOrder, testOrder, 0); CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *gmat, "RowOrder", (PetscObject) tOrder); CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *gmat, "ColOrder", (PetscObject) sOrder); CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "GMatCreateBoundaryRestriction" /*@ GMatCreateBoundaryRestriction - Creates a sparse matrix with the appropriate preallocation of nonzeros for a discretization of an operator on the boundary of the grid. Collective on Grid Input Parameter: . grid - The grid Output Parameter: . gmat - The grid matrix Level: beginner .keywords: grid matrix, boundary, create .seealso: GMatCreate(), GVecCreate() @*/ int GMatCreateBoundaryRestriction(Grid grid, GMat *gmat) { int ierr; PetscFunctionBegin; PetscValidHeaderSpecific(grid, GRID_COOKIE); PetscValidPointer(gmat); ierr = GridSetUp(grid); CHKERRQ(ierr); ierr = GridSetupBoundary(grid); CHKERRQ(ierr); ierr = (*grid->ops->gridcreategmat)(grid, grid->bdOrder, grid->order, PETSC_TRUE, gmat); CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *gmat, "RowOrder", (PetscObject) grid->order); CHKERRQ(ierr); ierr = PetscObjectCompose((PetscObject) *gmat, "ColOrder", (PetscObject) grid->bdOrder); CHKERRQ(ierr); PetscFunctionReturn(0); }