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 PetscCall(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 PetscCall(PetscOptionsGetInt(NULL,NULL,"-N",&appctx.param.N,NULL)); 129 PetscCall(PetscOptionsGetInt(NULL,NULL,"-E",&appctx.param.E,NULL)); 130 PetscCall(PetscOptionsGetReal(NULL,NULL,"-Tend",&appctx.param.Tend,NULL)); 131 PetscCall(PetscOptionsGetReal(NULL,NULL,"-mu",&appctx.param.mu,NULL)); 132 appctx.param.Le = appctx.param.L/appctx.param.E; 133 134 PetscCallMPI(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 PetscCall(PetscMalloc2(appctx.param.N,&appctx.SEMop.gll.nodes,appctx.param.N,&appctx.SEMop.gll.weights)); 141 PetscCall(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 PetscCall(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,lenglob,1,1,NULL,&appctx.da)); 152 PetscCall(DMSetFromOptions(appctx.da)); 153 PetscCall(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 PetscCall(DMCreateGlobalVector(appctx.da,&u)); 162 PetscCall(VecDuplicate(u,&appctx.dat.ic)); 163 PetscCall(VecDuplicate(u,&appctx.dat.true_solution)); 164 PetscCall(VecDuplicate(u,&appctx.dat.obj)); 165 PetscCall(VecDuplicate(u,&appctx.SEMop.grid)); 166 PetscCall(VecDuplicate(u,&appctx.SEMop.mass)); 167 PetscCall(VecDuplicate(u,&appctx.dat.curr_sol)); 168 169 PetscCall(DMDAGetCorners(appctx.da,&xs,NULL,NULL,&xm,NULL,NULL)); 170 PetscCall(DMDAVecGetArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1)); 171 PetscCall(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 PetscCall(DMDAVecRestoreArray(appctx.da,appctx.SEMop.grid,&wrk_ptr1)); 192 PetscCall(DMDAVecRestoreArray(appctx.da,appctx.SEMop.mass,&wrk_ptr2)); 193 194 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 195 Create matrix data structure; set matrix evaluation routine. 196 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 197 PetscCall(DMSetMatrixPreallocateOnly(appctx.da, PETSC_TRUE)); 198 PetscCall(DMCreateMatrix(appctx.da,&appctx.SEMop.stiff)); 199 PetscCall(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 PetscCall(RHSMatrixLaplaciangllDM(appctx.ts,0.0,u,appctx.SEMop.stiff,appctx.SEMop.stiff,&appctx)); 206 PetscCall(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 PetscCall(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 PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp)); 217 PetscCall(MatSetNullSpace(appctx.SEMop.stiff,nsp)); 218 PetscCall(MatSetNullSpace(appctx.SEMop.keptstiff,nsp)); 219 PetscCall(MatNullSpaceTest(nsp,appctx.SEMop.stiff,NULL)); 220 PetscCall(MatNullSpaceDestroy(&nsp)); 221 /* attach the null space to the matrix, this probably is not needed but does no harm */ 222 PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nsp)); 223 PetscCall(MatSetNullSpace(appctx.SEMop.grad,nsp)); 224 PetscCall(MatNullSpaceTest(nsp,appctx.SEMop.grad,NULL)); 225 PetscCall(MatNullSpaceDestroy(&nsp)); 226 227 /* Create the TS solver that solves the ODE and its adjoint; set its options */ 228 PetscCall(TSCreate(PETSC_COMM_WORLD,&appctx.ts)); 229 PetscCall(TSSetProblemType(appctx.ts,TS_NONLINEAR)); 230 PetscCall(TSSetType(appctx.ts,TSRK)); 231 PetscCall(TSSetDM(appctx.ts,appctx.da)); 232 PetscCall(TSSetTime(appctx.ts,0.0)); 233 PetscCall(TSSetTimeStep(appctx.ts,appctx.initial_dt)); 234 PetscCall(TSSetMaxSteps(appctx.ts,appctx.param.steps)); 235 PetscCall(TSSetMaxTime(appctx.ts,appctx.param.Tend)); 236 PetscCall(TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP)); 237 PetscCall(TSSetTolerances(appctx.ts,1e-7,NULL,1e-7,NULL)); 238 PetscCall(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 PetscCall(TSGetTimeStep(appctx.ts,&appctx.initial_dt)); 241 PetscCall(TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx)); 242 PetscCall(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 PetscCall(InitialConditions(appctx.dat.ic,&appctx)); 246 PetscCall(TrueSolution(appctx.dat.true_solution,&appctx)); 247 PetscCall(ComputeObjective(appctx.param.Tend,appctx.dat.obj,&appctx)); 248 249 PetscCall(TSSetSaveTrajectory(appctx.ts)); 250 PetscCall(TSSetFromOptions(appctx.ts)); 251 252 /* Create TAO solver and set desired solution method */ 253 PetscCall(TaoCreate(PETSC_COMM_WORLD,&tao)); 254 PetscCall(TaoSetMonitor(tao,MonitorError,&appctx,NULL)); 255 PetscCall(TaoSetType(tao,TAOBQNLS)); 256 PetscCall(TaoSetSolution(tao,appctx.dat.ic)); 257 /* Set routine for function and gradient evaluation */ 258 PetscCall(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&appctx)); 259 /* Check for any TAO command line options */ 260 PetscCall(TaoSetTolerances(tao,1e-8,PETSC_DEFAULT,PETSC_DEFAULT)); 261 PetscCall(TaoSetFromOptions(tao)); 262 PetscCall(TaoSolve(tao)); 263 264 PetscCall(TaoDestroy(&tao)); 265 PetscCall(MatDestroy(&appctx.SEMop.stiff)); 266 PetscCall(MatDestroy(&appctx.SEMop.keptstiff)); 267 PetscCall(MatDestroy(&appctx.SEMop.grad)); 268 PetscCall(VecDestroy(&u)); 269 PetscCall(VecDestroy(&appctx.dat.ic)); 270 PetscCall(VecDestroy(&appctx.dat.true_solution)); 271 PetscCall(VecDestroy(&appctx.dat.obj)); 272 PetscCall(VecDestroy(&appctx.SEMop.grid)); 273 PetscCall(VecDestroy(&appctx.SEMop.mass)); 274 PetscCall(VecDestroy(&appctx.dat.curr_sol)); 275 PetscCall(PetscFree2(appctx.SEMop.gll.nodes,appctx.SEMop.gll.weights)); 276 PetscCall(DMDestroy(&appctx.da)); 277 PetscCall(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 PetscCall(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 PetscCall(DMDAVecGetArray(appctx->da,u,&s)); 310 PetscCall(DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg)); 311 PetscCall(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 PetscCall(DMDAVecRestoreArray(appctx->da,u,&s)); 316 PetscCall(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 PetscCall(DMDAVecGetArray(appctx->da,u,&s)); 340 PetscCall(DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg)); 341 PetscCall(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 PetscCall(DMDAVecRestoreArray(appctx->da,u,&s)); 346 PetscCall(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 PetscCall(DMDAVecGetArray(appctx->da,obj,&s)); 367 PetscCall(DMDAVecGetArrayRead(appctx->da,appctx->SEMop.grid,(void*)&xg)); 368 PetscCall(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 PetscCall(DMDAVecRestoreArray(appctx->da,obj,&s)); 374 PetscCall(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 PetscCall(MatMult(appctx->SEMop.grad,globalin,globalout)); /* grad u */ 384 PetscCall(VecPointwiseMult(globalout,globalin,globalout)); /* u grad u */ 385 PetscCall(VecScale(globalout, -1.0)); 386 PetscCall(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 PetscCall(MatCopy(appctx->SEMop.grad,A,SAME_NONZERO_PATTERN)); 407 PetscCall(MatDiagonalScale(A,globalin,NULL)); 408 409 /* A = A + diag(Gu) */ 410 PetscCall(VecDuplicate(globalin,&Gglobalin)); 411 PetscCall(MatMult(appctx->SEMop.grad,globalin,Gglobalin)); 412 PetscCall(MatDiagonalSet(A,Gglobalin,ADD_VALUES)); 413 PetscCall(VecDestroy(&Gglobalin)); 414 415 /* A = K - A */ 416 PetscCall(MatScale(A,-1.0)); 417 PetscCall(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 PetscCall(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 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 462 PetscCall(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 PetscCall(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 PetscCall(MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES)); 476 } 477 PetscCall(PetscFree(rowsDM)); 478 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 479 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 480 PetscCall(VecReciprocal(appctx->SEMop.mass)); 481 PetscCall(MatDiagonalScale(A,appctx->SEMop.mass,0)); 482 PetscCall(VecReciprocal(appctx->SEMop.mass)); 483 484 PetscCall(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 PetscCall(PetscGaussLobattoLegendreElementAdvectionCreate(appctx->SEMop.gll.n,appctx->SEMop.gll.nodes,appctx->SEMop.gll.weights,&temp)); 516 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 517 518 PetscCall(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 PetscCall(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 PetscCall(MatSetValuesLocal(A,appctx->param.N,rowsDM,appctx->param.N,rowsDM,&temp[0][0],ADD_VALUES)); 529 } 530 PetscCall(PetscFree(rowsDM)); 531 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 532 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 533 534 PetscCall(VecReciprocal(appctx->SEMop.mass)); 535 PetscCall(MatDiagonalScale(A,appctx->SEMop.mass,0)); 536 PetscCall(VecReciprocal(appctx->SEMop.mass)); 537 PetscCall(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 PetscCall(TSSetTime(appctx->ts,0.0)); 585 PetscCall(TSSetStepNumber(appctx->ts,0)); 586 PetscCall(TSSetTimeStep(appctx->ts,appctx->initial_dt)); 587 PetscCall(VecCopy(IC,appctx->dat.curr_sol)); 588 589 PetscCall(TSSolve(appctx->ts,appctx->dat.curr_sol)); 590 591 PetscCall(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 PetscCall(VecDuplicate(G,&temp)); 597 PetscCall(VecPointwiseMult(temp,G,G)); 598 PetscCall(VecDot(temp,appctx->SEMop.mass,f)); 599 600 /* local error evaluation */ 601 PetscCall(VecWAXPY(temp,-1.0,appctx->dat.ic,appctx->dat.true_solution)); 602 PetscCall(VecPointwiseMult(temp,temp,temp)); 603 /* for error evaluation */ 604 PetscCall(VecDot(temp,appctx->SEMop.mass,&errex)); 605 PetscCall(VecDestroy(&temp)); 606 errex = PetscSqrtReal(errex); 607 608 /* 609 Compute initial conditions for the adjoint integration. See Notes above 610 */ 611 612 PetscCall(VecScale(G, -2.0)); 613 PetscCall(VecPointwiseMult(G,G,appctx->SEMop.mass)); 614 PetscCall(TSSetCostGradients(appctx->ts,1,&G,NULL)); 615 PetscCall(TSAdjointSolve(appctx->ts)); 616 PetscCall(VecPointwiseDivide(G,G,appctx->SEMop.mass)); 617 618 PetscCall(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 PetscCall(VecDuplicate(appctx->dat.ic,&temp)); 630 PetscCall(VecWAXPY(temp,-1.0,appctx->dat.ic,appctx->dat.true_solution)); 631 PetscCall(VecPointwiseMult(temp,temp,temp)); 632 PetscCall(VecDot(temp,appctx->SEMop.mass,&nrm)); 633 PetscCall(VecDestroy(&temp)); 634 nrm = PetscSqrtReal(nrm); 635 PetscCall(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