1 2 static char help[] ="Solves a simple data assimilation problem with one dimensional Burger's equation using TSAdjoint\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 Concepts: adjoints 13 Processors: n 14 */ 15 16 /* ------------------------------------------------------------------------ 17 18 This program uses the one-dimensional Burger's equation 19 u_t = mu*u_xx - u u_x, 20 on the domain 0 <= x <= 1, with periodic boundary conditions 21 22 to demonstrate solving a data assimilation problem of finding the initial conditions 23 to produce a given solution at a fixed time. 24 25 The operators are discretized with the spectral element method 26 27 See the paper PDE-CONSTRAINED OPTIMIZATION WITH SPECTRAL ELEMENTS USING PETSC AND TAO 28 by OANA MARIN, EMIL CONSTANTINESCU, AND BARRY SMITH for details on the exact solution 29 used 30 31 ------------------------------------------------------------------------- */ 32 33 #include <petsctao.h> 34 #include <petscts.h> 35 #include <petscdt.h> 36 #include <petscdraw.h> 37 #include <petscdmda.h> 38 39 /* 40 User-defined application context - contains data needed by the 41 application-provided call-back routines. 42 */ 43 44 typedef struct { 45 PetscInt n; /* number of nodes */ 46 PetscReal *nodes; /* GLL nodes */ 47 PetscReal *weights; /* GLL weights */ 48 } PetscGLL; 49 50 typedef struct { 51 PetscInt N; /* grid points per elements*/ 52 PetscInt E; /* number of elements */ 53 PetscReal tol_L2,tol_max; /* error norms */ 54 PetscInt steps; /* number of timesteps */ 55 PetscReal Tend; /* endtime */ 56 PetscReal mu; /* viscosity */ 57 PetscReal L; /* total length of domain */ 58 PetscReal Le; 59 PetscReal Tadj; 60 } PetscParam; 61 62 typedef struct { 63 Vec obj; /* desired end state */ 64 Vec grid; /* total grid */ 65 Vec grad; 66 Vec ic; 67 Vec curr_sol; 68 Vec true_solution; /* actual initial conditions for the final solution */ 69 } PetscData; 70 71 typedef struct { 72 Vec grid; /* total grid */ 73 Vec mass; /* mass matrix for total integration */ 74 Mat stiff; /* stifness matrix */ 75 Mat keptstiff; 76 Mat grad; 77 PetscGLL gll; 78 } PetscSEMOperators; 79 80 typedef struct { 81 DM da; /* distributed array data structure */ 82 PetscSEMOperators SEMop; 83 PetscParam param; 84 PetscData dat; 85 TS ts; 86 PetscReal initial_dt; 87 } AppCtx; 88 89 /* 90 User-defined routines 91 */ 92 extern PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*); 93 extern PetscErrorCode RHSMatrixLaplaciangllDM(TS,PetscReal,Vec,Mat,Mat,void*); 94 extern PetscErrorCode RHSMatrixAdvectiongllDM(TS,PetscReal,Vec,Mat,Mat,void*); 95 extern PetscErrorCode InitialConditions(Vec,AppCtx*); 96 extern PetscErrorCode TrueSolution(Vec,AppCtx*); 97 extern PetscErrorCode ComputeObjective(PetscReal,Vec,AppCtx*); 98 extern PetscErrorCode MonitorError(Tao,void*); 99 extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*); 100 extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*); 101 102 int main(int argc,char **argv) 103 { 104 AppCtx appctx; /* user-defined application context */ 105 Tao tao; 106 Vec u; /* approximate solution vector */ 107 PetscInt i, xs, xm, ind, j, lenglob; 108 PetscReal x, *wrk_ptr1, *wrk_ptr2; 109 MatNullSpace nsp; 110 PetscMPIInt size; 111 112 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 113 Initialize program and set problem parameters 114 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 115 PetscFunctionBegin; 116 117 CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 118 119 /*initialize parameters */ 120 appctx.param.N = 10; /* order of the spectral element */ 121 appctx.param.E = 10; /* number of elements */ 122 appctx.param.L = 4.0; /* length of the domain */ 123 appctx.param.mu = 0.01; /* diffusion coefficient */ 124 appctx.initial_dt = 5e-3; 125 appctx.param.steps = PETSC_MAX_INT; 126 appctx.param.Tend = 4; 127 128 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-N",&appctx.param.N,NULL)); 129 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-E",&appctx.param.E,NULL)); 130 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-Tend",&appctx.param.Tend,NULL)); 131 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-mu",&appctx.param.mu,NULL)); 132 appctx.param.Le = appctx.param.L/appctx.param.E; 133 134 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 135 PetscCheck((appctx.param.E % size) == 0,PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Number of elements must be divisible by number of processes"); 136 137 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 138 Create GLL data structures 139 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 140 CHKERRQ(PetscMalloc2(appctx.param.N,&appctx.SEMop.gll.nodes,appctx.param.N,&appctx.SEMop.gll.weights)); 141 CHKERRQ(PetscDTGaussLobattoLegendreQuadrature(appctx.param.N,PETSCGAUSSLOBATTOLEGENDRE_VIA_LINEAR_ALGEBRA,appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights)); 142 appctx.SEMop.gll.n = appctx.param.N; 143 lenglob = appctx.param.E*(appctx.param.N-1); 144 145 /* 146 Create distributed array (DMDA) to manage parallel grid and vectors 147 and to set up the ghost point communication pattern. There are E*(Nl-1)+1 148 total grid values spread equally among all the processors, except first and last 149 */ 150 151 CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,lenglob,1,1,NULL,&appctx.da)); 152 CHKERRQ(DMSetFromOptions(appctx.da)); 153 CHKERRQ(DMSetUp(appctx.da)); 154 155 /* 156 Extract global and local vectors from DMDA; we use these to store the 157 approximate solution. Then duplicate these for remaining vectors that 158 have the same types. 159 */ 160 161 CHKERRQ(DMCreateGlobalVector(appctx.da,&u)); 162 CHKERRQ(VecDuplicate(u,&appctx.dat.ic)); 163 CHKERRQ(VecDuplicate(u,&appctx.dat.true_solution)); 164 CHKERRQ(VecDuplicate(u,&appctx.dat.obj)); 165 CHKERRQ(VecDuplicate(u,&appctx.SEMop.grid)); 166 CHKERRQ(VecDuplicate(u,&appctx.SEMop.mass)); 167 CHKERRQ(VecDuplicate(u,&appctx.dat.curr_sol)); 168 169 CHKERRQ(DMDAGetCorners(appctx.da,&xs,NULL,NULL,&xm,NULL,NULL)); 170 CHKERRQ(DMDAVecGetArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1)); 171 CHKERRQ(DMDAVecGetArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2)); 172 173 /* Compute function over the locally owned part of the grid */ 174 175 xs=xs/(appctx.param.N-1); 176 xm=xm/(appctx.param.N-1); 177 178 /* 179 Build total grid and mass over entire mesh (multi-elemental) 180 */ 181 182 for (i=xs; i<xs+xm; i++) { 183 for (j=0; j<appctx.param.N-1; j++) { 184 x = (appctx.param.Le/2.0)*(appctx.SEMop.gll.nodes[j]+1.0)+appctx.param.Le*i; 185 ind=i*(appctx.param.N-1)+j; 186 wrk_ptr1[ind]=x; 187 wrk_ptr2[ind]=.5*appctx.param.Le*appctx.SEMop.gll.weights[j]; 188 if (j==0) wrk_ptr2[ind]+=.5*appctx.param.Le*appctx.SEMop.gll.weights[j]; 189 } 190 } 191 CHKERRQ(DMDAVecRestoreArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1)); 192 CHKERRQ(DMDAVecRestoreArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2)); 193 194 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 195 Create matrix data structure; set matrix evaluation routine. 196 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 197 CHKERRQ(DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE)); 198 CHKERRQ(DMCreateMatrix(appctx.da,&appctx.SEMop.stiff)); 199 CHKERRQ(DMCreateMatrix(appctx.da,&appctx.SEMop.grad)); 200 /* 201 For linear problems with a time-dependent f(u,t) in the equation 202 u_t = f(u,t), the user provides the discretized right-hand-side 203 as a time-dependent matrix. 204 */ 205 CHKERRQ(RHSMatrixLaplaciangllDM(appctx.ts,0.0,u,appctx.SEMop.stiff,appctx.SEMop.stiff,&appctx)); 206 CHKERRQ(RHSMatrixAdvectiongllDM(appctx.ts,0.0,u,appctx.SEMop.grad,appctx.SEMop.grad,&appctx)); 207 /* 208 For linear problems with a time-dependent f(u,t) in the equation 209 u_t = f(u,t), the user provides the discretized right-hand-side 210 as a time-dependent matrix. 211 */ 212 213 CHKERRQ(MatDuplicate(appctx.SEMop.stiff,MAT_COPY_VALUES,&appctx.SEMop.keptstiff)); 214 215 /* attach the null space to the matrix, this probably is not needed but does no harm */ 216 CHKERRQ(MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp)); 217 CHKERRQ(MatSetNullSpace(appctx.SEMop.stiff,nsp)); 218 CHKERRQ(MatSetNullSpace(appctx.SEMop.keptstiff,nsp)); 219 CHKERRQ(MatNullSpaceTest(nsp,appctx.SEMop.stiff,NULL)); 220 CHKERRQ(MatNullSpaceDestroy(&nsp)); 221 /* attach the null space to the matrix, this probably is not needed but does no harm */ 222 CHKERRQ(MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp)); 223 CHKERRQ(MatSetNullSpace(appctx.SEMop.grad,nsp)); 224 CHKERRQ(MatNullSpaceTest(nsp,appctx.SEMop.grad,NULL)); 225 CHKERRQ(MatNullSpaceDestroy(&nsp)); 226 227 /* Create the TS solver that solves the ODE and its adjoint; set its options */ 228 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&appctx.ts)); 229 CHKERRQ(TSSetProblemType(appctx.ts,TS_NONLINEAR)); 230 CHKERRQ(TSSetType(appctx.ts,TSRK)); 231 CHKERRQ(TSSetDM(appctx.ts,appctx.da)); 232 CHKERRQ(TSSetTime(appctx.ts,0.0)); 233 CHKERRQ(TSSetTimeStep(appctx.ts,appctx.initial_dt)); 234 CHKERRQ(TSSetMaxSteps(appctx.ts,appctx.param.steps)); 235 CHKERRQ(TSSetMaxTime(appctx.ts,appctx.param.Tend)); 236 CHKERRQ(TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP)); 237 CHKERRQ(TSSetTolerances(appctx.ts,1e-7,NULL,1e-7,NULL)); 238 CHKERRQ(TSSetFromOptions(appctx.ts)); 239 /* Need to save initial timestep user may have set with -ts_dt so it can be reset for each new TSSolve() */ 240 CHKERRQ(TSGetTimeStep(appctx.ts,&appctx.initial_dt)); 241 CHKERRQ(TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx)); 242 CHKERRQ(TSSetRHSJacobian(appctx.ts,appctx.SEMop.stiff,appctx.SEMop.stiff,RHSJacobian,&appctx)); 243 244 /* Set Objective and Initial conditions for the problem and compute Objective function (evolution of true_solution to final time */ 245 CHKERRQ(InitialConditions(appctx.dat.ic,&appctx)); 246 CHKERRQ(TrueSolution(appctx.dat.true_solution,&appctx)); 247 CHKERRQ(ComputeObjective(appctx.param.Tend,appctx.dat.obj,&appctx)); 248 249 CHKERRQ(TSSetSaveTrajectory(appctx.ts)); 250 CHKERRQ(TSSetFromOptions(appctx.ts)); 251 252 /* Create TAO solver and set desired solution method */ 253 CHKERRQ(TaoCreate(PETSC_COMM_WORLD,&tao)); 254 CHKERRQ(TaoSetMonitor(tao,MonitorError,&appctx,NULL)); 255 CHKERRQ(TaoSetType(tao,TAOBQNLS)); 256 CHKERRQ(TaoSetSolution(tao,appctx.dat.ic)); 257 /* Set routine for function and gradient evaluation */ 258 CHKERRQ(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&appctx)); 259 /* Check for any TAO command line options */ 260 CHKERRQ(TaoSetTolerances(tao,1e-8,PETSC_DEFAULT,PETSC_DEFAULT)); 261 CHKERRQ(TaoSetFromOptions(tao)); 262 CHKERRQ(TaoSolve(tao)); 263 264 CHKERRQ(TaoDestroy(&tao)); 265 CHKERRQ(MatDestroy(&appctx.SEMop.stiff)); 266 CHKERRQ(MatDestroy(&appctx.SEMop.keptstiff)); 267 CHKERRQ(MatDestroy(&appctx.SEMop.grad)); 268 CHKERRQ(VecDestroy(&u)); 269 CHKERRQ(VecDestroy(&appctx.dat.ic)); 270 CHKERRQ(VecDestroy(&appctx.dat.true_solution)); 271 CHKERRQ(VecDestroy(&appctx.dat.obj)); 272 CHKERRQ(VecDestroy(&appctx.SEMop.grid)); 273 CHKERRQ(VecDestroy(&appctx.SEMop.mass)); 274 CHKERRQ(VecDestroy(&appctx.dat.curr_sol)); 275 CHKERRQ(PetscFree2(appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights)); 276 CHKERRQ(DMDestroy(&appctx.da)); 277 CHKERRQ(TSDestroy(&appctx.ts)); 278 279 /* 280 Always call PetscFinalize() before exiting a program. This routine 281 - finalizes the PETSc libraries as well as MPI 282 - provides summary and diagnostic information if certain runtime 283 options are chosen (e.g., -log_summary). 284 */ 285 CHKERRQ(PetscFinalize()); 286 return 0; 287 } 288 289 /* --------------------------------------------------------------------- */ 290 /* 291 InitialConditions - Computes the initial conditions for the Tao optimization solve (these are also initial conditions for the first TSSolve() 292 293 The routine TrueSolution() computes the true solution for the Tao optimization solve which means they are the initial conditions for the objective function 294 295 Input Parameter: 296 u - uninitialized solution vector (global) 297 appctx - user-defined application context 298 299 Output Parameter: 300 u - vector with solution at initial time (global) 301 */ 302 PetscErrorCode InitialConditions(Vec u,AppCtx *appctx) 303 { 304 PetscScalar *s; 305 const PetscScalar *xg; 306 PetscInt i,xs,xn; 307 308 PetscFunctionBegin; 309 CHKERRQ(DMDAVecGetArray(appctx->da,u,&s)); 310 CHKERRQ(DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg)); 311 CHKERRQ(DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL)); 312 for (i=xs; i<xs+xn; i++) { 313 s[i]=2.0*appctx->param.mu*PETSC_PI*PetscSinScalar(PETSC_PI*xg[i])/(2.0+PetscCosScalar(PETSC_PI*xg[i]))+0.25*PetscExpReal(-4.0*PetscPowReal(xg[i]-2.0,2.0)); 314 } 315 CHKERRQ(DMDAVecRestoreArray(appctx->da,u,&s)); 316 CHKERRQ(DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg)); 317 PetscFunctionReturn(0); 318 } 319 320 /* 321 TrueSolution() computes the true solution for the Tao optimization solve which means they are the initial conditions for the objective function. 322 323 InitialConditions() computes the initial conditions for the beginning of the Tao iterations 324 325 Input Parameter: 326 u - uninitialized solution vector (global) 327 appctx - user-defined application context 328 329 Output Parameter: 330 u - vector with solution at initial time (global) 331 */ 332 PetscErrorCode TrueSolution(Vec u,AppCtx *appctx) 333 { 334 PetscScalar *s; 335 const PetscScalar *xg; 336 PetscInt i,xs,xn; 337 338 PetscFunctionBegin; 339 CHKERRQ(DMDAVecGetArray(appctx->da,u,&s)); 340 CHKERRQ(DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg)); 341 CHKERRQ(DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL)); 342 for (i=xs; i<xs+xn; i++) { 343 s[i]=2.0*appctx->param.mu*PETSC_PI*PetscSinScalar(PETSC_PI*xg[i])/(2.0+PetscCosScalar(PETSC_PI*xg[i])); 344 } 345 CHKERRQ(DMDAVecRestoreArray(appctx->da,u,&s)); 346 CHKERRQ(DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg)); 347 PetscFunctionReturn(0); 348 } 349 /* --------------------------------------------------------------------- */ 350 /* 351 Sets the desired profile for the final end time 352 353 Input Parameters: 354 t - final time 355 obj - vector storing the desired profile 356 appctx - user-defined application context 357 358 */ 359 PetscErrorCode ComputeObjective(PetscReal t,Vec obj,AppCtx *appctx) 360 { 361 PetscScalar *s; 362 const PetscScalar *xg; 363 PetscInt i, xs,xn; 364 365 PetscFunctionBegin; 366 CHKERRQ(DMDAVecGetArray(appctx->da,obj,&s)); 367 CHKERRQ(DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg)); 368 CHKERRQ(DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL)); 369 for (i=xs; i<xs+xn; i++) { 370 s[i]=2.0*appctx->param.mu*PETSC_PI*PetscSinScalar(PETSC_PI*xg[i])*PetscExpScalar(-PETSC_PI*PETSC_PI*t*appctx->param.mu)\ 371 /(2.0+PetscExpScalar(-PETSC_PI*PETSC_PI*t*appctx->param.mu)*PetscCosScalar(PETSC_PI*xg[i])); 372 } 373 CHKERRQ(DMDAVecRestoreArray(appctx->da,obj,&s)); 374 CHKERRQ(DMDAVecRestoreArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg)); 375 PetscFunctionReturn(0); 376 } 377 378 PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx) 379 { 380 AppCtx *appctx = (AppCtx*)ctx; 381 382 PetscFunctionBegin; 383 CHKERRQ(MatMult(appctx->SEMop.grad,globalin,globalout)); /* grad u */ 384 CHKERRQ(VecPointwiseMult(globalout,globalin,globalout)); /* u grad u */ 385 CHKERRQ(VecScale(globalout, -1.0)); 386 CHKERRQ(MatMultAdd(appctx->SEMop.keptstiff,globalin,globalout,globalout)); 387 PetscFunctionReturn(0); 388 } 389 390 /* 391 392 K is the discretiziation of the Laplacian 393 G is the discretization of the gradient 394 395 Computes Jacobian of K u + diag(u) G u which is given by 396 K + diag(u)G + diag(Gu) 397 */ 398 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec globalin,Mat A, Mat B,void *ctx) 399 { 400 AppCtx *appctx = (AppCtx*)ctx; 401 Vec Gglobalin; 402 403 PetscFunctionBegin; 404 /* A = diag(u) G */ 405 406 CHKERRQ(MatCopy(appctx->SEMop.grad,A,SAME_NONZERO_PATTERN)); 407 CHKERRQ(MatDiagonalScale(A,globalin,NULL)); 408 409 /* A = A + diag(Gu) */ 410 CHKERRQ(VecDuplicate(globalin,&Gglobalin)); 411 CHKERRQ(MatMult(appctx->SEMop.grad,globalin,Gglobalin)); 412 CHKERRQ(MatDiagonalSet(A,Gglobalin,ADD_VALUES)); 413 CHKERRQ(VecDestroy(&Gglobalin)); 414 415 /* A = K - A */ 416 CHKERRQ(MatScale(A,-1.0)); 417 CHKERRQ(MatAXPY(A,1.0,appctx->SEMop.keptstiff,SAME_NONZERO_PATTERN)); 418 PetscFunctionReturn(0); 419 } 420 421 /* --------------------------------------------------------------------- */ 422 423 /* 424 RHSMatrixLaplacian - User-provided routine to compute the right-hand-side 425 matrix for the heat equation. 426 427 Input Parameters: 428 ts - the TS context 429 t - current time (ignored) 430 X - current solution (ignored) 431 dummy - optional user-defined context, as set by TSetRHSJacobian() 432 433 Output Parameters: 434 AA - Jacobian matrix 435 BB - optionally different matrix from which the preconditioner is built 436 str - flag indicating matrix structure 437 438 */ 439 PetscErrorCode RHSMatrixLaplaciangllDM(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx) 440 { 441 PetscReal **temp; 442 PetscReal vv; 443 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 444 PetscInt i,xs,xn,l,j; 445 PetscInt *rowsDM; 446 447 PetscFunctionBegin; 448 /* 449 Creates the element stiffness matrix for the given gll 450 */ 451 CHKERRQ(PetscGaussLobattoLegendreElementLaplacianCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp)); 452 /* workaround for clang analyzer warning: Division by zero */ 453 PetscCheck(appctx->param.N > 1,PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Spectral element order should be > 1"); 454 455 /* scale by the size of the element */ 456 for (i=0; i<appctx->param.N; i++) { 457 vv=-appctx->param.mu*2.0/appctx->param.Le; 458 for (j=0; j<appctx->param.N; j++) temp[i][j]=temp[i][j]*vv; 459 } 460 461 CHKERRQ(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 462 CHKERRQ(DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL)); 463 464 xs = xs/(appctx->param.N-1); 465 xn = xn/(appctx->param.N-1); 466 467 CHKERRQ(PetscMalloc1(appctx->param.N,&rowsDM)); 468 /* 469 loop over local elements 470 */ 471 for (j=xs; j<xs+xn; j++) { 472 for (l=0; l<appctx->param.N; l++) { 473 rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l; 474 } 475 CHKERRQ(MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES)); 476 } 477 CHKERRQ(PetscFree(rowsDM)); 478 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 479 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 480 CHKERRQ(VecReciprocal(appctx->SEMop.mass)); 481 CHKERRQ(MatDiagonalScale(A,appctx->SEMop.mass,0)); 482 CHKERRQ(VecReciprocal(appctx->SEMop.mass)); 483 484 CHKERRQ(PetscGaussLobattoLegendreElementLaplacianDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp)); 485 PetscFunctionReturn(0); 486 } 487 488 /* 489 RHSMatrixAdvection - User-provided routine to compute the right-hand-side 490 matrix for the Advection equation. 491 492 Input Parameters: 493 ts - the TS context 494 t - current time 495 global_in - global input vector 496 dummy - optional user-defined context, as set by TSetRHSJacobian() 497 498 Output Parameters: 499 AA - Jacobian matrix 500 BB - optionally different preconditioning matrix 501 str - flag indicating matrix structure 502 503 */ 504 PetscErrorCode RHSMatrixAdvectiongllDM(TS ts,PetscReal t,Vec X,Mat A,Mat BB,void *ctx) 505 { 506 PetscReal **temp; 507 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 508 PetscInt xs,xn,l,j; 509 PetscInt *rowsDM; 510 511 PetscFunctionBegin; 512 /* 513 Creates the advection matrix for the given gll 514 */ 515 CHKERRQ(PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp)); 516 CHKERRQ(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 517 518 CHKERRQ(DMDAGetCorners(appctx->da,&xs,NULL,NULL,&xn,NULL,NULL)); 519 520 xs = xs/(appctx->param.N-1); 521 xn = xn/(appctx->param.N-1); 522 523 CHKERRQ(PetscMalloc1(appctx->param.N,&rowsDM)); 524 for (j=xs; j<xs+xn; j++) { 525 for (l=0; l<appctx->param.N; l++) { 526 rowsDM[l] = 1+(j-xs)*(appctx->param.N-1)+l; 527 } 528 CHKERRQ(MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES)); 529 } 530 CHKERRQ(PetscFree(rowsDM)); 531 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 532 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 533 534 CHKERRQ(VecReciprocal(appctx->SEMop.mass)); 535 CHKERRQ(MatDiagonalScale(A,appctx->SEMop.mass,0)); 536 CHKERRQ(VecReciprocal(appctx->SEMop.mass)); 537 CHKERRQ(PetscGaussLobattoLegendreElementAdvectionDestroy(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp)); 538 PetscFunctionReturn(0); 539 } 540 /* ------------------------------------------------------------------ */ 541 /* 542 FormFunctionGradient - Evaluates the function and corresponding gradient. 543 544 Input Parameters: 545 tao - the Tao context 546 IC - the input vector 547 ctx - optional user-defined context, as set when calling TaoSetObjectiveAndGradient() 548 549 Output Parameters: 550 f - the newly evaluated function 551 G - the newly evaluated gradient 552 553 Notes: 554 555 The forward equation is 556 M u_t = F(U) 557 which is converted to 558 u_t = M^{-1} F(u) 559 in the user code since TS has no direct way of providing a mass matrix. The Jacobian of this is 560 M^{-1} J 561 where J is the Jacobian of F. Now the adjoint equation is 562 M v_t = J^T v 563 but TSAdjoint does not solve this since it can only solve the transposed system for the 564 Jacobian the user provided. Hence TSAdjoint solves 565 w_t = J^T M^{-1} w (where w = M v) 566 since there is no way to indicate the mass matrix as a separate entity to TS. Thus one 567 must be careful in initializing the "adjoint equation" and using the result. This is 568 why 569 G = -2 M(u(T) - u_d) 570 below (instead of -2(u(T) - u_d) and why the result is 571 G = G/appctx->SEMop.mass (that is G = M^{-1}w) 572 below (instead of just the result of the "adjoint solve"). 573 574 */ 575 PetscErrorCode FormFunctionGradient(Tao tao,Vec IC,PetscReal *f,Vec G,void *ctx) 576 { 577 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 578 Vec temp; 579 PetscInt its; 580 PetscReal ff, gnorm, cnorm, xdiff,errex; 581 TaoConvergedReason reason; 582 583 PetscFunctionBegin; 584 CHKERRQ(TSSetTime(appctx->ts,0.0)); 585 CHKERRQ(TSSetStepNumber(appctx->ts,0)); 586 CHKERRQ(TSSetTimeStep(appctx->ts,appctx->initial_dt)); 587 CHKERRQ(VecCopy(IC,appctx->dat.curr_sol)); 588 589 CHKERRQ(TSSolve(appctx->ts,appctx->dat.curr_sol)); 590 591 CHKERRQ(VecWAXPY(G,-1.0,appctx->dat.curr_sol,appctx->dat.obj)); 592 593 /* 594 Compute the L2-norm of the objective function, cost function is f 595 */ 596 CHKERRQ(VecDuplicate(G,&temp)); 597 CHKERRQ(VecPointwiseMult(temp,G,G)); 598 CHKERRQ(VecDot(temp,appctx->SEMop.mass,f)); 599 600 /* local error evaluation */ 601 CHKERRQ(VecWAXPY(temp,-1.0,appctx->dat.ic,appctx->dat.true_solution)); 602 CHKERRQ(VecPointwiseMult(temp,temp,temp)); 603 /* for error evaluation */ 604 CHKERRQ(VecDot(temp,appctx->SEMop.mass,&errex)); 605 CHKERRQ(VecDestroy(&temp)); 606 errex = PetscSqrtReal(errex); 607 608 /* 609 Compute initial conditions for the adjoint integration. See Notes above 610 */ 611 612 CHKERRQ(VecScale(G, -2.0)); 613 CHKERRQ(VecPointwiseMult(G,G,appctx->SEMop.mass)); 614 CHKERRQ(TSSetCostGradients(appctx->ts,1,&G,NULL)); 615 CHKERRQ(TSAdjointSolve(appctx->ts)); 616 CHKERRQ(VecPointwiseDivide(G,G,appctx->SEMop.mass)); 617 618 CHKERRQ(TaoGetSolutionStatus(tao, &its, &ff, &gnorm, &cnorm, &xdiff, &reason)); 619 PetscFunctionReturn(0); 620 } 621 622 PetscErrorCode MonitorError(Tao tao,void *ctx) 623 { 624 AppCtx *appctx = (AppCtx*)ctx; 625 Vec temp; 626 PetscReal nrm; 627 628 PetscFunctionBegin; 629 CHKERRQ(VecDuplicate(appctx->dat.ic,&temp)); 630 CHKERRQ(VecWAXPY(temp,-1.0,appctx->dat.ic,appctx->dat.true_solution)); 631 CHKERRQ(VecPointwiseMult(temp,temp,temp)); 632 CHKERRQ(VecDot(temp,appctx->SEMop.mass,&nrm)); 633 CHKERRQ(VecDestroy(&temp)); 634 nrm = PetscSqrtReal(nrm); 635 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Error for initial conditions %g\n",(double)nrm)); 636 PetscFunctionReturn(0); 637 } 638 639 /*TEST 640 641 build: 642 requires: !complex 643 644 test: 645 args: -tao_max_it 5 -tao_gatol 1.e-4 646 requires: !single 647 648 test: 649 suffix: 2 650 nsize: 2 651 args: -tao_max_it 5 -tao_gatol 1.e-4 652 requires: !single 653 654 TEST*/ 655