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