1 #include "wash.h" 2 3 /* Subroutines for Pipe */ 4 /* -------------------------------------------------------*/ 5 6 /* 7 PipeCreate - Create Pipe object. 8 9 Input Parameters: 10 comm - MPI communicator 11 12 Output Parameter: 13 . pipe - location to put the PIPE context 14 */ 15 PetscErrorCode PipeCreate(MPI_Comm comm,Pipe *pipe) 16 { 17 PetscErrorCode ierr; 18 19 PetscFunctionBegin; 20 ierr = PetscNew(pipe);CHKERRQ(ierr); 21 PetscFunctionReturn(0); 22 } 23 24 /* 25 PipeDestroy - Destroy Pipe object. 26 27 Input Parameters: 28 pipe - Reference to pipe intended to be destroyed. 29 */ 30 PetscErrorCode PipeDestroy(Pipe *pipe) 31 { 32 PetscErrorCode ierr; 33 34 PetscFunctionBegin; 35 if (!*pipe) PetscFunctionReturn(0); 36 37 ierr = PipeDestroyJacobian(*pipe);CHKERRQ(ierr); 38 ierr = VecDestroy(&(*pipe)->x);CHKERRQ(ierr); 39 ierr = DMDestroy(&(*pipe)->da);CHKERRQ(ierr); 40 PetscFunctionReturn(0); 41 } 42 43 /* 44 PipeSetParameters - Set parameters for Pipe context 45 46 Input Parameter: 47 + pipe - PIPE object 48 . length - 49 . nnodes - 50 . D - 51 . a - 52 - fric - 53 */ 54 PetscErrorCode PipeSetParameters(Pipe pipe,PetscReal length,PetscReal D,PetscReal a,PetscReal fric) 55 { 56 PetscFunctionBegin; 57 pipe->length = length; 58 pipe->D = D; 59 pipe->a = a; 60 pipe->fric = fric; 61 PetscFunctionReturn(0); 62 } 63 64 /* 65 PipeSetUp - Set up pipe based on set parameters. 66 */ 67 PetscErrorCode PipeSetUp(Pipe pipe) 68 { 69 DMDALocalInfo info; 70 PetscErrorCode ierr; 71 72 PetscFunctionBegin; 73 ierr = DMDACreate1d(PETSC_COMM_SELF, DM_BOUNDARY_GHOSTED, pipe->nnodes, 2, 1, NULL, &pipe->da);CHKERRQ(ierr); 74 ierr = DMSetFromOptions(pipe->da);CHKERRQ(ierr); 75 ierr = DMSetUp(pipe->da);CHKERRQ(ierr); 76 ierr = DMDASetFieldName(pipe->da, 0, "Q");CHKERRQ(ierr); 77 ierr = DMDASetFieldName(pipe->da, 1, "H");CHKERRQ(ierr); 78 ierr = DMDASetUniformCoordinates(pipe->da, 0, pipe->length, 0, 0, 0, 0);CHKERRQ(ierr); 79 ierr = DMCreateGlobalVector(pipe->da, &(pipe->x));CHKERRQ(ierr); 80 81 ierr = DMDAGetLocalInfo(pipe->da, &info);CHKERRQ(ierr); 82 83 pipe->rad = pipe->D / 2; 84 pipe->A = PETSC_PI*pipe->rad*pipe->rad; 85 pipe->R = pipe->fric / (2*pipe->D*pipe->A); 86 PetscFunctionReturn(0); 87 } 88 89 /* 90 PipeCreateJacobian - Create Jacobian matrix structures for a Pipe. 91 92 Collective on Pipe 93 94 Input Parameter: 95 + pipe - the Pipe object 96 - Jin - array of three constructed Jacobian matrices to be reused. Set NULL if it is not available 97 98 Output Parameter: 99 . J - array of three empty Jacobian matrices 100 101 Level: beginner 102 */ 103 PetscErrorCode PipeCreateJacobian(Pipe pipe,Mat *Jin,Mat *J[]) 104 { 105 PetscErrorCode ierr; 106 Mat *Jpipe; 107 PetscInt M,rows[2],cols[2],*nz; 108 PetscScalar *aa; 109 110 PetscFunctionBegin; 111 if (Jin) { 112 *J = Jin; 113 pipe->jacobian = Jin; 114 ierr = PetscObjectReference((PetscObject)(Jin[0]));CHKERRQ(ierr); 115 PetscFunctionReturn(0); 116 } 117 ierr = PetscMalloc1(3,&Jpipe);CHKERRQ(ierr); 118 119 /* Jacobian for this pipe */ 120 ierr = DMSetMatrixStructureOnly(pipe->da,PETSC_TRUE);CHKERRQ(ierr); 121 ierr = DMCreateMatrix(pipe->da,&Jpipe[0]);CHKERRQ(ierr); 122 ierr = DMSetMatrixStructureOnly(pipe->da,PETSC_FALSE);CHKERRQ(ierr); 123 124 /* Jacobian for upstream vertex */ 125 ierr = MatGetSize(Jpipe[0],&M,NULL);CHKERRQ(ierr); 126 ierr = PetscCalloc2(M,&nz,4,&aa);CHKERRQ(ierr); 127 128 ierr = MatCreate(PETSC_COMM_SELF,&Jpipe[1]);CHKERRQ(ierr); 129 ierr = MatSetSizes(Jpipe[1],PETSC_DECIDE,PETSC_DECIDE,M,2);CHKERRQ(ierr); 130 ierr = MatSetFromOptions(Jpipe[1]);CHKERRQ(ierr); 131 ierr = MatSetOption(Jpipe[1],MAT_STRUCTURE_ONLY,PETSC_TRUE);CHKERRQ(ierr); 132 nz[0] = 2; nz[1] = 2; 133 rows[0] = 0; rows[1] = 1; 134 cols[0] = 0; cols[1] = 1; 135 ierr = MatSeqAIJSetPreallocation(Jpipe[1],0,nz);CHKERRQ(ierr); 136 ierr = MatSetValues(Jpipe[1],2,rows,2,cols,aa,INSERT_VALUES);CHKERRQ(ierr); 137 ierr = MatAssemblyBegin(Jpipe[1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 138 ierr = MatAssemblyEnd(Jpipe[1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 139 140 /* Jacobian for downstream vertex */ 141 ierr = MatCreate(PETSC_COMM_SELF,&Jpipe[2]);CHKERRQ(ierr); 142 ierr = MatSetSizes(Jpipe[2],PETSC_DECIDE,PETSC_DECIDE,M,2);CHKERRQ(ierr); 143 ierr = MatSetFromOptions(Jpipe[2]);CHKERRQ(ierr); 144 ierr = MatSetOption(Jpipe[2],MAT_STRUCTURE_ONLY,PETSC_TRUE);CHKERRQ(ierr); 145 nz[0] = 0; nz[1] = 0; nz[M-2] = 2; nz[M-1] = 2; 146 rows[0] = M - 2; rows[1] = M - 1; 147 ierr = MatSeqAIJSetPreallocation(Jpipe[2],0,nz);CHKERRQ(ierr); 148 ierr = MatSetValues(Jpipe[2],2,rows,2,cols,aa,INSERT_VALUES);CHKERRQ(ierr); 149 ierr = MatAssemblyBegin(Jpipe[2],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 150 ierr = MatAssemblyEnd(Jpipe[2],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 151 152 ierr = PetscFree2(nz,aa);CHKERRQ(ierr); 153 154 *J = Jpipe; 155 pipe->jacobian = Jpipe; 156 PetscFunctionReturn(0); 157 } 158 159 PetscErrorCode PipeDestroyJacobian(Pipe pipe) 160 { 161 PetscErrorCode ierr; 162 Mat *Jpipe = pipe->jacobian; 163 PetscInt i; 164 165 PetscFunctionBegin; 166 if (Jpipe) { 167 for (i=0; i<3; i++) { 168 ierr = MatDestroy(&Jpipe[i]);CHKERRQ(ierr); 169 } 170 } 171 ierr = PetscFree(Jpipe);CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 175 /* 176 JunctionCreateJacobian - Create Jacobian matrices for a vertex. 177 178 Collective on Pipe 179 180 Input Parameter: 181 + dm - the DMNetwork object 182 . v - vertex point 183 - Jin - Jacobian patterns created by JunctionCreateJacobianSample() for reuse 184 185 Output Parameter: 186 . J - array of Jacobian matrices (see dmnetworkimpl.h) 187 188 Level: beginner 189 */ 190 PetscErrorCode JunctionCreateJacobian(DM dm,PetscInt v,Mat *Jin,Mat *J[]) 191 { 192 PetscErrorCode ierr; 193 Mat *Jv; 194 PetscInt nedges,e,i,M,N,*rows,*cols; 195 PetscBool isSelf; 196 const PetscInt *edges,*cone; 197 PetscScalar *zeros; 198 199 PetscFunctionBegin; 200 /* Get array size of Jv */ 201 ierr = DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);CHKERRQ(ierr); 202 if (nedges <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"%d vertex, nedges %d\n",v,nedges); 203 204 /* two Jacobians for each connected edge: J(v,e) and J(v,vc); adding J(v,v), total 2*nedges+1 Jacobians */ 205 ierr = PetscCalloc1(2*nedges+1,&Jv);CHKERRQ(ierr); 206 207 /* Create dense zero block for this vertex: J[0] = Jacobian(v,v) */ 208 ierr = DMNetworkGetComponent(dm,v,-1,NULL,NULL,&M);CHKERRQ(ierr); 209 if (M !=2) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"M != 2",M); 210 ierr = PetscMalloc3(M,&rows,M,&cols,M*M,&zeros);CHKERRQ(ierr); 211 ierr = PetscArrayzero(zeros,M*M);CHKERRQ(ierr); 212 for (i=0; i<M; i++) rows[i] = i; 213 214 for (e=0; e<nedges; e++) { 215 /* create Jv[2*e+1] = Jacobian(v,e), e: supporting edge */ 216 ierr = DMNetworkGetConnectedVertices(dm,edges[e],&cone);CHKERRQ(ierr); 217 isSelf = (v == cone[0]) ? PETSC_TRUE:PETSC_FALSE; 218 219 if (Jin) { 220 if (isSelf) { 221 Jv[2*e+1] = Jin[0]; 222 } else { 223 Jv[2*e+1] = Jin[1]; 224 } 225 Jv[2*e+2] = Jin[2]; 226 ierr = PetscObjectReference((PetscObject)(Jv[2*e+1]));CHKERRQ(ierr); 227 ierr = PetscObjectReference((PetscObject)(Jv[2*e+2]));CHKERRQ(ierr); 228 } else { 229 /* create J(v,e) */ 230 ierr = MatCreate(PETSC_COMM_SELF,&Jv[2*e+1]);CHKERRQ(ierr); 231 ierr = DMNetworkGetComponent(dm,edges[e],-1,NULL,NULL,&N);CHKERRQ(ierr); 232 ierr = MatSetSizes(Jv[2*e+1],PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 233 ierr = MatSetFromOptions(Jv[2*e+1]);CHKERRQ(ierr); 234 ierr = MatSetOption(Jv[2*e+1],MAT_STRUCTURE_ONLY,PETSC_TRUE);CHKERRQ(ierr); 235 ierr = MatSeqAIJSetPreallocation(Jv[2*e+1],2,NULL);CHKERRQ(ierr); 236 if (N) { 237 if (isSelf) { /* coupling at upstream */ 238 for (i=0; i<2; i++) cols[i] = i; 239 } else { /* coupling at downstream */ 240 cols[0] = N-2; cols[1] = N-1; 241 } 242 ierr = MatSetValues(Jv[2*e+1],2,rows,2,cols,zeros,INSERT_VALUES);CHKERRQ(ierr); 243 } 244 ierr = MatAssemblyBegin(Jv[2*e+1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 245 ierr = MatAssemblyEnd(Jv[2*e+1],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 246 247 /* create Jv[2*e+2] = Jacobian(v,vc), vc: connected vertex. 248 In WashNetwork, v and vc are not connected, thus Jacobian(v,vc) is empty */ 249 ierr = MatCreate(PETSC_COMM_SELF,&Jv[2*e+2]);CHKERRQ(ierr); 250 ierr = MatSetSizes(Jv[2*e+2],PETSC_DECIDE,PETSC_DECIDE,M,M);CHKERRQ(ierr); /* empty matrix, sizes can be arbitrary */ 251 ierr = MatSetFromOptions(Jv[2*e+2]);CHKERRQ(ierr); 252 ierr = MatSetOption(Jv[2*e+2],MAT_STRUCTURE_ONLY,PETSC_TRUE);CHKERRQ(ierr); 253 ierr = MatSeqAIJSetPreallocation(Jv[2*e+2],1,NULL);CHKERRQ(ierr); 254 ierr = MatAssemblyBegin(Jv[2*e+2],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 255 ierr = MatAssemblyEnd(Jv[2*e+2],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 256 } 257 } 258 ierr = PetscFree3(rows,cols,zeros);CHKERRQ(ierr); 259 260 *J = Jv; 261 PetscFunctionReturn(0); 262 } 263 264 PetscErrorCode JunctionDestroyJacobian(DM dm,PetscInt v,Junction junc) 265 { 266 PetscErrorCode ierr; 267 Mat *Jv=junc->jacobian; 268 const PetscInt *edges; 269 PetscInt nedges,e; 270 271 PetscFunctionBegin; 272 if (!Jv) PetscFunctionReturn(0); 273 274 ierr = DMNetworkGetSupportingEdges(dm,v,&nedges,&edges);CHKERRQ(ierr); 275 for (e=0; e<nedges; e++) { 276 ierr = MatDestroy(&Jv[2*e+1]);CHKERRQ(ierr); 277 ierr = MatDestroy(&Jv[2*e+2]);CHKERRQ(ierr); 278 } 279 ierr = PetscFree(Jv);CHKERRQ(ierr); 280 PetscFunctionReturn(0); 281 } 282