1 2 static char help[] ="Solves one dimensional Burger's equation compares with exact solution\n\n"; 3 4 /* 5 6 Not yet tested in parallel 7 8 */ 9 /* 10 Concepts: TS^time-dependent nonlinear problems 11 Concepts: TS^Burger's equation 12 Processors: n 13 */ 14 15 /* ------------------------------------------------------------------------ 16 17 This program uses the one-dimensional Burger's equation 18 u_t = mu*u_xx - u u_x, 19 on the domain 0 <= x <= 1, with periodic boundary conditions 20 21 The operators are discretized with the spectral element method 22 23 See the paper PDE-CONSTRAINED OPTIMIZATION WITH SPECTRAL ELEMENTS USING PETSC AND TAO 24 by OANA MARIN, EMIL CONSTANTINESCU, AND BARRY SMITH for details on the exact solution 25 used 26 27 See src/tao/unconstrained/tutorials/burgers_spectral.c 28 29 ------------------------------------------------------------------------- */ 30 31 #include <petscts.h> 32 #include <petscdt.h> 33 #include <petscdraw.h> 34 #include <petscdmda.h> 35 36 /* 37 User-defined application context - contains data needed by the 38 application-provided call-back routines. 39 */ 40 41 typedef struct { 42 PetscInt n; /* number of nodes */ 43 PetscReal *nodes; /* GLL nodes */ 44 PetscReal *weights; /* GLL weights */ 45 } PetscGLL; 46 47 typedef struct { 48 PetscInt N; /* grid points per elements*/ 49 PetscInt E; /* number of elements */ 50 PetscReal tol_L2,tol_max; /* error norms */ 51 PetscInt steps; /* number of timesteps */ 52 PetscReal Tend; /* endtime */ 53 PetscReal mu; /* viscosity */ 54 PetscReal L; /* total length of domain */ 55 PetscReal Le; 56 PetscReal Tadj; 57 } PetscParam; 58 59 typedef struct { 60 Vec grid; /* total grid */ 61 Vec curr_sol; 62 } PetscData; 63 64 typedef struct { 65 Vec grid; /* total grid */ 66 Vec mass; /* mass matrix for total integration */ 67 Mat stiff; /* stifness matrix */ 68 Mat keptstiff; 69 Mat grad; 70 PetscGLL gll; 71 } PetscSEMOperators; 72 73 typedef struct { 74 DM da; /* distributed array data structure */ 75 PetscSEMOperators SEMop; 76 PetscParam param; 77 PetscData dat; 78 TS ts; 79 PetscReal initial_dt; 80 } AppCtx; 81 82 /* 83 User-defined routines 84 */ 85 extern PetscErrorCode RHSMatrixLaplaciangllDM(TS,PetscReal,Vec,Mat,Mat,void*); 86 extern PetscErrorCode RHSMatrixAdvectiongllDM(TS,PetscReal,Vec,Mat,Mat,void*); 87 extern PetscErrorCode TrueSolution(TS,PetscReal,Vec,AppCtx*); 88 extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*); 89 extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*); 90 91 int main(int argc,char **argv) 92 { 93 AppCtx appctx; /* user-defined application context */ 94 PetscInt i, xs, xm, ind, j, lenglob; 95 PetscReal x, *wrk_ptr1, *wrk_ptr2; 96 MatNullSpace nsp; 97 PetscMPIInt size; 98 99 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 100 Initialize program and set problem parameters 101 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 102 PetscFunctionBegin; 103 104 CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 105 106 /*initialize parameters */ 107 appctx.param.N = 10; /* order of the spectral element */ 108 appctx.param.E = 10; /* number of elements */ 109 appctx.param.L = 4.0; /* length of the domain */ 110 appctx.param.mu = 0.01; /* diffusion coefficient */ 111 appctx.initial_dt = 5e-3; 112 appctx.param.steps = PETSC_MAX_INT; 113 appctx.param.Tend = 4; 114 115 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-N",&appctx.param.N,NULL)); 116 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-E",&appctx.param.E,NULL)); 117 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-Tend",&appctx.param.Tend,NULL)); 118 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-mu",&appctx.param.mu,NULL)); 119 appctx.param.Le = appctx.param.L/appctx.param.E; 120 121 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 122 PetscCheck((appctx.param.E % size) == 0,PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of elements must be divisible by number of processes"); 123 124 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 125 Create GLL data structures 126 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 127 CHKERRQ(PetscMalloc2(appctx.param.N,&appctx.SEMop.gll.nodes,appctx.param.N,&appctx.SEMop.gll.weights)); 128 CHKERRQ(PetscDTGaussLobattoLegendreQuadrature(appctx.param.N,PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA,appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights)); 129 appctx.SEMop.gll.n = appctx.param.N; 130 lenglob = appctx.param.E*(appctx.param.N-1); 131 132 /* 133 Create distributed array (DMDA) to manage parallel grid and vectors 134 and to set up the ghost point communication pattern. There are E*(Nl-1)+1 135 total grid values spread equally among all the processors, except first and last 136 */ 137 138 CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,lenglob,1,1,NULL,&appctx.da)); 139 CHKERRQ(DMSetFromOptions(appctx.da)); 140 CHKERRQ(DMSetUp(appctx.da)); 141 142 /* 143 Extract global and local vectors from DMDA; we use these to store the 144 approximate solution. Then duplicate these for remaining vectors that 145 have the same types. 146 */ 147 148 CHKERRQ(DMCreateGlobalVector(appctx.da,&appctx.dat.curr_sol)); 149 CHKERRQ(VecDuplicate(appctx.dat.curr_sol,&appctx.SEMop.grid)); 150 CHKERRQ(VecDuplicate(appctx.dat.curr_sol,&appctx.SEMop.mass)); 151 152 CHKERRQ(DMDAGetCorners(appctx.da,&xs,NULL,NULL,&xm,NULL,NULL)); 153 CHKERRQ(DMDAVecGetArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1)); 154 CHKERRQ(DMDAVecGetArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2)); 155 156 /* Compute function over the locally owned part of the grid */ 157 158 xs=xs/(appctx.param.N-1); 159 xm=xm/(appctx.param.N-1); 160 161 /* 162 Build total grid and mass over entire mesh (multi-elemental) 163 */ 164 165 for (i=xs; i<xs+xm; i++) { 166 for (j=0; j<appctx.param.N-1; j++) { 167 x = (appctx.param.Le/2.0)*(appctx.SEMop.gll.nodes[j]+1.0)+appctx.param.Le*i; 168 ind=i*(appctx.param.N-1)+j; 169 wrk_ptr1[ind]=x; 170 wrk_ptr2[ind]=.5*appctx.param.Le*appctx.SEMop.gll.weights[j]; 171 if (j==0) wrk_ptr2[ind]+=.5*appctx.param.Le*appctx.SEMop.gll.weights[j]; 172 } 173 } 174 CHKERRQ(DMDAVecRestoreArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1)); 175 CHKERRQ(DMDAVecRestoreArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2)); 176 177 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 178 Create matrix data structure; set matrix evaluation routine. 179 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 180 CHKERRQ(DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE)); 181 CHKERRQ(DMCreateMatrix(appctx.da,&appctx.SEMop.stiff)); 182 CHKERRQ(DMCreateMatrix(appctx.da,&appctx.SEMop.grad)); 183 /* 184 For linear problems with a time-dependent f(u,t) in the equation 185 u_t = f(u,t), the user provides the discretized right-hand-side 186 as a time-dependent matrix. 187 */ 188 CHKERRQ(RHSMatrixLaplaciangllDM(appctx.ts,0.0,appctx.dat.curr_sol,appctx.SEMop.stiff,appctx.SEMop.stiff,&appctx)); 189 CHKERRQ(RHSMatrixAdvectiongllDM(appctx.ts,0.0,appctx.dat.curr_sol,appctx.SEMop.grad,appctx.SEMop.grad,&appctx)); 190 /* 191 For linear problems with a time-dependent f(u,t) in the equation 192 u_t = f(u,t), the user provides the discretized right-hand-side 193 as a time-dependent matrix. 194 */ 195 196 CHKERRQ(MatDuplicate(appctx.SEMop.stiff,MAT_COPY_VALUES,&appctx.SEMop.keptstiff)); 197 198 /* attach the null space to the matrix, this probably is not needed but does no harm */ 199 CHKERRQ(MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp)); 200 CHKERRQ(MatSetNullSpace(appctx.SEMop.stiff,nsp)); 201 CHKERRQ(MatSetNullSpace(appctx.SEMop.keptstiff,nsp)); 202 CHKERRQ(MatNullSpaceTest(nsp,appctx.SEMop.stiff,NULL)); 203 CHKERRQ(MatNullSpaceDestroy(&nsp)); 204 /* attach the null space to the matrix, this probably is not needed but does no harm */ 205 CHKERRQ(MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp)); 206 CHKERRQ(MatSetNullSpace(appctx.SEMop.grad,nsp)); 207 CHKERRQ(MatNullSpaceTest(nsp,appctx.SEMop.grad,NULL)); 208 CHKERRQ(MatNullSpaceDestroy(&nsp)); 209 210 /* Create the TS solver that solves the ODE and its adjoint; set its options */ 211 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&appctx.ts)); 212 CHKERRQ(TSSetProblemType(appctx.ts,TS_NONLINEAR)); 213 CHKERRQ(TSSetType(appctx.ts,TSRK)); 214 CHKERRQ(TSSetDM(appctx.ts,appctx.da)); 215 CHKERRQ(TSSetTime(appctx.ts,0.0)); 216 CHKERRQ(TSSetTimeStep(appctx.ts,appctx.initial_dt)); 217 CHKERRQ(TSSetMaxSteps(appctx.ts,appctx.param.steps)); 218 CHKERRQ(TSSetMaxTime(appctx.ts,appctx.param.Tend)); 219 CHKERRQ(TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP)); 220 CHKERRQ(TSSetTolerances(appctx.ts,1e-7,NULL,1e-7,NULL)); 221 CHKERRQ(TSSetSaveTrajectory(appctx.ts)); 222 CHKERRQ(TSSetFromOptions(appctx.ts)); 223 CHKERRQ(TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx)); 224 CHKERRQ(TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,RHSJacobian,&appctx)); 225 226 /* Set Initial conditions for the problem */ 227 CHKERRQ(TrueSolution(appctx.ts,0,appctx.dat.curr_sol,&appctx)); 228 229 CHKERRQ(TSSetSolutionFunction(appctx.ts,(PetscErrorCode (*)(TS,PetscReal,Vec,void *))TrueSolution,&appctx)); 230 CHKERRQ(TSSetTime(appctx.ts,0.0)); 231 CHKERRQ(TSSetStepNumber(appctx.ts,0)); 232 233 CHKERRQ(TSSolve(appctx.ts,appctx.dat.curr_sol)); 234 235 CHKERRQ(MatDestroy(&appctx.SEMop.stiff)); 236 CHKERRQ(MatDestroy(&appctx.SEMop.keptstiff)); 237 CHKERRQ(MatDestroy(&appctx.SEMop.grad)); 238 CHKERRQ(VecDestroy(&appctx.SEMop.grid)); 239 CHKERRQ(VecDestroy(&appctx.SEMop.mass)); 240 CHKERRQ(VecDestroy(&appctx.dat.curr_sol)); 241 CHKERRQ(PetscFree2(appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights)); 242 CHKERRQ(DMDestroy(&appctx.da)); 243 CHKERRQ(TSDestroy(&appctx.ts)); 244 245 /* 246 Always call PetscFinalize() before exiting a program. This routine 247 - finalizes the PETSc libraries as well as MPI 248 - provides summary and diagnostic information if certain runtime 249 options are chosen (e.g., -log_summary). 250 */ 251 CHKERRQ(PetscFinalize()); 252 return 0; 253 } 254 255 /* 256 TrueSolution() computes the true solution for the PDE 257 258 Input Parameter: 259 u - uninitialized solution vector (global) 260 appctx - user-defined application context 261 262 Output Parameter: 263 u - vector with solution at initial time (global) 264 */ 265 PetscErrorCode TrueSolution(TS ts, PetscReal t, Vec u,AppCtx *appctx) 266 { 267 PetscScalar *s; 268 const PetscScalar *xg; 269 PetscInt i,xs,xn; 270 271 CHKERRQ(DMDAVecGetArray(appctx->da,u,&s)); 272 CHKERRQ(DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg)); 273 CHKERRQ(DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL)); 274 for (i=xs; i<xs+xn; i++) { 275 s[i]=2.0*appctx->param.mu*PETSC_PI*PetscSinScalar(PETSC_PI*xg[i])*PetscExpReal(-appctx->param.mu*PETSC_PI*PETSC_PI*t)/(2.0+PetscCosScalar(PETSC_PI*xg[i])*PetscExpReal(-appctx->param.mu*PETSC_PI*PETSC_PI*t)); 276 } 277 CHKERRQ(DMDAVecRestoreArray(appctx->da,u,&s)); 278 CHKERRQ(DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg)); 279 return 0; 280 } 281 282 PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx) 283 { 284 AppCtx *appctx = (AppCtx*)ctx; 285 286 PetscFunctionBegin; 287 CHKERRQ(MatMult(appctx->SEMop.grad,globalin,globalout)); /* grad u */ 288 CHKERRQ(VecPointwiseMult(globalout,globalin,globalout)); /* u grad u */ 289 CHKERRQ(VecScale(globalout, -1.0)); 290 CHKERRQ(MatMultAdd(appctx->SEMop.keptstiff,globalin,globalout,globalout)); 291 PetscFunctionReturn(0); 292 } 293 294 /* 295 296 K is the discretiziation of the Laplacian 297 G is the discretization of the gradient 298 299 Computes Jacobian of K u + diag(u) G u which is given by 300 K + diag(u)G + diag(Gu) 301 */ 302 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec globalin,Mat A, Mat B,void *ctx) 303 { 304 AppCtx *appctx = (AppCtx*)ctx; 305 Vec Gglobalin; 306 307 PetscFunctionBegin; 308 /* A = diag(u) G */ 309 310 CHKERRQ(MatCopy(appctx->SEMop.grad,A,SAME_NONZERO_PATTERN)); 311 CHKERRQ(MatDiagonalScale(A,globalin,NULL)); 312 313 /* A = A + diag(Gu) */ 314 CHKERRQ(VecDuplicate(globalin,&Gglobalin)); 315 CHKERRQ(MatMult(appctx->SEMop.grad,globalin,Gglobalin)); 316 CHKERRQ(MatDiagonalSet(A,Gglobalin,ADD_VALUES)); 317 CHKERRQ(VecDestroy(&Gglobalin)); 318 319 /* A = K - A */ 320 CHKERRQ(MatScale(A,-1.0)); 321 CHKERRQ(MatAXPY(A,0.0,appctx->SEMop.keptstiff,SAME_NONZERO_PATTERN)); 322 PetscFunctionReturn(0); 323 } 324 325 /* --------------------------------------------------------------------- */ 326 327 #include "petscblaslapack.h" 328 /* 329 Matrix free operation of 1d Laplacian and Grad for GLL spectral elements 330 */ 331 PetscErrorCode MatMult_Laplacian(Mat A,Vec x,Vec y) 332 { 333 AppCtx *appctx; 334 PetscReal **temp,vv; 335 PetscInt i,j,xs,xn; 336 Vec xlocal,ylocal; 337 const PetscScalar *xl; 338 PetscScalar *yl; 339 PetscBLASInt _One = 1,n; 340 PetscScalar _DOne = 1; 341 342 CHKERRQ(MatShellGetContext(A,&appctx)); 343 CHKERRQ(DMGetLocalVector(appctx->da,&xlocal)); 344 CHKERRQ(DMGlobalToLocalBegin(appctx->da,x,INSERT_VALUES,xlocal)); 345 CHKERRQ(DMGlobalToLocalEnd(appctx->da,x,INSERT_VALUES,xlocal)); 346 CHKERRQ(DMGetLocalVector(appctx->da,&ylocal)); 347 CHKERRQ(VecSet(ylocal,0.0)); 348 CHKERRQ(PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp)); 349 for (i=0; i<appctx->param.N; i++) { 350 vv =-appctx->param.mu*2.0/appctx->param.Le; 351 for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv; 352 } 353 CHKERRQ(DMDAVecGetArrayRead(appctx->da,xlocal,(void*)&xl)); 354 CHKERRQ(DMDAVecGetArray(appctx->da,ylocal,&yl)); 355 CHKERRQ(DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL)); 356 CHKERRQ(PetscBLASIntCast(appctx->param.N,&n)); 357 for (j=xs; j<xs+xn; j += appctx->param.N-1) { 358 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&_DOne,&temp[0][0],&n,&xl[j],&_One,&_DOne,&yl[j],&_One)); 359 } 360 CHKERRQ(DMDAVecRestoreArrayRead(appctx->da,xlocal,(void*)&xl)); 361 CHKERRQ(DMDAVecRestoreArray(appctx->da,ylocal,&yl)); 362 CHKERRQ(PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp)); 363 CHKERRQ(VecSet(y,0.0)); 364 CHKERRQ(DMLocalToGlobalBegin(appctx->da,ylocal,ADD_VALUES,y)); 365 CHKERRQ(DMLocalToGlobalEnd(appctx->da,ylocal,ADD_VALUES,y)); 366 CHKERRQ(DMRestoreLocalVector(appctx->da,&xlocal)); 367 CHKERRQ(DMRestoreLocalVector(appctx->da,&ylocal)); 368 CHKERRQ(VecPointwiseDivide(y,y,appctx->SEMop.mass)); 369 return 0; 370 } 371 372 PetscErrorCode MatMult_Advection(Mat A,Vec x,Vec y) 373 { 374 AppCtx *appctx; 375 PetscReal **temp; 376 PetscInt j,xs,xn; 377 Vec xlocal,ylocal; 378 const PetscScalar *xl; 379 PetscScalar *yl; 380 PetscBLASInt _One = 1,n; 381 PetscScalar _DOne = 1; 382 383 CHKERRQ(MatShellGetContext(A,&appctx)); 384 CHKERRQ(DMGetLocalVector(appctx->da,&xlocal)); 385 CHKERRQ(DMGlobalToLocalBegin(appctx->da,x,INSERT_VALUES,xlocal)); 386 CHKERRQ(DMGlobalToLocalEnd(appctx->da,x,INSERT_VALUES,xlocal)); 387 CHKERRQ(DMGetLocalVector(appctx->da,&ylocal)); 388 CHKERRQ(VecSet(ylocal,0.0)); 389 CHKERRQ(PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp)); 390 CHKERRQ(DMDAVecGetArrayRead(appctx->da,xlocal,(void*)&xl)); 391 CHKERRQ(DMDAVecGetArray(appctx->da,ylocal,&yl)); 392 CHKERRQ(DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL)); 393 CHKERRQ(PetscBLASIntCast(appctx->param.N,&n)); 394 for (j=xs; j<xs+xn; j += appctx->param.N-1) { 395 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&_DOne,&temp[0][0],&n,&xl[j],&_One,&_DOne,&yl[j],&_One)); 396 } 397 CHKERRQ(DMDAVecRestoreArrayRead(appctx->da,xlocal,(void*)&xl)); 398 CHKERRQ(DMDAVecRestoreArray(appctx->da,ylocal,&yl)); 399 CHKERRQ(PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp)); 400 CHKERRQ(VecSet(y,0.0)); 401 CHKERRQ(DMLocalToGlobalBegin(appctx->da,ylocal,ADD_VALUES,y)); 402 CHKERRQ(DMLocalToGlobalEnd(appctx->da,ylocal,ADD_VALUES,y)); 403 CHKERRQ(DMRestoreLocalVector(appctx->da,&xlocal)); 404 CHKERRQ(DMRestoreLocalVector(appctx->da,&ylocal)); 405 CHKERRQ(VecPointwiseDivide(y,y,appctx->SEMop.mass)); 406 CHKERRQ(VecScale(y,-1.0)); 407 return 0; 408 } 409 410 /* 411 RHSMatrixLaplacian - User-provided routine to compute the right-hand-side 412 matrix for the Laplacian operator 413 414 Input Parameters: 415 ts - the TS context 416 t - current time (ignored) 417 X - current solution (ignored) 418 dummy - optional user-defined context, as set by TSetRHSJacobian() 419 420 Output Parameters: 421 AA - Jacobian matrix 422 BB - optionally different matrix from which the preconditioner is built 423 str - flag indicating matrix structure 424 425 */ 426 PetscErrorCode RHSMatrixLaplaciangllDM(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx) 427 { 428 PetscReal **temp; 429 PetscReal vv; 430 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 431 PetscInt i,xs,xn,l,j; 432 PetscInt *rowsDM; 433 PetscBool flg = PETSC_FALSE; 434 435 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-gll_mf",&flg,NULL)); 436 437 if (!flg) { 438 /* 439 Creates the element stiffness matrix for the given gll 440 */ 441 CHKERRQ(PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp)); 442 /* workaround for clang analyzer warning: Division by zero */ 443 PetscCheck(appctx->param.N > 1,PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Spectral element order should be > 1"); 444 445 /* scale by the size of the element */ 446 for (i=0; i<appctx->param.N; i++) { 447 vv=-appctx->param.mu*2.0/appctx->param.Le; 448 for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv; 449 } 450 451 CHKERRQ(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 452 CHKERRQ(DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL)); 453 454 xs = xs/(appctx->param.N-1); 455 xn = xn/(appctx->param.N-1); 456 457 CHKERRQ(PetscMalloc1(appctx->param.N,&rowsDM)); 458 /* 459 loop over local elements 460 */ 461 for (j=xs; j<xs+xn; j++) { 462 for (l=0; l<appctx->param.N; l++) { 463 rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l; 464 } 465 CHKERRQ(MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES)); 466 } 467 CHKERRQ(PetscFree(rowsDM)); 468 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 469 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 470 CHKERRQ(VecReciprocal(appctx->SEMop.mass)); 471 CHKERRQ(MatDiagonalScale(A,appctx->SEMop.mass,0)); 472 CHKERRQ(VecReciprocal(appctx->SEMop.mass)); 473 474 CHKERRQ(PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp)); 475 } else { 476 CHKERRQ(MatSetType(A,MATSHELL)); 477 CHKERRQ(MatSetUp(A)); 478 CHKERRQ(MatShellSetContext(A,appctx)); 479 CHKERRQ(MatShellSetOperation(A,MATOP_MULT,(void (*)(void))MatMult_Laplacian)); 480 } 481 return 0; 482 } 483 484 /* 485 RHSMatrixAdvection - User-provided routine to compute the right-hand-side 486 matrix for the Advection (gradient) operator. 487 488 Input Parameters: 489 ts - the TS context 490 t - current time 491 global_in - global input vector 492 dummy - optional user-defined context, as set by TSetRHSJacobian() 493 494 Output Parameters: 495 AA - Jacobian matrix 496 BB - optionally different preconditioning matrix 497 str - flag indicating matrix structure 498 499 */ 500 PetscErrorCode RHSMatrixAdvectiongllDM(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx) 501 { 502 PetscReal **temp; 503 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 504 PetscInt xs,xn,l,j; 505 PetscInt *rowsDM; 506 PetscBool flg = PETSC_FALSE; 507 508 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-gll_mf",&flg,NULL)); 509 510 if (!flg) { 511 /* 512 Creates the advection matrix for the given gll 513 */ 514 CHKERRQ(PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp)); 515 CHKERRQ(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 516 CHKERRQ(DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL)); 517 xs = xs/(appctx->param.N-1); 518 xn = xn/(appctx->param.N-1); 519 520 CHKERRQ(PetscMalloc1(appctx->param.N,&rowsDM)); 521 for (j=xs; j<xs+xn; j++) { 522 for (l=0; l<appctx->param.N; l++) { 523 rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l; 524 } 525 CHKERRQ(MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES)); 526 } 527 CHKERRQ(PetscFree(rowsDM)); 528 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 529 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 530 531 CHKERRQ(VecReciprocal(appctx->SEMop.mass)); 532 CHKERRQ(MatDiagonalScale(A,appctx->SEMop.mass,0)); 533 CHKERRQ(VecReciprocal(appctx->SEMop.mass)); 534 535 CHKERRQ(PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp)); 536 } else { 537 CHKERRQ(MatSetType(A,MATSHELL)); 538 CHKERRQ(MatSetUp(A)); 539 CHKERRQ(MatShellSetContext(A,appctx)); 540 CHKERRQ(MatShellSetOperation(A,MATOP_MULT,(void (*)(void))MatMult_Advection)); 541 } 542 return 0; 543 } 544 545 /*TEST 546 547 build: 548 requires: !complex 549 550 test: 551 suffix: 1 552 requires: !single 553 554 test: 555 suffix: 2 556 nsize: 5 557 requires: !single 558 559 test: 560 suffix: 3 561 requires: !single 562 args: -ts_view -ts_type beuler -gll_mf -pc_type none -ts_max_steps 5 -ts_monitor_error 563 564 test: 565 suffix: 4 566 requires: !single 567 args: -ts_view -ts_type beuler -pc_type none -ts_max_steps 5 -ts_monitor_error 568 569 TEST*/ 570