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