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