[petsc-users] DM misuse causes massive memory leak?
Jed Brown
jed at jedbrown.org
Wed Dec 29 16:55:09 CST 2021
"Ferrand, Jesus A." <FERRANJ2 at my.erau.edu> writes:
> Dear PETSc Team:
>
> I have a question about DM and PetscSection. Say I import a mesh (for FEM purposes) and create a DMPlex for it. I then use PetscSections to set degrees of freedom per "point" (by point I mean vertices, lines, faces, and cells). I then use PetscSectionGetStorageSize() to get the size of the global stiffness matrix (K) needed for my FEM problem. One last detail, this K I populate inside a rather large loop using an element stiffness matrix function of my own. Instead of using DMCreateMatrix(), I manually created a Mat using MatCreate(), MatSetType(), MatSetSizes(), and MatSetUp(). I come to find that said loop is painfully slow when I use the manually created matrix, but 20x faster when I use the Mat coming out of DMCreateMatrix().
The sparse matrix hasn't been preallocated, which forces the data structure to do a lot of copies (as bad as O(n^2) complexity). DMCreateMatrix() preallocates for you.
https://petsc.org/release/docs/manual/performance/#memory-allocation-for-sparse-matrix-assembly
https://petsc.org/release/docs/manual/mat/#sec-matsparse
> My question is then: Is the manual Mat a noob mistake and is it somehow creating a memory leak with K? Just in case it's something else I'm attaching the code. The loop that populates K is between lines 221 and 278. Anything related to DM, DMPlex, and PetscSection is between lines 117 and 180.
>
> Machine Type: HP Laptop
> C-compiler: Gnu C
> OS: Ubuntu 20.04
> PETSc version: 3.16.0
> MPI Implementation: MPICH
>
> Hope you all had a Merry Christmas and that you will have a happy and productive New Year. :D
>
>
> Sincerely:
>
> J.A. Ferrand
>
> Embry-Riddle Aeronautical University - Daytona Beach FL
>
> M.Sc. Aerospace Engineering | May 2022
>
> B.Sc. Aerospace Engineering
>
> B.Sc. Computational Mathematics
>
>
>
> Sigma Gamma Tau
>
> Tau Beta Pi
>
> Honors Program
>
>
>
> Phone: (386)-843-1829
>
> Email(s): ferranj2 at my.erau.edu
>
> jesus.ferrand at gmail.com
> //REFERENCE: https://github.com/FreeFem/FreeFem-sources/blob/master/plugin/mpi/PETSc-code.hpp
> #include <petsc.h>
> static char help[] = "Imports a Gmsh mesh with boundary conditions and solves the elasticity equation.\n"
> "Option prefix = opt_.\n";
>
> struct preKE{//Preallocation before computing KE
> Mat matB,
> matBTCB;
> //matKE;
> PetscInt x_insert[3],
> y_insert[3],
> z_insert[3],
> m,//Looping variables.
> sizeKE,//size of the element stiffness matrix.
> N,//Number of nodes in element.
> x_in,y_in,z_in; //LI to index B matrix.
> PetscReal J[3][3],//Jacobian matrix.
> invJ[3][3],//Inverse of the Jacobian matrix.
> detJ,//Determinant of the Jacobian.
> dX[3],
> dY[3],
> dZ[3],
> minor00,
> minor01,
> minor02,//Determinants of minors in a 3x3 matrix.
> dPsidX, dPsidY, dPsidZ,//Shape function derivatives w.r.t global coordinates.
> weight,//Multiplier of quadrature weights.
> *dPsidXi,//Derivatives of shape functions w.r.t. Xi.
> *dPsidEta,//Derivatives of shape functions w.r.t. Eta.
> *dPsidZeta;//Derivatives of shape functions w.r.t Zeta.
> PetscErrorCode ierr;
> };//end struct.
>
> //Function declarations.
> extern PetscErrorCode tetra4(PetscScalar*, PetscScalar*, PetscScalar*,struct preKE*, Mat*, Mat*);
> extern PetscErrorCode ConstitutiveMatrix(Mat*,const char*,PetscInt);
> extern PetscErrorCode InitializeKEpreallocation(struct preKE*,const char*);
>
> PetscErrorCode PetscViewerVTKWriteFunction(PetscObject vec,PetscViewer viewer){
> PetscErrorCode ierr;
> ierr = VecView((Vec)vec,viewer); CHKERRQ(ierr);
> return ierr;
> }
>
>
>
>
> int main(int argc, char **args){
> //DEFINITIONS OF PETSC's DMPLEX LINGO:
> //POINT: A topology element (cell, face, edge, or vertex).
> //CHART: It an interval from 0 to the number of "points." (the range of admissible linear indices)
> //STRATUM: A subset of the "chart" which corresponds to all "points" at a given "level."
> //LEVEL: This is either a "depth" or a "height".
> //HEIGHT: Dimensionality of an element measured from 0D to 3D. Heights: cell = 0, face = 1, edge = 2, vertex = 3.
> //DEPTH: Dimensionality of an element measured from 3D to 0D. Depths: cell = 3, face = 2, edge = 1, vertex = 0;
> //CLOSURE: *of an element is the collection of all other elements that define it.I.e., the closure of a surface is the collection of vertices and edges that make it up.
> //STAR:
> //STANDARD LABELS: These are default tags that DMPlex has for its topology. ("depth")
> PetscErrorCode ierr;//Error tracking variable.
> DM dm;//Distributed memory object (useful for managing grids.)
> DMLabel physicalgroups;//Identifies user-specified tags in gmsh (to impose BC's).
> DMPolytopeType celltype;//When looping through cells, determines its type (tetrahedron, pyramid, hexahedron, etc.)
> PetscSection s;
> KSP ksp;//Krylov Sub-Space (linear solver object)
> Mat K,//Global stiffness matrix (Square, assume unsymmetric).
> KE,//Element stiffness matrix (Square, assume unsymmetric).
> matC;//Constitutive matrix.
> Vec XYZ,//Coordinate vector, contains spatial locations of mesh's vertices (NOTE: This vector self-destroys!).
> U,//Displacement vector.
> F;//Load Vector.
> PetscViewer XYZviewer,//Viewer object to output mesh to ASCII format.
> XYZpUviewer; //Viewer object to output displacements to ASCII format.
> PetscBool interpolate = PETSC_TRUE,//Instructs Gmsh importer whether to generate faces and edges (Needed when using P2 or higher elements).
> useCone = PETSC_TRUE,//Instructs "DMPlexGetTransitiveClosure()" whether to extract the closure or the star.
> dirichletBC = PETSC_FALSE,//For use when assembling the K matrix.
> neumannBC = PETSC_FALSE,//For use when assembling the F vector.
> saveASCII = PETSC_FALSE,//Whether to save results in ASCII format.
> saveVTK = PETSC_FALSE;//Whether to save results as VTK format.
> PetscInt nc,//number of cells. (PETSc lingo for "elements")
> nv,//number of vertices. (PETSc lingo for "nodes")
> nf,//number of faces. (PETSc lingo for "surfaces")
> ne,//number of edges. (PETSc lingo for "lines")
> pStart,//starting LI of global elements.
> pEnd,//ending LI of all elements.
> cStart,//starting LI for cells global arrangement.
> cEnd,//ending LI for cells in global arrangement.
> vStart,//starting LI for vertices in global arrangement.
> vEnd,//ending LI for vertices in global arrangement.
> fStart,//starting LI for faces in global arrangement.
> fEnd,//ending LI for faces in global arrangement.
> eStart,//starting LI for edges in global arrangement.
> eEnd,//ending LI for edges in global arrangement.
> sizeK,//Size of the element stiffness matrix.
> ii,jj,kk,//Dedicated looping variables.
> indexXYZ,//Variable to access the elements of XYZ vector.
> indexK,//Variable to access the elements of the U and F vectors (can reference rows and colums of K matrix.)
> *closure = PETSC_NULL,//Pointer to the closure elements of a cell.
> size_closure,//Size of the closure of a cell.
> dim,//Dimension of the mesh.
> //*edof,//Linear indices of dof's inside the K matrix.
> dof = 3,//Degrees of freedom per node.
> cells=0, edges=0, vertices=0, faces=0,//Topology counters when looping through cells.
> pinXcode=10, pinZcode=11,forceZcode=12;//Gmsh codes to extract relevant "Face Sets."
> PetscReal //*x_el,//Pointer to a vector that will store the x-coordinates of an element's vertices.
> //*y_el,//Pointer to a vector that will store the y-coordinates of an element's vertices.
> //*z_el,//Pointer to a vector that will store the z-coordinates of an element's vertices.
> *xyz_el,//Pointer to xyz array in the XYZ vector.
> traction = -10,
> *KEdata,
> t1,t2; //time keepers.
> const char *gmshfile = "TopOptmeshfine2.msh";//Name of gmsh file to import.
>
> ierr = PetscInitialize(&argc,&args,NULL,help); if(ierr) return ierr; //And the machine shall work....
>
> //MESH IMPORT=================================================================
> //IMPORTANT NOTE: Gmsh only creates "cells" and "vertices," it does not create the "faces" or "edges."
> //Gmsh probably can generate them, must figure out how to.
> t1 = MPI_Wtime();
> ierr = DMPlexCreateGmshFromFile(PETSC_COMM_WORLD,gmshfile,interpolate,&dm); CHKERRQ(ierr);//Read Gmsh file and generate the DMPlex.
> ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr);//1-D, 2-D, or 3-D
> ierr = DMPlexGetChart(dm, &pStart, &pEnd); CHKERRQ(ierr);//Extracts linear indices of cells, vertices, faces, and edges.
> ierr = DMGetCoordinatesLocal(dm,&XYZ); CHKERRQ(ierr);//Extracts coordinates from mesh.(Contiguous storage: [x0,y0,z0,x1,y1,z1,...])
> ierr = VecGetArray(XYZ,&xyz_el); CHKERRQ(ierr);//Get pointer to vector's data.
> t2 = MPI_Wtime();
> PetscPrintf(PETSC_COMM_WORLD,"Mesh Import time: %10f\n",t2-t1);
> DMView(dm,PETSC_VIEWER_STDOUT_WORLD);
>
> //IMPORTANT NOTE: PETSc assumes that vertex-cell meshes are 2D even if they were 3D, so its ordering changes.
> //Cells remain at height 0, but vertices move to height 1 from height 3. To prevent this from becoming an issue
> //the "interpolate" variable is set to PETSC_TRUE to tell the mesh importer to generate faces and edges.
> //PETSc, therefore, technically does additional meshing. Gotta figure out how to get this from Gmsh directly.
> ierr = DMPlexGetDepthStratum(dm,3, &cStart, &cEnd);//Get LI of cells.
> ierr = DMPlexGetDepthStratum(dm,2, &fStart, &fEnd);//Get LI of faces
> ierr = DMPlexGetDepthStratum(dm,1, &eStart, &eEnd);//Get LI of edges.
> ierr = DMPlexGetDepthStratum(dm,0, &vStart, &vEnd);//Get LI of vertices.
> ierr = DMGetStratumSize(dm,"depth", 3, &nc);//Get number of cells.
> ierr = DMGetStratumSize(dm,"depth", 2, &nf);//Get number of faces.
> ierr = DMGetStratumSize(dm,"depth", 1, &ne);//Get number of edges.
> ierr = DMGetStratumSize(dm,"depth", 0, &nv);//Get number of vertices.
> /*
> PetscPrintf(PETSC_COMM_WORLD,"global start = %10d\t global end = %10d\n",pStart,pEnd);
> PetscPrintf(PETSC_COMM_WORLD,"#cells = %10d\t i = %10d\t i < %10d\n",nc,cStart,cEnd);
> PetscPrintf(PETSC_COMM_WORLD,"#faces = %10d\t i = %10d\t i < %10d\n",nf,fStart,fEnd);
> PetscPrintf(PETSC_COMM_WORLD,"#edges = %10d\t i = %10d\t i < %10d\n",ne,eStart,eEnd);
> PetscPrintf(PETSC_COMM_WORLD,"#vertices = %10d\t i = %10d\t i < %10d\n",nv,vStart,vEnd);
> */
> //MESH IMPORT=================================================================
>
> //NOTE: This section extremely hardcoded right now.
> //Current setup would only support P1 meshes.
> //MEMORY ALLOCATION ==========================================================
> ierr = PetscSectionCreate(PETSC_COMM_WORLD, &s); CHKERRQ(ierr);
> //The chart is akin to a contiguous memory storage allocation. Each chart entry is associated
> //with a "thing," could be a vertex, face, cell, or edge, or anything else.
> ierr = PetscSectionSetChart(s, pStart, pEnd); CHKERRQ(ierr);
> //For each "thing" in the chart, additional room can be made. This is helpful for associating
> //nodes to multiple degrees of freedom. These commands help associate nodes with
> for(ii = cStart; ii < cEnd; ii++){//Cell loop.
> ierr = PetscSectionSetDof(s, ii, 0);CHKERRQ(ierr);}//NOTE: Currently no dof's associated with cells.
> for(ii = fStart; ii < fEnd; ii++){//Face loop.
> ierr = PetscSectionSetDof(s, ii, 0);CHKERRQ(ierr);}//NOTE: Currently no dof's associated with faces.
> for(ii = vStart; ii < vEnd; ii++){//Vertex loop.
> ierr = PetscSectionSetDof(s, ii, dof);CHKERRQ(ierr);}//Sets x, y, and z displacements as dofs.
> for(ii = eStart; ii < eEnd; ii++){//Edge loop
> ierr = PetscSectionSetDof(s, ii, 0);CHKERRQ(ierr);}//NOTE: Currently no dof's associated with edges.
> ierr = PetscSectionSetUp(s); CHKERRQ(ierr);
> ierr = PetscSectionGetStorageSize(s,&sizeK);CHKERRQ(ierr);//Determine the size of the global stiffness matrix.
> ierr = DMSetLocalSection(dm,s); CHKERRQ(ierr);//Associate the PetscSection with the DM object.
> //PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec)
> //ierr = DMCreateGlobalVector(dm,&U); CHKERRQ(ierr);
> PetscSectionDestroy(&s);
> //PetscPrintf(PETSC_COMM_WORLD,"sizeK = %10d\n",sizeK);
>
> //OBJECT SETUP================================================================
> //Global stiffness matrix.
> //PetscErrorCode DMCreateMatrix(DM dm,Mat *mat)
>
> //This makes the loop fast.
> ierr = DMCreateMatrix(dm,&K);
>
> //This makes the loop uber slow.
> //ierr = MatCreate(PETSC_COMM_WORLD,&K); CHKERRQ(ierr);
> //ierr = MatSetType(K,MATAIJ); CHKERRQ(ierr);// Global stiffness matrix set to some sparse type.
> //ierr = MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,sizeK,sizeK); CHKERRQ(ierr);
> //ierr = MatSetUp(K); CHKERRQ(ierr);
>
> //Displacement vector.
> ierr = VecCreate(PETSC_COMM_WORLD,&U); CHKERRQ(ierr);
> ierr = VecSetType(U,VECSTANDARD); CHKERRQ(ierr);// Global stiffness matrix set to some sparse type.
> ierr = VecSetSizes(U,PETSC_DECIDE,sizeK); CHKERRQ(ierr);
>
> //Load vector.
> ierr = VecCreate(PETSC_COMM_WORLD,&F); CHKERRQ(ierr);
> ierr = VecSetType(F,VECSTANDARD); CHKERRQ(ierr);// Global stiffness matrix set to some sparse type.
> ierr = VecSetSizes(F,PETSC_DECIDE,sizeK); CHKERRQ(ierr);
> //OBJECT SETUP================================================================
>
> //WARNING: This loop is currently hardcoded for P1 elements only! Must Figure
> //out a clever way to modify to accomodate Pn (n>1) elements.
>
> //BEGIN GLOBAL STIFFNESS MATRIX BUILDER=======================================
> t1 = MPI_Wtime();
>
> //PREALLOCATIONS==============================================================
> ierr = ConstitutiveMatrix(&matC,"isotropic",0); CHKERRQ(ierr);
> struct preKE preKEtetra4;
> ierr = InitializeKEpreallocation(&preKEtetra4,"tetra4"); CHKERRQ(ierr);
> ierr = MatCreate(PETSC_COMM_WORLD,&KE); CHKERRQ(ierr); //SEQUENTIAL
> ierr = MatSetSizes(KE,PETSC_DECIDE,PETSC_DECIDE,12,12); CHKERRQ(ierr); //SEQUENTIAL
> ierr = MatSetType(KE,MATDENSE); CHKERRQ(ierr); //SEQUENTIAL
> ierr = MatSetUp(KE); CHKERRQ(ierr);
> PetscReal x_tetra4[4], y_tetra4[4],z_tetra4[4],
> x_hex8[8], y_hex8[8],z_hex8[8],
> *x,*y,*z;
> PetscInt *EDOF,edof_tetra4[12],edof_hex8[24];
> DMPolytopeType previous = DM_POLYTOPE_UNKNOWN;
> //PREALLOCATIONS==============================================================
>
>
>
> for(ii=cStart;ii<cEnd;ii++){//loop through cells.
> ierr = DMPlexGetTransitiveClosure(dm, ii, useCone, &size_closure, &closure); CHKERRQ(ierr);
> ierr = DMPlexGetCellType(dm, ii, &celltype); CHKERRQ(ierr);
> //IMPORTANT NOTE: MOST OF THIS LOOP SHOULD BE INCLUDED IN THE KE3D function.
> if(previous != celltype){
> //PetscPrintf(PETSC_COMM_WORLD,"run \n");
> if(celltype == DM_POLYTOPE_TETRAHEDRON){
> x = x_tetra4;
> y = y_tetra4;
> z = z_tetra4;
> EDOF = edof_tetra4;
> }//end if.
> else if(celltype == DM_POLYTOPE_HEXAHEDRON){
> x = x_hex8;
> y = y_hex8;
> z = z_hex8;
> EDOF = edof_hex8;
> }//end else if.
> }
> previous = celltype;
>
> //PetscPrintf(PETSC_COMM_WORLD,"Cell # %4i\t",ii);
> cells=0;
> edges=0;
> vertices=0;
> faces=0;
> kk = 0;
> for(jj=0;jj<(2*size_closure);jj+=2){//Scan the closure of the current cell.
> //Use information from the DM's strata to determine composition of cell_ii.
> if(vStart <= closure[jj] && closure[jj]< vEnd){//Check for vertices.
> //PetscPrintf(PETSC_COMM_WORLD,"%5i\t",closure[jj]);
> indexXYZ = dim*(closure[jj]-vStart);//Linear index of x-coordinate in the xyz_el array.
>
> *(x+vertices) = xyz_el[indexXYZ];
> *(y+vertices) = xyz_el[indexXYZ+1];//Extract Y-coordinates of the current vertex.
> *(z+vertices) = xyz_el[indexXYZ+2];//Extract Y-coordinates of the current vertex.
> *(EDOF + kk) = indexXYZ;
> *(EDOF + kk+1) = indexXYZ+1;
> *(EDOF + kk+2) = indexXYZ+2;
> kk+=3;
> vertices++;//Update vertex counter.
> }//end if
> else if(eStart <= closure[jj] && closure[jj]< eEnd){//Check for edge ID's
> edges++;
> }//end else ifindexK
> else if(fStart <= closure[jj] && closure[jj]< fEnd){//Check for face ID's
> faces++;
> }//end else if
> else if(cStart <= closure[jj] && closure[jj]< cEnd){//Check for cell ID's
> cells++;
> }//end else if
> }//end "jj" loop.
> ierr = tetra4(x,y,z,&preKEtetra4,&matC,&KE); CHKERRQ(ierr); //Generate the element stiffness matrix for this cell.
> ierr = MatDenseGetArray(KE,&KEdata); CHKERRQ(ierr);
> ierr = MatSetValues(K,12,EDOF,12,EDOF,KEdata,ADD_VALUES); CHKERRQ(ierr);//WARNING: HARDCODED FOR TETRAHEDRAL P1 ELEMENTS ONLY !!!!!!!!!!!!!!!!!!!!!!!
> ierr = MatDenseRestoreArray(KE,&KEdata); CHKERRQ(ierr);
> ierr = DMPlexRestoreTransitiveClosure(dm, ii,useCone, &size_closure, &closure); CHKERRQ(ierr);
> }//end "ii" loop.
> ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> //ierr = MatView(K,PETSC_VIEWER_DRAW_WORLD); CHKERRQ(ierr);
> //END GLOBAL STIFFNESS MATRIX BUILDER===========================================
> t2 = MPI_Wtime();
> PetscPrintf(PETSC_COMM_WORLD,"K build time: %10f\n",t2-t1);
>
>
>
>
>
>
>
>
> t1 = MPI_Wtime();
> //BEGIN BOUNDARY CONDITION ENFORCEMENT==========================================
> IS TrianglesIS, physicalsurfaceID;//, VerticesIS;
> PetscInt numsurfvals,
> //numRows,
> dof_offset,numTri;
> const PetscInt *surfvals,
> //*pinZID,
> *TriangleID;
> PetscScalar diag =1;
> PetscReal area,force;
> //NOTE: Petsc can read/assign labels. Eeach label may posses multiple "values."
> //These values act as tags within a tag.
> //IMPORTANT NOTE: The below line needs a safety. If a mesh that does not feature
> //face sets is imported, the code in its current state will crash!!!. This is currently
> //hardcoded for the test mesh.
> ierr = DMGetLabel(dm, "Face Sets", &physicalgroups); CHKERRQ(ierr);//Inspects Physical surface groups defined by gmsh (if any).
> ierr = DMLabelGetValueIS(physicalgroups, &physicalsurfaceID); CHKERRQ(ierr);//Gets the physical surface ID's defined in gmsh (as specified in the .geo file).
> ierr = ISGetIndices(physicalsurfaceID,&surfvals); CHKERRQ(ierr);//Get a pointer to the actual surface values.
> ierr = DMLabelGetNumValues(physicalgroups, &numsurfvals); CHKERRQ(ierr);//Gets the number of different values that the label assigns.
> for(ii=0;ii<numsurfvals;ii++){//loop through the values under the label.
> //PetscPrintf(PETSC_COMM_WORLD,"Values = %5i\n",surfvals[ii]);
> //PROBLEM: The surface values are hardcoded in the gmsh file. We need to adopt standard "codes"
> //that we can give to users when they make their meshes so that this code recognizes the Type
> // of boundary conditions that are to be imposed.
> if(surfvals[ii] == pinXcode){
> dof_offset = 0;
> dirichletBC = PETSC_TRUE;
> }//end if.
> else if(surfvals[ii] == pinZcode){
> dof_offset = 2;
> dirichletBC = PETSC_TRUE;
> }//end else if.
> else if(surfvals[ii] == forceZcode){
> dof_offset = 2;
> neumannBC = PETSC_TRUE;
> }//end else if.
>
> ierr = DMLabelGetStratumIS(physicalgroups, surfvals[ii], &TrianglesIS); CHKERRQ(ierr);//Get the ID's (as an IS) of the surfaces belonging to value 11.
> //PROBLEM: DMPlexGetConeRecursiveVertices returns an array with repeated node ID's. For each repetition, the lines that enforce BC's unnecessarily re-run.
> ierr = ISGetSize(TrianglesIS,&numTri); CHKERRQ(ierr);
> ierr = ISGetIndices(TrianglesIS,&TriangleID); CHKERRQ(ierr);//Get a pointer to the actual surface values.
> for(kk=0;kk<numTri;kk++){
> ierr = DMPlexGetTransitiveClosure(dm, TriangleID[kk], useCone, &size_closure, &closure); CHKERRQ(ierr);
> if(neumannBC){
> ierr = DMPlexComputeCellGeometryFVM(dm, TriangleID[kk], &area,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
> force = traction*area/3;//WARNING: The 3 here is hardcoded for a purely tetrahedral mesh only!!!!!!!!!!
> }
> for(jj=0;jj<(2*size_closure);jj+=2){
> //PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
> if(vStart <= closure[jj] && closure[jj]< vEnd){//Check for vertices.
> indexK = dof*(closure[jj] - vStart) + dof_offset; //Compute the dof ID's in the K matrix.
> if(dirichletBC){//Boundary conditions requiring an edit of K matrix.
> ierr = MatZeroRows(K,1,&indexK,diag,NULL,NULL); CHKERRQ(ierr);
> }//end if.
> else if(neumannBC){//Boundary conditions requiring an edit of RHS vector.
> ierr = VecSetValue(F,indexK,force,ADD_VALUES); CHKERRQ(ierr);
> }// end else if.
> }//end if.
> }//end "jj" loop.
> ierr = DMPlexRestoreTransitiveClosure(dm, closure[jj],useCone, &size_closure, &closure); CHKERRQ(ierr);
> }//end "kk" loop.
> ierr = ISRestoreIndices(TrianglesIS,&TriangleID); CHKERRQ(ierr);
>
> /*
> ierr = DMPlexGetConeRecursiveVertices(dm, TrianglesIS, &VerticesIS); CHKERRQ(ierr);//Get the ID's (as an IS) of the vertices that make up the surfaces of value 11.
> ierr = ISGetSize(VerticesIS,&numRows); CHKERRQ(ierr);//Get number of flagged vertices (this includes repeated indices for faces that share nodes).
> ierr = ISGetIndices(VerticesIS,&pinZID); CHKERRQ(ierr);//Get a pointer to the actual surface values.
> if(dirichletBC){//Boundary conditions requiring an edit of K matrix.
> for(kk=0;kk<numRows;kk++){
> indexK = 3*(pinZID[kk] - vStart) + dof_offset; //Compute the dof ID's in the K matrix. (NOTE: the 3* ishardcoded for 3 degrees of freedom, tie this to a variable in the FUTURE.)
> ierr = MatZeroRows(K,1,&indexK,diag,NULL,NULL); CHKERRQ(ierr);
> }//end "kk" loop.
> }//end if.
> else if(neumannBC){//Boundary conditions requiring an edit of RHS vector.
> for(kk=0;kk<numRows;kk++){
> indexK = 3*(pinZID[kk] - vStart) + dof_offset;
> ierr = VecSetValue(F,indexK,traction,INSERT_VALUES); CHKERRQ(ierr);
> }//end "kk" loop.
> }// end else if.
> ierr = ISRestoreIndices(VerticesIS,&pinZID); CHKERRQ(ierr);
> */
> dirichletBC = PETSC_FALSE;
> neumannBC = PETSC_FALSE;
> }//end "ii" loop.
> ierr = ISRestoreIndices(physicalsurfaceID,&surfvals); CHKERRQ(ierr);
> //ierr = ISRestoreIndices(VerticesIS,&pinZID); CHKERRQ(ierr);
> ierr = ISDestroy(&physicalsurfaceID); CHKERRQ(ierr);
> //ierr = ISDestroy(&VerticesIS); CHKERRQ(ierr);
> ierr = ISDestroy(&TrianglesIS); CHKERRQ(ierr);
> //END BOUNDARY CONDITION ENFORCEMENT============================================
> t2 = MPI_Wtime();
> PetscPrintf(PETSC_COMM_WORLD,"BC imposition time: %10f\n",t2-t1);
>
> /*
> PetscInt kk = 0;
> for(ii=vStart;ii<vEnd;ii++){
> kk++;
> PetscPrintf(PETSC_COMM_WORLD,"Vertex #%4i\t x = %10.9f\ty = %10.9f\tz = %10.9f\n",ii,xyz_el[3*kk],xyz_el[3*kk+1],xyz_el[3*kk+2]);
> }// end "ii" loop.
> */
>
> t1 = MPI_Wtime();
> //SOLVER========================================================================
> ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr);
> ierr = KSPSetOperators(ksp,K,K); CHKERRQ(ierr);
> ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
> ierr = KSPSolve(ksp,F,U); CHKERRQ(ierr);
> t2 = MPI_Wtime();
> //ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
> //SOLVER========================================================================
> t2 = MPI_Wtime();
> PetscPrintf(PETSC_COMM_WORLD,"Solver time: %10f\n",t2-t1);
> ierr = VecRestoreArray(XYZ,&xyz_el); CHKERRQ(ierr);//Get pointer to vector's data.
>
> //BEGIN MAX/MIN DISPLACEMENTS===================================================
> IS ISux,ISuy,ISuz;
> Vec UX,UY,UZ;
> PetscReal UXmax,UYmax,UZmax,UXmin,UYmin,UZmin;
> ierr = ISCreateStride(PETSC_COMM_WORLD,nv,0,3,&ISux); CHKERRQ(ierr);
> ierr = ISCreateStride(PETSC_COMM_WORLD,nv,1,3,&ISuy); CHKERRQ(ierr);
> ierr = ISCreateStride(PETSC_COMM_WORLD,nv,2,3,&ISuz); CHKERRQ(ierr);
>
> //PetscErrorCode VecGetSubVector(Vec X,IS is,Vec *Y)
> ierr = VecGetSubVector(U,ISux,&UX); CHKERRQ(ierr);
> ierr = VecGetSubVector(U,ISuy,&UY); CHKERRQ(ierr);
> ierr = VecGetSubVector(U,ISuz,&UZ); CHKERRQ(ierr);
>
> //PetscErrorCode VecMax(Vec x,PetscInt *p,PetscReal *val)
> ierr = VecMax(UX,PETSC_NULL,&UXmax); CHKERRQ(ierr);
> ierr = VecMax(UY,PETSC_NULL,&UYmax); CHKERRQ(ierr);
> ierr = VecMax(UZ,PETSC_NULL,&UZmax); CHKERRQ(ierr);
>
> ierr = VecMin(UX,PETSC_NULL,&UXmin); CHKERRQ(ierr);
> ierr = VecMin(UY,PETSC_NULL,&UYmin); CHKERRQ(ierr);
> ierr = VecMin(UZ,PETSC_NULL,&UZmin); CHKERRQ(ierr);
>
> PetscPrintf(PETSC_COMM_WORLD,"%10f\t <= ux <= %10f\n",UXmin,UXmax);
> PetscPrintf(PETSC_COMM_WORLD,"%10f\t <= uy <= %10f\n",UYmin,UYmax);
> PetscPrintf(PETSC_COMM_WORLD,"%10f\t <= uz <= %10f\n",UZmin,UZmax);
>
>
>
>
> //BEGIN OUTPUT SOLUTION=========================================================
> if(saveASCII){
> ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"XYZ.txt",&XYZviewer);
> ierr = VecView(XYZ,XYZviewer); CHKERRQ(ierr);
> ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"U.txt",&XYZpUviewer);
> ierr = VecView(U,XYZpUviewer); CHKERRQ(ierr);
> PetscViewerDestroy(&XYZviewer); PetscViewerDestroy(&XYZpUviewer);
>
> }//end if.
> if(saveVTK){
> const char *meshfile = "starting_mesh.vtk",
> *deformedfile = "deformed_mesh.vtk";
> ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,meshfile,FILE_MODE_WRITE,&XYZviewer); CHKERRQ(ierr);
> //PetscErrorCode DMSetAuxiliaryVec(DM dm, DMLabel label, PetscInt value, Vec aux)
> DMLabel UXlabel,UYlabel, UZlabel;
> //PetscErrorCode DMLabelCreate(MPI_Comm comm, const char name[], DMLabel *label)
> ierr = DMLabelCreate(PETSC_COMM_WORLD, "X-Displacement", &UXlabel); CHKERRQ(ierr);
> ierr = DMLabelCreate(PETSC_COMM_WORLD, "Y-Displacement", &UYlabel); CHKERRQ(ierr);
> ierr = DMLabelCreate(PETSC_COMM_WORLD, "Z-Displacement", &UZlabel); CHKERRQ(ierr);
> ierr = DMSetAuxiliaryVec(dm,UXlabel, 1, UX); CHKERRQ(ierr);
> ierr = DMSetAuxiliaryVec(dm,UYlabel, 1, UY); CHKERRQ(ierr);
> ierr = DMSetAuxiliaryVec(dm,UZlabel, 1, UZ); CHKERRQ(ierr);
> //PetscErrorCode PetscViewerVTKAddField(PetscViewer viewer,PetscObject dm,PetscErrorCode (*PetscViewerVTKWriteFunction)(PetscObject,PetscViewer),PetscInt fieldnum,PetscViewerVTKFieldType fieldtype,PetscBool checkdm,PetscObject vec)
>
>
>
> //ierr = PetscViewerVTKAddField(XYZviewer, dm,PetscErrorCode (*PetscViewerVTKWriteFunction)(Vec,PetscViewer),PETSC_DEFAULT,PETSC_VTK_POINT_FIELD,PETSC_FALSE,UX);
> ierr = PetscViewerVTKAddField(XYZviewer, (PetscObject)dm,&PetscViewerVTKWriteFunction,PETSC_DEFAULT,PETSC_VTK_POINT_FIELD,PETSC_FALSE,(PetscObject)UX);
>
>
> ierr = DMPlexVTKWriteAll((PetscObject)dm, XYZviewer); CHKERRQ(ierr);
> ierr = VecAXPY(XYZ,1,U); CHKERRQ(ierr);//Add displacement field to the mesh coordinates to deform.
> ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,deformedfile,FILE_MODE_WRITE,&XYZpUviewer); CHKERRQ(ierr);
> ierr = DMPlexVTKWriteAll((PetscObject)dm, XYZpUviewer); CHKERRQ(ierr);//
> PetscViewerDestroy(&XYZviewer); PetscViewerDestroy(&XYZpUviewer);
>
> }//end else if.
> else{
> ierr = PetscPrintf(PETSC_COMM_WORLD,"No output format specified! Files not saved.\n"); CHKERRQ(ierr);
> }//end else.
>
>
> //END OUTPUT SOLUTION===========================================================
> VecDestroy(&UX); ISDestroy(&ISux);
> VecDestroy(&UY); ISDestroy(&ISuy);
> VecDestroy(&UZ); ISDestroy(&ISuz);
> //END MAX/MIN DISPLACEMENTS=====================================================
>
> //CLEANUP=====================================================================
> DMDestroy(&dm);
> KSPDestroy(&ksp);
> MatDestroy(&K); MatDestroy(&KE); MatDestroy(&matC); //MatDestroy(preKEtetra4.matB); MatDestroy(preKEtetra4.matBTCB);
> VecDestroy(&U); VecDestroy(&F);
>
> //DMLabelDestroy(&physicalgroups);//Destroyig the DM destroys the label.
> //CLEANUP=====================================================================
> //PetscErrorCode PetscMallocDump(FILE *fp)
> //ierr = PetscMallocDump(NULL);
> return PetscFinalize();//And the machine shall rest....
> }//end main.
>
> PetscErrorCode tetra4(PetscScalar* X,PetscScalar* Y, PetscScalar* Z,struct preKE *P, Mat* matC, Mat* KE){
> //INPUTS:
> //X: Global X coordinates of the elemental nodes.
> //Y: Global Y coordinates of the elemental nodes.
> //Z: Global Z coordinates of the elemental nodes.
> //J: Jacobian matrix.
> //invJ: Inverse Jacobian matrix.
> PetscErrorCode ierr;
> //For current quadrature point, get dPsi/dXi_i Xi_i = {Xi,Eta,Zeta}
> /*
> P->dPsidXi[0] = +1.; P->dPsidEta[0] = 0.0; P->dPsidZeta[0] = 0.0;
> P->dPsidXi[1] = 0.0; P->dPsidEta[1] = +1.; P->dPsidZeta[1] = 0.0;
> P->dPsidXi[2] = 0.0; P->dPsidEta[2] = 0.0; P->dPsidZeta[2] = +1.;
> P->dPsidXi[3] = -1.; P->dPsidEta[3] = -1.; P->dPsidZeta[3] = -1.;
> */
> //Populate the Jacobian matrix.
> P->J[0][0] = X[0] - X[3];
> P->J[0][1] = Y[0] - Y[3];
> P->J[0][2] = Z[0] - Z[3];
> P->J[1][0] = X[1] - X[3];
> P->J[1][1] = Y[1] - Y[3];
> P->J[1][2] = Z[1] - Z[3];
> P->J[2][0] = X[2] - X[3];
> P->J[2][1] = Y[2] - Y[3];
> P->J[2][2] = Z[2] - Z[3];
>
> //Determinant of the 3x3 Jacobian. (Expansion along 1st row).
> P->minor00 = P->J[1][1]*P->J[2][2] - P->J[2][1]*P->J[1][2];//Reuse when finding InvJ.
> P->minor01 = P->J[1][0]*P->J[2][2] - P->J[2][0]*P->J[1][2];//Reuse when finding InvJ.
> P->minor02 = P->J[1][0]*P->J[2][1] - P->J[2][0]*P->J[1][1];//Reuse when finding InvJ.
> P->detJ = P->J[0][0]*P->minor00 - P->J[0][1]*P->minor01 + P->J[0][2]*P->minor02;
> //Inverse of the 3x3 Jacobian
> P->invJ[0][0] = +P->minor00/P->detJ;//Reuse precomputed minor.
> P->invJ[0][1] = -(P->J[0][1]*P->J[2][2] - P->J[0][2]*P->J[2][1])/P->detJ;
> P->invJ[0][2] = +(P->J[0][1]*P->J[1][2] - P->J[1][1]*P->J[0][2])/P->detJ;
> P->invJ[1][0] = -P->minor01/P->detJ;//Reuse precomputed minor.
> P->invJ[1][1] = +(P->J[0][0]*P->J[2][2] - P->J[0][2]*P->J[2][0])/P->detJ;
> P->invJ[1][2] = -(P->J[0][0]*P->J[1][2] - P->J[1][0]*P->J[0][2])/P->detJ;
> P->invJ[2][0] = +P->minor02/P->detJ;//Reuse precomputed minor.
> P->invJ[2][1] = -(P->J[0][0]*P->J[2][1] - P->J[0][1]*P->J[2][0])/P->detJ;
> P->invJ[2][2] = +(P->J[0][0]*P->J[1][1] - P->J[0][1]*P->J[1][0])/P->detJ;
>
> //*****************STRAIN MATRIX (B)**************************************
> for(P->m=0;P->m<P->N;P->m++){//Scan all shape functions.
>
> P->x_in = 0 + P->m*3;//Every 3rd column starting at 0
> P->y_in = P->x_in +1;//Every 3rd column starting at 1
> P->z_in = P->y_in +1;//Every 3rd column starting at 2
>
> P->dX[0] = P->invJ[0][0]*P->dPsidXi[P->m] + P->invJ[0][1]*P->dPsidEta[P->m] + P->invJ[0][2]*P->dPsidZeta[P->m];
> P->dY[0] = P->invJ[1][0]*P->dPsidXi[P->m] + P->invJ[1][1]*P->dPsidEta[P->m] + P->invJ[1][2]*P->dPsidZeta[P->m];
> P->dZ[0] = P->invJ[2][0]*P->dPsidXi[P->m] + P->invJ[2][1]*P->dPsidEta[P->m] + P->invJ[2][2]*P->dPsidZeta[P->m];
>
> P->dX[1] = P->dZ[0]; P->dX[2] = P->dY[0];
> P->dY[1] = P->dZ[0]; P->dY[2] = P->dX[0];
> P->dZ[1] = P->dX[0]; P->dZ[2] = P->dY[0];
>
> ierr = MatSetValues(P->matB,3,P->x_insert,1,&(P->x_in),P->dX,INSERT_VALUES); CHKERRQ(ierr);
> ierr = MatSetValues(P->matB,3,P->y_insert,1,&(P->y_in),P->dY,INSERT_VALUES); CHKERRQ(ierr);
> ierr = MatSetValues(P->matB,3,P->z_insert,1,&(P->z_in),P->dZ,INSERT_VALUES); CHKERRQ(ierr);
>
> }//end "m" loop.
> ierr = MatAssemblyBegin(P->matB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> ierr = MatAssemblyEnd(P->matB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> //*****************STRAIN MATRIX (B)**************************************
>
> //Compute the matrix product B^t*C*B, scale it by the quadrature weights and add to KE.
> P->weight = -P->detJ/6;
>
> ierr = MatZeroEntries(*KE); CHKERRQ(ierr);
> ierr = MatPtAP(*matC,P->matB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&(P->matBTCB));CHKERRQ(ierr);
> ierr = MatScale(P->matBTCB,P->weight); CHKERRQ(ierr);
> ierr = MatAssemblyBegin(P->matBTCB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> ierr = MatAssemblyEnd(P->matBTCB,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> ierr = MatAXPY(*KE,1,P->matBTCB,DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);//Add contribution of current quadrature point to KE.
>
> //ierr = MatPtAP(*matC,P->matB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,KE);CHKERRQ(ierr);
> //ierr = MatScale(*KE,P->weight); CHKERRQ(ierr);
>
> ierr = MatAssemblyBegin(*KE,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> ierr = MatAssemblyEnd(*KE,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
>
> //Cleanup
> return ierr;
> }//end tetra4.
>
> PetscErrorCode ConstitutiveMatrix(Mat *matC,const char* type,PetscInt materialID){
> PetscErrorCode ierr;
> PetscBool isotropic = PETSC_FALSE,
> orthotropic = PETSC_FALSE;
> //PetscErrorCode PetscStrcmp(const char a[],const char b[],PetscBool *flg)
> ierr = PetscStrcmp(type,"isotropic",&isotropic);
> ierr = PetscStrcmp(type,"orthotropic",&orthotropic);
> ierr = MatCreate(PETSC_COMM_WORLD,matC); CHKERRQ(ierr);
> ierr = MatSetSizes(*matC,PETSC_DECIDE,PETSC_DECIDE,6,6); CHKERRQ(ierr);
> ierr = MatSetType(*matC,MATAIJ); CHKERRQ(ierr);
> ierr = MatSetUp(*matC); CHKERRQ(ierr);
>
> if(isotropic){
> PetscReal E,nu, M,L,vals[3];
> switch(materialID){
> case 0://Hardcoded properties for isotropic material #0
> E = 200;
> nu = 1./3;
> break;
> case 1://Hardcoded properties for isotropic material #1
> E = 96;
> nu = 1./3;
> break;
> }//end switch.
> M = E/(2*(1+nu)),//Lame's constant 1 ("mu").
> L = E*nu/((1+nu)*(1-2*nu));//Lame's constant 2 ("lambda").
> //PetscErrorCode MatSetValues(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv)
> PetscInt idxn[3] = {0,1,2};
> vals[0] = L+2*M; vals[1] = L; vals[2] = vals[1];
> ierr = MatSetValues(*matC,1,&idxn[0],3,idxn,vals,INSERT_VALUES); CHKERRQ(ierr);
> vals[1] = vals[0]; vals[0] = vals[2];
> ierr = MatSetValues(*matC,1,&idxn[1],3,idxn,vals,INSERT_VALUES); CHKERRQ(ierr);
> vals[2] = vals[1]; vals[1] = vals[0];
> ierr = MatSetValues(*matC,1,&idxn[2],3,idxn,vals,INSERT_VALUES); CHKERRQ(ierr);
> ierr = MatSetValue(*matC,3,3,M,INSERT_VALUES); CHKERRQ(ierr);
> ierr = MatSetValue(*matC,4,4,M,INSERT_VALUES); CHKERRQ(ierr);
> ierr = MatSetValue(*matC,5,5,M,INSERT_VALUES); CHKERRQ(ierr);
> }//end if.
> /*
> else if(orthotropic){
> switch(materialID){
> case 0:
> break;
> case 1:
> break;
> }//end switch.
> }//end else if.
> */
> ierr = MatAssemblyBegin(*matC,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> ierr = MatAssemblyEnd(*matC,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
> //MatView(*matC,0);
> return ierr;
> }//End ConstitutiveMatrix
>
> PetscErrorCode InitializeKEpreallocation(struct preKE *P,const char* type){
> PetscErrorCode ierr;
> PetscBool istetra4 = PETSC_FALSE,
> ishex8 = PETSC_FALSE;
> ierr = PetscStrcmp(type,"tetra4",&istetra4); CHKERRQ(ierr);
> ierr = PetscStrcmp(type,"hex8",&ishex8); CHKERRQ(ierr);
> if(istetra4){
> P->sizeKE = 12;
> P->N = 4;
> }//end if.
> else if(ishex8){
> P->sizeKE = 24;
> P->N = 8;
> }//end else if.
>
>
> P->x_insert[0] = 0; P->x_insert[1] = 3; P->x_insert[2] = 5;
> P->y_insert[0] = 1; P->y_insert[1] = 4; P->y_insert[2] = 5;
> P->z_insert[0] = 2; P->z_insert[1] = 3; P->z_insert[2] = 4;
> //Allocate memory for the differentiated shape function vectors.
> ierr = PetscMalloc1(P->N,&(P->dPsidXi)); CHKERRQ(ierr);
> ierr = PetscMalloc1(P->N,&(P->dPsidEta)); CHKERRQ(ierr);
> ierr = PetscMalloc1(P->N,&(P->dPsidZeta)); CHKERRQ(ierr);
>
> P->dPsidXi[0] = +1.; P->dPsidEta[0] = 0.0; P->dPsidZeta[0] = 0.0;
> P->dPsidXi[1] = 0.0; P->dPsidEta[1] = +1.; P->dPsidZeta[1] = 0.0;
> P->dPsidXi[2] = 0.0; P->dPsidEta[2] = 0.0; P->dPsidZeta[2] = +1.;
> P->dPsidXi[3] = -1.; P->dPsidEta[3] = -1.; P->dPsidZeta[3] = -1.;
>
>
> //Strain matrix.
> ierr = MatCreate(PETSC_COMM_WORLD,&(P->matB)); CHKERRQ(ierr);
> ierr = MatSetSizes(P->matB,PETSC_DECIDE,PETSC_DECIDE,6,P->sizeKE); CHKERRQ(ierr);//Hardcoded
> ierr = MatSetType(P->matB,MATAIJ); CHKERRQ(ierr);
> ierr = MatSetUp(P->matB); CHKERRQ(ierr);
>
> //Contribution matrix.
> ierr = MatCreate(PETSC_COMM_WORLD,&(P->matBTCB)); CHKERRQ(ierr);
> ierr = MatSetSizes(P->matBTCB,PETSC_DECIDE,PETSC_DECIDE,P->sizeKE,P->sizeKE); CHKERRQ(ierr);
> ierr = MatSetType(P->matBTCB,MATAIJ); CHKERRQ(ierr);
> ierr = MatSetUp(P->matBTCB); CHKERRQ(ierr);
>
> //Element stiffness matrix.
> //ierr = MatCreateSeqDense(PETSC_COMM_SELF,12,12,NULL,&KE); CHKERRQ(ierr); //PARALLEL
>
> return ierr;
> }
More information about the petsc-users
mailing list