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; 127 nz[1] = 2; 128 rows[0] = 0; 129 rows[1] = 1; 130 cols[0] = 0; 131 cols[1] = 1; 132 PetscCall(MatSeqAIJSetPreallocation(Jpipe[1], 0, nz)); 133 PetscCall(MatSetValues(Jpipe[1], 2, rows, 2, cols, aa, INSERT_VALUES)); 134 PetscCall(MatAssemblyBegin(Jpipe[1], MAT_FINAL_ASSEMBLY)); 135 PetscCall(MatAssemblyEnd(Jpipe[1], MAT_FINAL_ASSEMBLY)); 136 137 /* Jacobian for downstream vertex */ 138 PetscCall(MatCreate(PETSC_COMM_SELF, &Jpipe[2])); 139 PetscCall(MatSetSizes(Jpipe[2], PETSC_DECIDE, PETSC_DECIDE, M, 2)); 140 PetscCall(MatSetFromOptions(Jpipe[2])); 141 PetscCall(MatSetOption(Jpipe[2], MAT_STRUCTURE_ONLY, PETSC_TRUE)); 142 nz[0] = 0; 143 nz[1] = 0; 144 nz[M - 2] = 2; 145 nz[M - 1] = 2; 146 rows[0] = M - 2; 147 rows[1] = M - 1; 148 PetscCall(MatSeqAIJSetPreallocation(Jpipe[2], 0, nz)); 149 PetscCall(MatSetValues(Jpipe[2], 2, rows, 2, cols, aa, INSERT_VALUES)); 150 PetscCall(MatAssemblyBegin(Jpipe[2], MAT_FINAL_ASSEMBLY)); 151 PetscCall(MatAssemblyEnd(Jpipe[2], MAT_FINAL_ASSEMBLY)); 152 153 PetscCall(PetscFree2(nz, aa)); 154 155 *J = Jpipe; 156 pipe->jacobian = Jpipe; 157 PetscFunctionReturn(0); 158 } 159 160 PetscErrorCode PipeDestroyJacobian(Pipe pipe) 161 { 162 Mat *Jpipe = pipe->jacobian; 163 PetscInt i; 164 165 PetscFunctionBegin; 166 if (Jpipe) { 167 for (i = 0; i < 3; i++) PetscCall(MatDestroy(&Jpipe[i])); 168 } 169 PetscCall(PetscFree(Jpipe)); 170 PetscFunctionReturn(0); 171 } 172 173 /* 174 JunctionCreateJacobian - Create Jacobian matrices for a vertex. 175 176 Collective on Pipe 177 178 Input Parameter: 179 + dm - the DMNetwork object 180 . v - vertex point 181 - Jin - Jacobian patterns created by JunctionCreateJacobianSample() for reuse 182 183 Output Parameter: 184 . J - array of Jacobian matrices (see dmnetworkimpl.h) 185 186 Level: beginner 187 */ 188 PetscErrorCode JunctionCreateJacobian(DM dm, PetscInt v, Mat *Jin, Mat *J[]) 189 { 190 Mat *Jv; 191 PetscInt nedges, e, i, M, N, *rows, *cols; 192 PetscBool isSelf; 193 const PetscInt *edges, *cone; 194 PetscScalar *zeros; 195 196 PetscFunctionBegin; 197 /* Get array size of Jv */ 198 PetscCall(DMNetworkGetSupportingEdges(dm, v, &nedges, &edges)); 199 PetscCheck(nedges > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "%" PetscInt_FMT " vertex, nedges %" PetscInt_FMT, v, nedges); 200 201 /* two Jacobians for each connected edge: J(v,e) and J(v,vc); adding J(v,v), total 2*nedges+1 Jacobians */ 202 PetscCall(PetscCalloc1(2 * nedges + 1, &Jv)); 203 204 /* Create dense zero block for this vertex: J[0] = Jacobian(v,v) */ 205 PetscCall(DMNetworkGetComponent(dm, v, -1, NULL, NULL, &M)); 206 PetscCheck(M == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "%" PetscInt_FMT " != 2", M); 207 PetscCall(PetscMalloc3(M, &rows, M, &cols, M * M, &zeros)); 208 PetscCall(PetscArrayzero(zeros, M * M)); 209 for (i = 0; i < M; i++) rows[i] = i; 210 211 for (e = 0; e < nedges; e++) { 212 /* create Jv[2*e+1] = Jacobian(v,e), e: supporting edge */ 213 PetscCall(DMNetworkGetConnectedVertices(dm, edges[e], &cone)); 214 isSelf = (v == cone[0]) ? PETSC_TRUE : PETSC_FALSE; 215 216 if (Jin) { 217 if (isSelf) { 218 Jv[2 * e + 1] = Jin[0]; 219 } else { 220 Jv[2 * e + 1] = Jin[1]; 221 } 222 Jv[2 * e + 2] = Jin[2]; 223 PetscCall(PetscObjectReference((PetscObject)(Jv[2 * e + 1]))); 224 PetscCall(PetscObjectReference((PetscObject)(Jv[2 * e + 2]))); 225 } else { 226 /* create J(v,e) */ 227 PetscCall(MatCreate(PETSC_COMM_SELF, &Jv[2 * e + 1])); 228 PetscCall(DMNetworkGetComponent(dm, edges[e], -1, NULL, NULL, &N)); 229 PetscCall(MatSetSizes(Jv[2 * e + 1], PETSC_DECIDE, PETSC_DECIDE, M, N)); 230 PetscCall(MatSetFromOptions(Jv[2 * e + 1])); 231 PetscCall(MatSetOption(Jv[2 * e + 1], MAT_STRUCTURE_ONLY, PETSC_TRUE)); 232 PetscCall(MatSeqAIJSetPreallocation(Jv[2 * e + 1], 2, NULL)); 233 if (N) { 234 if (isSelf) { /* coupling at upstream */ 235 for (i = 0; i < 2; i++) cols[i] = i; 236 } else { /* coupling at downstream */ 237 cols[0] = N - 2; 238 cols[1] = N - 1; 239 } 240 PetscCall(MatSetValues(Jv[2 * e + 1], 2, rows, 2, cols, zeros, INSERT_VALUES)); 241 } 242 PetscCall(MatAssemblyBegin(Jv[2 * e + 1], MAT_FINAL_ASSEMBLY)); 243 PetscCall(MatAssemblyEnd(Jv[2 * e + 1], MAT_FINAL_ASSEMBLY)); 244 245 /* create Jv[2*e+2] = Jacobian(v,vc), vc: connected vertex. 246 In WashNetwork, v and vc are not connected, thus Jacobian(v,vc) is empty */ 247 PetscCall(MatCreate(PETSC_COMM_SELF, &Jv[2 * e + 2])); 248 PetscCall(MatSetSizes(Jv[2 * e + 2], PETSC_DECIDE, PETSC_DECIDE, M, M)); /* empty matrix, sizes can be arbitrary */ 249 PetscCall(MatSetFromOptions(Jv[2 * e + 2])); 250 PetscCall(MatSetOption(Jv[2 * e + 2], MAT_STRUCTURE_ONLY, PETSC_TRUE)); 251 PetscCall(MatSeqAIJSetPreallocation(Jv[2 * e + 2], 1, NULL)); 252 PetscCall(MatAssemblyBegin(Jv[2 * e + 2], MAT_FINAL_ASSEMBLY)); 253 PetscCall(MatAssemblyEnd(Jv[2 * e + 2], MAT_FINAL_ASSEMBLY)); 254 } 255 } 256 PetscCall(PetscFree3(rows, cols, zeros)); 257 258 *J = Jv; 259 PetscFunctionReturn(0); 260 } 261 262 PetscErrorCode JunctionDestroyJacobian(DM dm, PetscInt v, Junction junc) 263 { 264 Mat *Jv = junc->jacobian; 265 const PetscInt *edges; 266 PetscInt nedges, e; 267 268 PetscFunctionBegin; 269 if (!Jv) PetscFunctionReturn(0); 270 271 PetscCall(DMNetworkGetSupportingEdges(dm, v, &nedges, &edges)); 272 for (e = 0; e < nedges; e++) { 273 PetscCall(MatDestroy(&Jv[2 * e + 1])); 274 PetscCall(MatDestroy(&Jv[2 * e + 2])); 275 } 276 PetscCall(PetscFree(Jv)); 277 PetscFunctionReturn(0); 278 } 279