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