1 static char help[] = "Time-dependent PDE in 2d for calculating joint PDF. \n"; 2 /* 3 p_t = -x_t*p_x -y_t*p_y + f(t)*p_yy 4 xmin < x < xmax, ymin < y < ymax; 5 x_t = (y - ws) y_t = (ws/2H)*(Pm - Pmax*sin(x)) 6 7 Boundary conditions: -bc_type 0 => Zero dirichlet boundary 8 -bc_type 1 => Steady state boundary condition 9 Steady state boundary condition found by setting p_t = 0 10 */ 11 12 #include <petscdm.h> 13 #include <petscdmda.h> 14 #include <petscts.h> 15 16 /* 17 User-defined data structures and routines 18 */ 19 typedef struct { 20 PetscScalar ws; /* Synchronous speed */ 21 PetscScalar H; /* Inertia constant */ 22 PetscScalar D; /* Damping constant */ 23 PetscScalar Pmax; /* Maximum power output of generator */ 24 PetscScalar PM_min; /* Mean mechanical power input */ 25 PetscScalar lambda; /* correlation time */ 26 PetscScalar q; /* noise strength */ 27 PetscScalar mux; /* Initial average angle */ 28 PetscScalar sigmax; /* Standard deviation of initial angle */ 29 PetscScalar muy; /* Average speed */ 30 PetscScalar sigmay; /* standard deviation of initial speed */ 31 PetscScalar rho; /* Cross-correlation coefficient */ 32 PetscScalar t0; /* Initial time */ 33 PetscScalar tmax; /* Final time */ 34 PetscScalar xmin; /* left boundary of angle */ 35 PetscScalar xmax; /* right boundary of angle */ 36 PetscScalar ymin; /* bottom boundary of speed */ 37 PetscScalar ymax; /* top boundary of speed */ 38 PetscScalar dx; /* x step size */ 39 PetscScalar dy; /* y step size */ 40 PetscInt bc; /* Boundary conditions */ 41 PetscScalar disper_coe; /* Dispersion coefficient */ 42 DM da; 43 } AppCtx; 44 45 PetscErrorCode Parameter_settings(AppCtx*); 46 PetscErrorCode ini_bou(Vec,AppCtx*); 47 PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*); 48 PetscErrorCode IJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*); 49 PetscErrorCode PostStep(TS); 50 51 int main(int argc, char **argv) 52 { 53 Vec x; /* Solution vector */ 54 TS ts; /* Time-stepping context */ 55 AppCtx user; /* Application context */ 56 Mat J; 57 PetscViewer viewer; 58 59 PetscFunctionBeginUser; 60 PetscCall(PetscInitialize(&argc,&argv,"petscopt_ex6", help)); 61 /* Get physics and time parameters */ 62 PetscCall(Parameter_settings(&user)); 63 /* Create a 2D DA with dof = 1 */ 64 PetscCall(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&user.da)); 65 PetscCall(DMSetFromOptions(user.da)); 66 PetscCall(DMSetUp(user.da)); 67 /* Set x and y coordinates */ 68 PetscCall(DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,0.0,1.0)); 69 70 /* Get global vector x from DM */ 71 PetscCall(DMCreateGlobalVector(user.da,&x)); 72 73 PetscCall(ini_bou(x,&user)); 74 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,"ini_x",FILE_MODE_WRITE,&viewer)); 75 PetscCall(VecView(x,viewer)); 76 PetscCall(PetscViewerDestroy(&viewer)); 77 78 /* Get Jacobian matrix structure from the da */ 79 PetscCall(DMSetMatType(user.da,MATAIJ)); 80 PetscCall(DMCreateMatrix(user.da,&J)); 81 82 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 83 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 84 PetscCall(TSSetIFunction(ts,NULL,IFunction,&user)); 85 PetscCall(TSSetIJacobian(ts,J,J,IJacobian,&user)); 86 PetscCall(TSSetApplicationContext(ts,&user)); 87 PetscCall(TSSetMaxTime(ts,user.tmax)); 88 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 89 PetscCall(TSSetTime(ts,user.t0)); 90 PetscCall(TSSetTimeStep(ts,.005)); 91 PetscCall(TSSetFromOptions(ts)); 92 PetscCall(TSSetPostStep(ts,PostStep)); 93 PetscCall(TSSolve(ts,x)); 94 95 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,"fin_x",FILE_MODE_WRITE,&viewer)); 96 PetscCall(VecView(x,viewer)); 97 PetscCall(PetscViewerDestroy(&viewer)); 98 99 PetscCall(VecDestroy(&x)); 100 PetscCall(MatDestroy(&J)); 101 PetscCall(DMDestroy(&user.da)); 102 PetscCall(TSDestroy(&ts)); 103 PetscCall(PetscFinalize()); 104 return 0; 105 } 106 107 PetscErrorCode PostStep(TS ts) 108 { 109 Vec X; 110 AppCtx *user; 111 PetscScalar sum; 112 PetscReal t; 113 114 PetscFunctionBegin; 115 PetscCall(TSGetApplicationContext(ts,&user)); 116 PetscCall(TSGetTime(ts,&t)); 117 PetscCall(TSGetSolution(ts,&X)); 118 PetscCall(VecSum(X,&sum)); 119 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"sum(p)*dw*dtheta at t = %3.2f = %3.6f\n",(double)t,(double)(sum*user->dx*user->dy))); 120 PetscFunctionReturn(0); 121 } 122 123 PetscErrorCode ini_bou(Vec X,AppCtx* user) 124 { 125 DM cda; 126 DMDACoor2d **coors; 127 PetscScalar **p; 128 Vec gc; 129 PetscInt i,j; 130 PetscInt xs,ys,xm,ym,M,N; 131 PetscScalar xi,yi; 132 PetscScalar sigmax=user->sigmax,sigmay=user->sigmay; 133 PetscScalar rho =user->rho; 134 PetscScalar mux =user->mux,muy=user->muy; 135 PetscMPIInt rank; 136 137 PetscFunctionBeginUser; 138 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 139 PetscCall(DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL)); 140 user->dx = (user->xmax - user->xmin)/(M-1); user->dy = (user->ymax - user->ymin)/(N-1); 141 PetscCall(DMGetCoordinateDM(user->da,&cda)); 142 PetscCall(DMGetCoordinates(user->da,&gc)); 143 PetscCall(DMDAVecGetArray(cda,gc,&coors)); 144 PetscCall(DMDAVecGetArray(user->da,X,&p)); 145 PetscCall(DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0)); 146 for (i=xs; i < xs+xm; i++) { 147 for (j=ys; j < ys+ym; j++) { 148 xi = coors[j][i].x; yi = coors[j][i].y; 149 if (i == 0 || j == 0 || i == M-1 || j == N-1) p[j][i] = 0.0; 150 else p[j][i] = (0.5/(PETSC_PI*sigmax*sigmay*PetscSqrtScalar(1.0-rho*rho)))*PetscExpScalar(-0.5/(1-rho*rho)*(PetscPowScalar((xi-mux)/sigmax,2) + PetscPowScalar((yi-muy)/sigmay,2) - 2*rho*(xi-mux)*(yi-muy)/(sigmax*sigmay))); 151 } 152 } 153 /* p[N/2+N%2][M/2+M%2] = 1/(user->dx*user->dy); */ 154 155 PetscCall(DMDAVecRestoreArray(cda,gc,&coors)); 156 PetscCall(DMDAVecRestoreArray(user->da,X,&p)); 157 PetscFunctionReturn(0); 158 } 159 160 /* First advection term */ 161 PetscErrorCode adv1(PetscScalar **p,PetscScalar y,PetscInt i,PetscInt j,PetscInt M,PetscScalar *p1,AppCtx *user) 162 { 163 PetscScalar f; 164 /* PetscScalar v1,v2,v3,v4,v5,s1,s2,s3; */ 165 PetscFunctionBegin; 166 /* if (i > 2 && i < M-2) { 167 v1 = (y-user->ws)*(p[j][i-2] - p[j][i-3])/user->dx; 168 v2 = (y-user->ws)*(p[j][i-1] - p[j][i-2])/user->dx; 169 v3 = (y-user->ws)*(p[j][i] - p[j][i-1])/user->dx; 170 v4 = (y-user->ws)*(p[j][i+1] - p[j][i])/user->dx; 171 v5 = (y-user->ws)*(p[j][i+1] - p[j][i+2])/user->dx; 172 173 s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3; 174 s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4; 175 s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5; 176 177 *p1 = 0.1*s1 + 0.6*s2 + 0.3*s3; 178 } else *p1 = 0.0; */ 179 f = (y - user->ws); 180 *p1 = f*(p[j][i+1] - p[j][i-1])/(2*user->dx); 181 PetscFunctionReturn(0); 182 } 183 184 /* Second advection term */ 185 PetscErrorCode adv2(PetscScalar **p,PetscScalar x,PetscInt i,PetscInt j,PetscInt N,PetscScalar *p2,AppCtx *user) 186 { 187 PetscScalar f; 188 /* PetscScalar v1,v2,v3,v4,v5,s1,s2,s3; */ 189 PetscFunctionBegin; 190 /* if (j > 2 && j < N-2) { 191 v1 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-2][i] - p[j-3][i])/user->dy; 192 v2 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-1][i] - p[j-2][i])/user->dy; 193 v3 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j][i] - p[j-1][i])/user->dy; 194 v4 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+1][i] - p[j][i])/user->dy; 195 v5 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+2][i] - p[j+1][i])/user->dy; 196 197 s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3; 198 s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4; 199 s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5; 200 201 *p2 = 0.1*s1 + 0.6*s2 + 0.3*s3; 202 } else *p2 = 0.0; */ 203 f = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(x)); 204 *p2 = f*(p[j+1][i] - p[j-1][i])/(2*user->dy); 205 PetscFunctionReturn(0); 206 } 207 208 /* Diffusion term */ 209 PetscErrorCode diffuse(PetscScalar **p,PetscInt i,PetscInt j,PetscReal t,PetscScalar *p_diff,AppCtx * user) 210 { 211 PetscFunctionBeginUser; 212 213 *p_diff = user->disper_coe*((p[j-1][i] - 2*p[j][i] + p[j+1][i])/(user->dy*user->dy)); 214 PetscFunctionReturn(0); 215 } 216 217 PetscErrorCode BoundaryConditions(PetscScalar **p,DMDACoor2d **coors,PetscInt i,PetscInt j,PetscInt M, PetscInt N,PetscScalar **f,AppCtx *user) 218 { 219 PetscScalar fwc,fthetac; 220 PetscScalar w=coors[j][i].y,theta=coors[j][i].x; 221 222 PetscFunctionBeginUser; 223 if (user->bc == 0) { /* Natural boundary condition */ 224 f[j][i] = p[j][i]; 225 } else { /* Steady state boundary condition */ 226 fthetac = user->ws/(2*user->H)*(user->PM_min - user->Pmax*PetscSinScalar(theta)); 227 fwc = (w*w/2.0 - user->ws*w); 228 if (i == 0 && j == 0) { /* left bottom corner */ 229 f[j][i] = fwc*(p[j][i+1] - p[j][i])/user->dx + fthetac*p[j][i] - user->disper_coe*(p[j+1][i] - p[j][i])/user->dy; 230 } else if (i == 0 && j == N-1) { /* right bottom corner */ 231 f[j][i] = fwc*(p[j][i+1] - p[j][i])/user->dx + fthetac*p[j][i] - user->disper_coe*(p[j][i] - p[j-1][i])/user->dy; 232 } else if (i == M-1 && j == 0) { /* left top corner */ 233 f[j][i] = fwc*(p[j][i] - p[j][i-1])/user->dx + fthetac*p[j][i] - user->disper_coe*(p[j+1][i] - p[j][i])/user->dy; 234 } else if (i == M-1 && j == N-1) { /* right top corner */ 235 f[j][i] = fwc*(p[j][i] - p[j][i-1])/user->dx + fthetac*p[j][i] - user->disper_coe*(p[j][i] - p[j-1][i])/user->dy; 236 } else if (i == 0) { /* Bottom edge */ 237 f[j][i] = fwc*(p[j][i+1] - p[j][i])/(user->dx) + fthetac*p[j][i] - user->disper_coe*(p[j+1][i] - p[j-1][i])/(2*user->dy); 238 } else if (i == M-1) { /* Top edge */ 239 f[j][i] = fwc*(p[j][i] - p[j][i-1])/(user->dx) + fthetac*p[j][i] - user->disper_coe*(p[j+1][i] - p[j-1][i])/(2*user->dy); 240 } else if (j == 0) { /* Left edge */ 241 f[j][i] = fwc*(p[j][i+1] - p[j][i-1])/(2*user->dx) + fthetac*p[j][i] - user->disper_coe*(p[j+1][i] - p[j][i])/(user->dy); 242 } else if (j == N-1) { /* Right edge */ 243 f[j][i] = fwc*(p[j][i+1] - p[j][i-1])/(2*user->dx) + fthetac*p[j][i] - user->disper_coe*(p[j][i] - p[j-1][i])/(user->dy); 244 } 245 } 246 PetscFunctionReturn(0); 247 } 248 249 PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx) 250 { 251 AppCtx *user=(AppCtx*)ctx; 252 DM cda; 253 DMDACoor2d **coors; 254 PetscScalar **p,**f,**pdot; 255 PetscInt i,j; 256 PetscInt xs,ys,xm,ym,M,N; 257 Vec localX,gc,localXdot; 258 PetscScalar p_adv1,p_adv2,p_diff; 259 260 PetscFunctionBeginUser; 261 PetscCall(DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL)); 262 PetscCall(DMGetCoordinateDM(user->da,&cda)); 263 PetscCall(DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0)); 264 265 PetscCall(DMGetLocalVector(user->da,&localX)); 266 PetscCall(DMGetLocalVector(user->da,&localXdot)); 267 268 PetscCall(DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX)); 269 PetscCall(DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX)); 270 PetscCall(DMGlobalToLocalBegin(user->da,Xdot,INSERT_VALUES,localXdot)); 271 PetscCall(DMGlobalToLocalEnd(user->da,Xdot,INSERT_VALUES,localXdot)); 272 273 PetscCall(DMGetCoordinatesLocal(user->da,&gc)); 274 275 PetscCall(DMDAVecGetArrayRead(cda,gc,&coors)); 276 PetscCall(DMDAVecGetArrayRead(user->da,localX,&p)); 277 PetscCall(DMDAVecGetArrayRead(user->da,localXdot,&pdot)); 278 PetscCall(DMDAVecGetArray(user->da,F,&f)); 279 280 user->disper_coe = PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda)); 281 for (i=xs; i < xs+xm; i++) { 282 for (j=ys; j < ys+ym; j++) { 283 if (i == 0 || j == 0 || i == M-1 || j == N-1) { 284 PetscCall(BoundaryConditions(p,coors,i,j,M,N,f,user)); 285 } else { 286 PetscCall(adv1(p,coors[j][i].y,i,j,M,&p_adv1,user)); 287 PetscCall(adv2(p,coors[j][i].x,i,j,N,&p_adv2,user)); 288 PetscCall(diffuse(p,i,j,t,&p_diff,user)); 289 f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i]; 290 } 291 } 292 } 293 PetscCall(DMDAVecRestoreArrayRead(user->da,localX,&p)); 294 PetscCall(DMDAVecRestoreArrayRead(user->da,localX,&pdot)); 295 PetscCall(DMRestoreLocalVector(user->da,&localX)); 296 PetscCall(DMRestoreLocalVector(user->da,&localXdot)); 297 PetscCall(DMDAVecRestoreArray(user->da,F,&f)); 298 PetscCall(DMDAVecRestoreArrayRead(cda,gc,&coors)); 299 300 PetscFunctionReturn(0); 301 } 302 303 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ctx) 304 { 305 AppCtx *user=(AppCtx*)ctx; 306 DM cda; 307 DMDACoor2d **coors; 308 PetscInt i,j; 309 PetscInt xs,ys,xm,ym,M,N; 310 Vec gc; 311 PetscScalar val[5],xi,yi; 312 MatStencil row,col[5]; 313 PetscScalar c1,c3,c5; 314 315 PetscFunctionBeginUser; 316 PetscCall(DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL)); 317 PetscCall(DMGetCoordinateDM(user->da,&cda)); 318 PetscCall(DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0)); 319 320 PetscCall(DMGetCoordinatesLocal(user->da,&gc)); 321 PetscCall(DMDAVecGetArrayRead(cda,gc,&coors)); 322 for (i=xs; i < xs+xm; i++) { 323 for (j=ys; j < ys+ym; j++) { 324 PetscInt nc = 0; 325 xi = coors[j][i].x; yi = coors[j][i].y; 326 row.i = i; row.j = j; 327 if (i == 0 || j == 0 || i == M-1 || j == N-1) { 328 if (user->bc == 0) { 329 col[nc].i = i; col[nc].j = j; val[nc++] = 1.0; 330 } else { 331 PetscScalar fthetac,fwc; 332 fthetac = user->ws/(2*user->H)*(user->PM_min - user->Pmax*PetscSinScalar(xi)); 333 fwc = (yi*yi/2.0 - user->ws*yi); 334 if (i==0 && j==0) { 335 col[nc].i = i+1; col[nc].j = j; val[nc++] = fwc/user->dx; 336 col[nc].i = i; col[nc].j = j+1; val[nc++] = -user->disper_coe/user->dy; 337 col[nc].i = i; col[nc].j = j; val[nc++] = -fwc/user->dx + fthetac + user->disper_coe/user->dy; 338 } else if (i==0 && j == N-1) { 339 col[nc].i = i+1; col[nc].j = j; val[nc++] = fwc/user->dx; 340 col[nc].i = i; col[nc].j = j-1; val[nc++] = user->disper_coe/user->dy; 341 col[nc].i = i; col[nc].j = j; val[nc++] = -fwc/user->dx + fthetac - user->disper_coe/user->dy; 342 } else if (i== M-1 && j == 0) { 343 col[nc].i = i-1; col[nc].j = j; val[nc++] = -fwc/user->dx; 344 col[nc].i = i; col[nc].j = j+1; val[nc++] = -user->disper_coe/user->dy; 345 col[nc].i = i; col[nc].j = j; val[nc++] = fwc/user->dx + fthetac + user->disper_coe/user->dy; 346 } else if (i == M-1 && j == N-1) { 347 col[nc].i = i-1; col[nc].j = j; val[nc++] = -fwc/user->dx; 348 col[nc].i = i; col[nc].j = j-1; val[nc++] = user->disper_coe/user->dy; 349 col[nc].i = i; col[nc].j = j; val[nc++] = fwc/user->dx + fthetac - user->disper_coe/user->dy; 350 } else if (i==0) { 351 col[nc].i = i+1; col[nc].j = j; val[nc++] = fwc/user->dx; 352 col[nc].i = i; col[nc].j = j+1; val[nc++] = -user->disper_coe/(2*user->dy); 353 col[nc].i = i; col[nc].j = j-1; val[nc++] = user->disper_coe/(2*user->dy); 354 col[nc].i = i; col[nc].j = j; val[nc++] = -fwc/user->dx + fthetac; 355 } else if (i == M-1) { 356 col[nc].i = i-1; col[nc].j = j; val[nc++] = -fwc/user->dx; 357 col[nc].i = i; col[nc].j = j+1; val[nc++] = -user->disper_coe/(2*user->dy); 358 col[nc].i = i; col[nc].j = j-1; val[nc++] = user->disper_coe/(2*user->dy); 359 col[nc].i = i; col[nc].j = j; val[nc++] = fwc/user->dx + fthetac; 360 } else if (j==0) { 361 col[nc].i = i+1; col[nc].j = j; val[nc++] = fwc/(2*user->dx); 362 col[nc].i = i-1; col[nc].j = j; val[nc++] = -fwc/(2*user->dx); 363 col[nc].i = i; col[nc].j = j+1; val[nc++] = -user->disper_coe/user->dy; 364 col[nc].i = i; col[nc].j = j; val[nc++] = user->disper_coe/user->dy + fthetac; 365 } else if (j == N-1) { 366 col[nc].i = i+1; col[nc].j = j; val[nc++] = fwc/(2*user->dx); 367 col[nc].i = i-1; col[nc].j = j; val[nc++] = -fwc/(2*user->dx); 368 col[nc].i = i; col[nc].j = j-1; val[nc++] = user->disper_coe/user->dy; 369 col[nc].i = i; col[nc].j = j; val[nc++] = -user->disper_coe/user->dy + fthetac; 370 } 371 } 372 } else { 373 c1 = (yi-user->ws)/(2*user->dx); 374 c3 = (user->ws/(2.0*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(xi))/(2*user->dy); 375 c5 = (PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda)))/(user->dy*user->dy); 376 col[nc].i = i-1; col[nc].j = j; val[nc++] = c1; 377 col[nc].i = i+1; col[nc].j = j; val[nc++] = -c1; 378 col[nc].i = i; col[nc].j = j-1; val[nc++] = c3 + c5; 379 col[nc].i = i; col[nc].j = j+1; val[nc++] = -c3 + c5; 380 col[nc].i = i; col[nc].j = j; val[nc++] = -2*c5 -a; 381 } 382 PetscCall(MatSetValuesStencil(Jpre,1,&row,nc,col,val,INSERT_VALUES)); 383 } 384 } 385 PetscCall(DMDAVecRestoreArrayRead(cda,gc,&coors)); 386 387 PetscCall(MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY)); 388 PetscCall(MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY)); 389 if (J != Jpre) { 390 PetscCall(MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY)); 391 PetscCall(MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY)); 392 } 393 PetscFunctionReturn(0); 394 } 395 396 PetscErrorCode Parameter_settings(AppCtx *user) 397 { 398 PetscBool flg; 399 400 PetscFunctionBeginUser; 401 402 /* Set default parameters */ 403 user->ws = 1.0; 404 user->H = 5.0; user->Pmax = 2.1; 405 user->PM_min = 1.0; user->lambda = 0.1; 406 user->q = 1.0; user->mux = PetscAsinScalar(user->PM_min/user->Pmax); 407 user->sigmax = 0.1; 408 user->sigmay = 0.1; user->rho = 0.0; 409 user->t0 = 0.0; user->tmax = 2.0; 410 user->xmin = -1.0; user->xmax = 10.0; 411 user->ymin = -1.0; user->ymax = 10.0; 412 user->bc = 0; 413 414 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-ws",&user->ws,&flg)); 415 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-Inertia",&user->H,&flg)); 416 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-Pmax",&user->Pmax,&flg)); 417 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-PM_min",&user->PM_min,&flg)); 418 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-lambda",&user->lambda,&flg)); 419 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-q",&user->q,&flg)); 420 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-mux",&user->mux,&flg)); 421 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-sigmax",&user->sigmax,&flg)); 422 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-muy",&user->muy,&flg)); 423 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-sigmay",&user->sigmay,&flg)); 424 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-rho",&user->rho,&flg)); 425 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-t0",&user->t0,&flg)); 426 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-tmax",&user->tmax,&flg)); 427 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-xmin",&user->xmin,&flg)); 428 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-xmax",&user->xmax,&flg)); 429 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-ymin",&user->ymin,&flg)); 430 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-ymax",&user->ymax,&flg)); 431 PetscCall(PetscOptionsGetInt(NULL,NULL,"-bc_type",&user->bc,&flg)); 432 user->muy = user->ws; 433 PetscFunctionReturn(0); 434 } 435 436 /*TEST 437 438 build: 439 requires: !complex 440 441 test: 442 args: -nox -ts_max_steps 2 443 localrunfiles: petscopt_ex6 444 timeoutfactor: 4 445 446 TEST*/ 447