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