HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
SolvePETSc.cpp
Go to the documentation of this file.
1 
9 #ifdef with_petsc
10 
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <math.h>
14 #include <string.h>
15 #include <basic.h>
16 #include <io_cpp.h>
17 #include <petscinterface.h>
18 #include <simulation_object.h>
19 
20 #ifdef with_librom
21 #include <librom_interface.h>
22 #endif
23 
24 #undef __FUNCT__
25 #define __FUNCT__ "SolvePETSc"
26 
27 extern "C" int CalculateError (void*,void*);
28 int OutputSolution (void*,int,double);
29 extern "C" void ResetFilenameIndex(char*, int);
30 #ifdef with_librom
31 extern "C" int CalculateROMDiff(void*,void*);
32 int OutputROMSolution(void*,int,double);
33 #endif
34 
50 int SolvePETSc( void* s,
51  int nsims,
52  int rank,
53  int nproc )
54 {
55  SimulationObject* sim = (SimulationObject*) s;
56 
57  DM dm; /* data management object */
58  TS ts; /* time integration object */
59  Vec Y,Z; /* PETSc solution vectors */
60  Mat A, B; /* Jacobian and preconditioning matrices */
61  MatFDColoring fdcoloring; /* coloring for sparse Jacobian computation */
62  TSType time_scheme; /* time integration method */
63  TSProblemType ptype; /* problem type - nonlinear or linear */
64 
65  int flag_mat_a = 0,
66  flag_mat_b = 0,
67  flag_fdcoloring = 0,
68  iAuxSize = 0, i;
69 
70  PetscFunctionBegin;
71 
72  /* Register custom time-integration methods, if specified */
73  PetscRegisterTIMethods(rank);
74  if(!rank) printf("Setting up PETSc time integration... \n");
75 
76  /* create and set a PETSc context */
77  PETScContext context;
78 
79  context.rank = rank;
80  context.nproc = nproc;
81 
82  context.simobj = sim;
83  context.nsims = nsims;
84 
85  /* default: everything explicit */
86  context.flag_hyperbolic = _EXPLICIT_;
87  context.flag_hyperbolic_f = _EXPLICIT_;
88  context.flag_hyperbolic_df = _EXPLICIT_;
89  context.flag_parabolic = _EXPLICIT_;
90  context.flag_source = _EXPLICIT_;
91 
92  context.tic = 0;
93  context.flag_is_linear = 0;
94  context.globalDOF.clear();
95  context.points.clear();
96  context.ti_runtime = 0.0;
97  context.waqt = 0.0;
98  context.dt = sim[0].solver.dt;
99  context.stage_times.clear();
100  context.stage_index = 0;
101 
102 #ifdef with_librom
103  if (!rank) printf("Setting up libROM interface.\n");
104  context.rom_interface = new libROMInterface( sim,
105  nsims,
106  rank,
107  nproc,
108  sim[0].solver.dt );
109  context.rom_mode = ((libROMInterface*)context.rom_interface)->mode();
110  context.op_times_arr.clear();
111 #endif
112 
113 #ifdef with_librom
114  if ( (context.rom_mode == _ROM_MODE_TRAIN_)
115  || (context.rom_mode == _ROM_MODE_INITIAL_GUESS_ )
116  || (context.rom_mode == _ROM_MODE_NONE_ ) ) {
117 
118  if (context.rom_mode == _ROM_MODE_INITIAL_GUESS_) {
119  ((libROMInterface*)context.rom_interface)->loadROM();
120  ((libROMInterface*)context.rom_interface)->projectInitialSolution(sim);
121  }
122 #endif
123  PetscCreatePointList(&context);
124 
125  /* create and initialize PETSc solution vector and other parameters */
126  /* PETSc solution vector does not have ghost points */
127  VecCreate(MPI_COMM_WORLD,&Y);
128  VecSetSizes(Y,context.ndofs,PETSC_DECIDE);
129  VecSetUp(Y);
130 
131  /* copy initial solution to PETSc's vector */
132  for (int ns = 0; ns < nsims; ns++) {
133  TransferVecToPETSc( sim[ns].solver.u,
134  Y,
135  &context,
136  ns,
137  context.offsets[ns] );
138  }
139 
140  /* Create the global DOF mapping for all the grid points */
141  PetscGlobalDOF(&context);
142 
143  /* Define and initialize the time-integration object */
144  TSCreate(MPI_COMM_WORLD,&ts);
145  TSSetMaxSteps(ts,sim[0].solver.n_iter);
146  TSSetMaxTime(ts,sim[0].solver.dt*sim[0].solver.n_iter);
147  TSSetTimeStep(ts,sim[0].solver.dt);
148  TSSetTime(ts,context.waqt);
149  TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
150  TSSetType(ts,TSBEULER);
151 
152  /* set default time step adaptivity to none */
153  TSAdapt adapt;
154  TSAdaptType adapt_type = TSADAPTNONE;
155  TSGetAdapt(ts,&adapt);
156  TSAdaptSetType(adapt,adapt_type);
157 
158  /* set options from input */
159  TSSetFromOptions(ts);
160 
161  /* create DM */
162  DMShellCreate(MPI_COMM_WORLD, &dm);
163  DMShellSetGlobalVector(dm, Y);
164  TSSetDM(ts, dm);
165 
166 #ifdef with_librom
167  TSAdaptGetType(adapt,&adapt_type);
168  if (strcmp(adapt_type, TSADAPTNONE)) {
169  if (!rank) printf("Warning: libROM interface not yet implemented for adaptive timestepping.\n");
170  }
171 #endif
172 
173  /* Define the right and left -hand side functions for each time-integration scheme */
174  TSGetType(ts,&time_scheme);
175  TSGetProblemType(ts,&ptype);
176 
177  if (!strcmp(time_scheme,TSARKIMEX)) {
178 
179  /* implicit - explicit time integration */
180 
181  TSSetRHSFunction(ts,nullptr,PetscRHSFunctionIMEX,&context);
182  TSSetIFunction (ts,nullptr,PetscIFunctionIMEX, &context);
183 
184  SNES snes;
185  KSP ksp;
186  PC pc;
187  SNESType snestype;
188  TSGetSNES(ts,&snes);
189  SNESGetType(snes,&snestype);
190 
191 #ifdef with_librom
192  if (context.rom_mode == _ROM_MODE_INITIAL_GUESS_) {
193  SNESSetComputeInitialGuess(snes, PetscSetInitialGuessROM, &context);
194  }
195 #endif
196 
197  context.flag_use_precon = 0;
198  PetscOptionsGetBool( nullptr,nullptr,
199  "-with_pc",
200  (PetscBool*)(&context.flag_use_precon),
201  nullptr );
202 
203  char precon_mat_type_c_st[_MAX_STRING_SIZE_] = "default";
204  PetscOptionsGetString( nullptr,
205  nullptr,
206  "-pc_matrix_type",
207  precon_mat_type_c_st,
208  _MAX_STRING_SIZE_,
209  nullptr );
210  context.precon_matrix_type = std::string(precon_mat_type_c_st);
211 
212  if (context.flag_use_precon) {
213 
214  if (context.precon_matrix_type == "default") {
215 
216  /* Matrix-free representation of the Jacobian */
217  flag_mat_a = 1;
218  MatCreateShell( MPI_COMM_WORLD,
219  context.ndofs,
220  context.ndofs,
221  PETSC_DETERMINE,
222  PETSC_DETERMINE,
223  &context,
224  &A);
225  if ((!strcmp(snestype,SNESKSPONLY)) || (ptype == TS_LINEAR)) {
226  /* linear problem */
227  context.flag_is_linear = 1;
228  MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunctionIMEX_Linear);
229  SNESSetType(snes,SNESKSPONLY);
230  } else {
231  /* nonlinear problem */
232  context.flag_is_linear = 0;
233  context.jfnk_eps = 1e-7;
234  PetscOptionsGetReal(NULL,NULL,"-jfnk_epsilon",&context.jfnk_eps,NULL);
235  MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunctionIMEX_JFNK);
236  }
237  MatSetUp(A);
238  /* check if Jacobian of the physical model is defined */
239  for (int ns = 0; ns < nsims; ns++) {
240  if ((!sim[ns].solver.JFunction) && (!sim[ns].solver.KFunction)) {
241  if (!rank) {
242  fprintf(stderr,"Error in SolvePETSc(): solver->JFunction or solver->KFunction ");
243  fprintf(stderr,"(point-wise Jacobians for hyperbolic or parabolic terms) must ");
244  fprintf(stderr,"be defined for preconditioning.\n");
245  }
246  PetscFunctionReturn(1);
247  }
248  }
249  /* Set up preconditioner matrix */
250  flag_mat_b = 1;
251  MatCreateAIJ( MPI_COMM_WORLD,
252  context.ndofs,
253  context.ndofs,
254  PETSC_DETERMINE,
255  PETSC_DETERMINE,
256  (sim[0].solver.ndims*2+1)*sim[0].solver.nvars, NULL,
257  2*sim[0].solver.ndims*sim[0].solver.nvars, NULL,
258  &B );
259  MatSetBlockSize(B,sim[0].solver.nvars);
260  /* Set the IJacobian function for TS */
261  TSSetIJacobian(ts,A,B,PetscIJacobianIMEX,&context);
262 
263  } else if (context.precon_matrix_type == "fd") {
264 
265  flag_mat_a = 1;
266  MatCreateSNESMF(snes,&A);
267  flag_mat_b = 1;
268  MatCreateAIJ( MPI_COMM_WORLD,
269  context.ndofs,
270  context.ndofs,
271  PETSC_DETERMINE,
272  PETSC_DETERMINE,
273  (sim[0].solver.ndims*2+1)*sim[0].solver.nvars, NULL,
274  2*sim[0].solver.ndims*sim[0].solver.nvars, NULL,
275  &B);
276  MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
277  /* Set the Jacobian function for SNES */
278  SNESSetJacobian(snes, A, B, SNESComputeJacobianDefault, NULL);
279 
280  } else if (context.precon_matrix_type == "colored_fd") {
281 
282  int stencil_width = 1;
283  PetscOptionsGetInt( NULL,
284  NULL,
285  "-pc_matrix_colored_fd_stencil_width",
286  &stencil_width,
287  NULL );
288 
289  flag_mat_a = 1;
290  MatCreateSNESMF(snes,&A);
291  flag_mat_b = 1;
292  MatCreateAIJ( MPI_COMM_WORLD,
293  context.ndofs,
294  context.ndofs,
295  PETSC_DETERMINE,
296  PETSC_DETERMINE,
297  (sim[0].solver.ndims*2+1)*sim[0].solver.nvars, NULL,
298  2*sim[0].solver.ndims*sim[0].solver.nvars, NULL,
299  &B);
300  MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
301  if (!rank) {
302  printf("PETSc: Setting Jacobian non-zero pattern (stencil width %d).\n",
303  stencil_width );
304  }
305  PetscJacobianMatNonzeroEntriesImpl(B, stencil_width, &context);
306 
307  /* Set the Jacobian function for SNES */
308  SNESSetJacobian(snes, A, B, SNESComputeJacobianDefaultColor, NULL);
309 
310  } else {
311 
312  if (!rank) {
313  fprintf( stderr,"Invalid input for \"-pc_matrix_type\": %s.\n",
314  context.precon_matrix_type.c_str());
315  }
316  PetscFunctionReturn(0);
317 
318  }
319 
320  /* set PC side to right */
321  SNESGetKSP(snes,&ksp);
322  KSPSetPCSide(ksp, PC_RIGHT);
323 
324  } else {
325 
326  /* Matrix-free representation of the Jacobian */
327  flag_mat_a = 1;
328  MatCreateShell( MPI_COMM_WORLD,
329  context.ndofs,
330  context.ndofs,
331  PETSC_DETERMINE,
332  PETSC_DETERMINE,
333  &context,
334  &A);
335  if ((!strcmp(snestype,SNESKSPONLY)) || (ptype == TS_LINEAR)) {
336  /* linear problem */
337  context.flag_is_linear = 1;
338  MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunctionIMEX_Linear);
339  SNESSetType(snes,SNESKSPONLY);
340  } else {
341  /* nonlinear problem */
342  context.flag_is_linear = 0;
343  context.jfnk_eps = 1e-7;
344  PetscOptionsGetReal(NULL,NULL,"-jfnk_epsilon",&context.jfnk_eps,NULL);
345  MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunctionIMEX_JFNK);
346  }
347  MatSetUp(A);
348  /* Set the RHSJacobian function for TS */
349  TSSetIJacobian(ts,A,A,PetscIJacobianIMEX,&context);
350  /* Set PC (preconditioner) to none */
351  SNESGetKSP(snes,&ksp);
352  KSPGetPC(ksp,&pc);
353  PCSetType(pc,PCNONE);
354  }
355 
356  /* read the implicit/explicit flags for each of the terms for IMEX schemes */
357  /* default -> hyperbolic - explicit, parabolic and source - implicit */
358  PetscBool flag = PETSC_FALSE;
359 
360  context.flag_hyperbolic = _EXPLICIT_;
361  context.flag_hyperbolic_f = _EXPLICIT_;
362  context.flag_hyperbolic_df = _IMPLICIT_;
363  context.flag_parabolic = _IMPLICIT_;
364  context.flag_source = _IMPLICIT_;
365 
366  if (!strcmp(sim[0].solver.SplitHyperbolicFlux,"yes")) {
367 
368  flag = PETSC_FALSE;
369  PetscOptionsGetBool(nullptr,nullptr,"-hyperbolic_f_explicit",&flag,nullptr);
370  if (flag == PETSC_TRUE) context.flag_hyperbolic_f = _EXPLICIT_;
371  flag = PETSC_FALSE;
372  PetscOptionsGetBool(nullptr,nullptr,"-hyperbolic_f_implicit",&flag,nullptr);
373  if (flag == PETSC_TRUE) context.flag_hyperbolic_f = _IMPLICIT_;
374 
375  flag = PETSC_FALSE;
376  PetscOptionsGetBool(nullptr,nullptr,"-hyperbolic_df_explicit",&flag,nullptr);
377  if (flag == PETSC_TRUE) context.flag_hyperbolic_df = _EXPLICIT_;
378  flag = PETSC_FALSE;
379  PetscOptionsGetBool(nullptr,nullptr,"-hyperbolic_df_implicit",&flag,nullptr);
380  if (flag == PETSC_TRUE) context.flag_hyperbolic_df = _IMPLICIT_;
381 
382  } else {
383 
384  flag = PETSC_FALSE;
385  PetscOptionsGetBool(nullptr,nullptr,"-hyperbolic_explicit",&flag,nullptr);
386  if (flag == PETSC_TRUE) context.flag_hyperbolic = _EXPLICIT_;
387  flag = PETSC_FALSE;
388  PetscOptionsGetBool(nullptr,nullptr,"-hyperbolic_implicit",&flag,nullptr);
389  if (flag == PETSC_TRUE) context.flag_hyperbolic = _IMPLICIT_;
390 
391  }
392 
393  flag = PETSC_FALSE;
394  PetscOptionsGetBool(nullptr,nullptr,"-parabolic_explicit",&flag,nullptr);
395  if (flag == PETSC_TRUE) context.flag_parabolic = _EXPLICIT_;
396  flag = PETSC_FALSE;
397  PetscOptionsGetBool(nullptr,nullptr,"-parabolic_implicit",&flag,nullptr);
398  if (flag == PETSC_TRUE) context.flag_parabolic = _IMPLICIT_;
399 
400  flag = PETSC_FALSE;
401  PetscOptionsGetBool(nullptr,nullptr,"-source_explicit",&flag,nullptr);
402  if (flag == PETSC_TRUE) context.flag_source = _EXPLICIT_;
403  flag = PETSC_FALSE;
404  PetscOptionsGetBool(nullptr,nullptr,"-source_implicit",&flag,nullptr);
405  if (flag == PETSC_TRUE) context.flag_source = _IMPLICIT_;
406 
407  flag = PETSC_FALSE;
408  PetscOptionsGetBool(nullptr,nullptr,"-ts_arkimex_fully_implicit",&flag,nullptr);
409  if (flag == PETSC_TRUE) {
410  context.flag_hyperbolic_f = _IMPLICIT_;
411  context.flag_hyperbolic_df = _IMPLICIT_;
412  context.flag_hyperbolic = _IMPLICIT_;
413  context.flag_parabolic = _IMPLICIT_;
414  context.flag_source = _IMPLICIT_;
415  }
416 
417  /* print out a summary of the treatment of each term */
418  if (!rank) {
419  printf("Implicit-Explicit time-integration:-\n");
420  if (!strcmp(sim[0].solver.SplitHyperbolicFlux,"yes")) {
421  if (context.flag_hyperbolic_f == _EXPLICIT_) printf("Hyperbolic (f-df) term: Explicit\n");
422  else printf("Hyperbolic (f-df) term: Implicit\n");
423  if (context.flag_hyperbolic_df == _EXPLICIT_) printf("Hyperbolic (df) term: Explicit\n");
424  else printf("Hyperbolic (df) term: Implicit\n");
425  } else {
426  if (context.flag_hyperbolic == _EXPLICIT_) printf("Hyperbolic term: Explicit\n");
427  else printf("Hyperbolic term: Implicit\n");
428  }
429  if (context.flag_parabolic == _EXPLICIT_) printf("Parabolic term: Explicit\n");
430  else printf("Parabolic term: Implicit\n");
431  if (context.flag_source == _EXPLICIT_) printf("Source term: Explicit\n");
432  else printf("Source term: Implicit\n");
433  }
434 
435  } else if ( (!strcmp(time_scheme,TSEULER))
436  || (!strcmp(time_scheme,TSRK ))
437  || (!strcmp(time_scheme,TSSSP )) ) {
438 
439  /* Explicit time integration */
440  TSSetRHSFunction(ts,nullptr,PetscRHSFunctionExpl,&context);
441 
442  } else if ( (!strcmp(time_scheme,TSCN))
443  || (!strcmp(time_scheme,TSBEULER )) ) {
444 
445 
446  /* Implicit time integration */
447 
448  TSSetIFunction(ts,nullptr,PetscIFunctionImpl,&context);
449 
450  SNES snes;
451  KSP ksp;
452  PC pc;
453  SNESType snestype;
454  TSGetSNES(ts,&snes);
455  SNESGetType(snes,&snestype);
456 
457 #ifdef with_librom
458  if (context.rom_mode == _ROM_MODE_INITIAL_GUESS_) {
459  SNESSetComputeInitialGuess(snes, PetscSetInitialGuessROM, &context);
460  }
461 #endif
462 
463  context.flag_use_precon = 0;
464  PetscOptionsGetBool( nullptr,
465  nullptr,
466  "-with_pc",
467  (PetscBool*)(&context.flag_use_precon),
468  nullptr );
469 
470  char precon_mat_type_c_st[_MAX_STRING_SIZE_] = "default";
471  PetscOptionsGetString( nullptr,
472  nullptr,
473  "-pc_matrix_type",
474  precon_mat_type_c_st,
475  _MAX_STRING_SIZE_,
476  nullptr );
477  context.precon_matrix_type = std::string(precon_mat_type_c_st);
478 
479  if (context.flag_use_precon) {
480 
481  if (context.precon_matrix_type == "default") {
482 
483  /* Matrix-free representation of the Jacobian */
484  flag_mat_a = 1;
485  MatCreateShell( MPI_COMM_WORLD,
486  context.ndofs,
487  context.ndofs,
488  PETSC_DETERMINE,
489  PETSC_DETERMINE,
490  &context,
491  &A);
492  if ((!strcmp(snestype,SNESKSPONLY)) || (ptype == TS_LINEAR)) {
493  /* linear problem */
494  context.flag_is_linear = 1;
495  MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunction_Linear);
496  SNESSetType(snes,SNESKSPONLY);
497  } else {
498  /* nonlinear problem */
499  context.flag_is_linear = 0;
500  context.jfnk_eps = 1e-7;
501  PetscOptionsGetReal(NULL,NULL,"-jfnk_epsilon",&context.jfnk_eps,NULL);
502  MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunction_JFNK);
503  }
504  MatSetUp(A);
505  /* check if Jacobian of the physical model is defined */
506  for (int ns = 0; ns < nsims; ns++) {
507  if ((!sim[ns].solver.JFunction) && (!sim[ns].solver.KFunction)) {
508  if (!rank) {
509  fprintf(stderr,"Error in SolvePETSc(): solver->JFunction or solver->KFunction ");
510  fprintf(stderr,"(point-wise Jacobians for hyperbolic or parabolic terms) must ");
511  fprintf(stderr,"be defined for preconditioning.\n");
512  }
513  PetscFunctionReturn(1);
514  }
515  }
516  /* Set up preconditioner matrix */
517  flag_mat_b = 1;
518  MatCreateAIJ( MPI_COMM_WORLD,
519  context.ndofs,
520  context.ndofs,
521  PETSC_DETERMINE,
522  PETSC_DETERMINE,
523  (sim[0].solver.ndims*2+1)*sim[0].solver.nvars, NULL,
524  2*sim[0].solver.ndims*sim[0].solver.nvars, NULL,
525  &B );
526  MatSetBlockSize(B,sim[0].solver.nvars);
527  /* Set the IJacobian function for TS */
528  TSSetIJacobian(ts,A,B,PetscIJacobian,&context);
529 
530  } else if (context.precon_matrix_type == "fd") {
531 
532  flag_mat_a = 1;
533  MatCreateSNESMF(snes,&A);
534  flag_mat_b = 1;
535  MatCreateAIJ( MPI_COMM_WORLD,
536  context.ndofs,
537  context.ndofs,
538  PETSC_DETERMINE,
539  PETSC_DETERMINE,
540  (sim[0].solver.ndims*2+1)*sim[0].solver.nvars, NULL,
541  2*sim[0].solver.ndims*sim[0].solver.nvars, NULL,
542  &B);
543  MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
544  /* Set the Jacobian function for SNES */
545  SNESSetJacobian(snes, A, B, SNESComputeJacobianDefault, NULL);
546 
547  } else if (context.precon_matrix_type == "colored_fd") {
548 
549  int stencil_width = 1;
550  PetscOptionsGetInt( NULL,
551  NULL,
552  "-pc_matrix_colored_fd_stencil_width",
553  &stencil_width,
554  NULL );
555 
556  flag_mat_a = 1;
557  MatCreateSNESMF(snes,&A);
558  flag_mat_b = 1;
559  MatCreateAIJ( MPI_COMM_WORLD,
560  context.ndofs,
561  context.ndofs,
562  PETSC_DETERMINE,
563  PETSC_DETERMINE,
564  (sim[0].solver.ndims*2+1)*sim[0].solver.nvars, NULL,
565  2*sim[0].solver.ndims*sim[0].solver.nvars, NULL,
566  &B);
567  MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
568  if (!rank) {
569  printf("PETSc: Setting Jacobian non-zero pattern (stencil width %d).\n",
570  stencil_width );
571  }
572  PetscJacobianMatNonzeroEntriesImpl(B, stencil_width, &context);
573 
574  /* Set the Jacobian function for SNES */
575  SNESSetJacobian(snes, A, B, SNESComputeJacobianDefaultColor, NULL);
576 
577  } else {
578 
579  if (!rank) {
580  fprintf( stderr,"Invalid input for \"-pc_matrix_type\": %s.\n",
581  context.precon_matrix_type.c_str());
582  }
583  PetscFunctionReturn(0);
584 
585  }
586 
587  /* set PC side to right */
588  SNESGetKSP(snes,&ksp);
589  KSPSetPCSide(ksp, PC_RIGHT);
590 
591  } else {
592 
593  /* Matrix-free representation of the Jacobian */
594  flag_mat_a = 1;
595  MatCreateShell( MPI_COMM_WORLD,
596  context.ndofs,
597  context.ndofs,
598  PETSC_DETERMINE,
599  PETSC_DETERMINE,
600  &context,
601  &A);
602  if ((!strcmp(snestype,SNESKSPONLY)) || (ptype == TS_LINEAR)) {
603  /* linear problem */
604  context.flag_is_linear = 1;
605  MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunction_Linear);
606  SNESSetType(snes,SNESKSPONLY);
607  } else {
608  /* nonlinear problem */
609  context.flag_is_linear = 0;
610  context.jfnk_eps = 1e-7;
611  PetscOptionsGetReal(NULL,NULL,"-jfnk_epsilon",&context.jfnk_eps,NULL);
612  MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscJacobianFunction_JFNK);
613  }
614  MatSetUp(A);
615  /* Set the RHSJacobian function for TS */
616  TSSetIJacobian(ts,A,A,PetscIJacobian,&context);
617  /* Set PC (preconditioner) to none */
618  SNESGetKSP(snes,&ksp);
619  KSPGetPC(ksp,&pc);
620  PCSetType(pc,PCNONE);
621  }
622 
623  } else {
624 
625  if (!rank) {
626  fprintf(stderr, "Time integration type %s is not yet supported.\n", time_scheme);
627  }
628  PetscFunctionReturn(0);
629 
630  }
631 
632  /* Set pre/post-stage and post-timestep function */
633  TSSetPreStep (ts,PetscPreTimeStep );
634  TSSetPreStage(ts,PetscPreStage );
635  TSSetPostStage(ts,PetscPostStage );
636  TSSetPostStep(ts,PetscPostTimeStep);
637  /* Set solution vector for TS */
638  TSSetSolution(ts,Y);
639  /* Set it all up */
640  TSSetUp(ts);
641  /* Set application context */
642  TSSetApplicationContext(ts,&context);
643 
644  if (!rank) {
645  if (context.flag_is_linear) printf("SolvePETSc(): Problem type is linear.\n");
646  else printf("SolvePETSc(): Problem type is nonlinear.\n");
647  }
648 
649  if (!rank) printf("** Starting PETSc time integration **\n");
650  context.ti_runtime = 0.0;
651  TSSolve(ts,Y);
652  if (!rank) {
653  printf("** Completed PETSc time integration (Final time: %f), total wctime: %f (seconds) **\n",
654  context.waqt, context.ti_runtime );
655  }
656 
657  /* Get the number of time steps */
658  for (int ns = 0; ns < nsims; ns++) {
659  TSGetStepNumber(ts,&(sim[ns].solver.n_iter));
660  }
661 
662  /* get and write to file any auxiliary solutions */
663  char aux_fname_root[4] = "ts0";
664  TSGetSolutionComponents(ts,&iAuxSize,NULL);
665  if (iAuxSize) {
666  if (iAuxSize > 10) iAuxSize = 10;
667  if (!rank) printf("Number of auxiliary solutions from time integration: %d\n",iAuxSize);
668  VecDuplicate(Y,&Z);
669  for (i=0; i<iAuxSize; i++) {
670  TSGetSolutionComponents(ts,&i,&Z);
671  for (int ns = 0; ns < nsims; ns++) {
672  TransferVecFromPETSc(sim[ns].solver.u,Z,&context,ns,context.offsets[ns]);
673  WriteArray( sim[ns].solver.ndims,
674  sim[ns].solver.nvars,
675  sim[ns].solver.dim_global,
676  sim[ns].solver.dim_local,
677  sim[ns].solver.ghosts,
678  sim[ns].solver.x,
679  sim[ns].solver.u,
680  &(sim[ns].solver),
681  &(sim[ns].mpi),
682  aux_fname_root );
683  }
684  aux_fname_root[2]++;
685  }
686  VecDestroy(&Z);
687  }
688 
689  /* if available, get error estimates */
690  PetscTimeError(ts);
691 
692  /* copy final solution from PETSc's vector */
693  for (int ns = 0; ns < nsims; ns++) {
694  TransferVecFromPETSc(sim[ns].solver.u,Y,&context,ns,context.offsets[ns]);
695  }
696 
697  /* clean up */
698  VecDestroy(&Y);
699  if (flag_mat_a) { MatDestroy(&A); }
700  if (flag_mat_b) { MatDestroy(&B); }
701  if (flag_fdcoloring) { MatFDColoringDestroy(&fdcoloring); }
702  TSDestroy(&ts);
703  DMDestroy(&dm);
704 
705  /* write a final solution file, if last iteration did not write one */
706  if (context.tic) {
707  for (int ns = 0; ns < nsims; ns++) {
708  HyPar* solver = &(sim[ns].solver);
709  MPIVariables* mpi = &(sim[ns].mpi);
710  if (solver->PhysicsOutput) {
711  solver->PhysicsOutput(solver,mpi, context.waqt);
712  }
713  CalculateError(solver,mpi);
714  }
715  OutputSolution(sim, nsims, context.waqt);
716  }
717  /* calculate error if exact solution has been provided */
718  for (int ns = 0; ns < nsims; ns++) {
719  CalculateError(&(sim[ns].solver), &(sim[ns].mpi));
720  }
721 
722  PetscCleanup(&context);
723 
724 #ifdef with_librom
725  context.op_times_arr.push_back(context.waqt);
726 
727  for (int ns = 0; ns < nsims; ns++) {
728  ResetFilenameIndex( sim[ns].solver.filename_index,
729  sim[ns].solver.index_length );
730  }
731 
732  if (((libROMInterface*)context.rom_interface)->mode() == _ROM_MODE_TRAIN_) {
733 
734  ((libROMInterface*)context.rom_interface)->train();
735  if (!rank) printf("libROM: total training wallclock time: %f (seconds).\n",
736  ((libROMInterface*)context.rom_interface)->trainWallclockTime() );
737 
738  double total_rom_predict_time = 0;
739  for (int iter = 0; iter < context.op_times_arr.size(); iter++) {
740 
741  double waqt = context.op_times_arr[iter];
742 
743  ((libROMInterface*)context.rom_interface)->predict(sim, waqt);
744  if (!rank) printf( "libROM: Predicted solution at time %1.4e using ROM, wallclock time: %f.\n",
745  waqt, ((libROMInterface*)context.rom_interface)->predictWallclockTime() );
746  total_rom_predict_time += ((libROMInterface*)context.rom_interface)->predictWallclockTime();
747 
748  /* calculate diff between ROM and PDE solutions */
749  if (iter == (context.op_times_arr.size()-1)) {
750  if (!rank) printf("libROM: Calculating diff between PDE and ROM solutions.\n");
751  for (int ns = 0; ns < nsims; ns++) {
752  CalculateROMDiff( &(sim[ns].solver),
753  &(sim[ns].mpi) );
754  }
755  }
756  /* write the ROM solution to file */
757  OutputROMSolution(sim, nsims, waqt);
758 
759  }
760 
761  if (!rank) {
762  printf( "libROM: total prediction/query wallclock time: %f (seconds).\n",
763  total_rom_predict_time );
764  }
765 
766  ((libROMInterface*)context.rom_interface)->saveROM();
767 
768  } else {
769 
770  for (int ns = 0; ns < nsims; ns++) {
771  sim[ns].solver.rom_diff_norms[0]
772  = sim[ns].solver.rom_diff_norms[1]
773  = sim[ns].solver.rom_diff_norms[2]
774  = -1;
775  }
776 
777  }
778 
779  } else if (context.rom_mode == _ROM_MODE_PREDICT_) {
780 
781  for (int ns = 0; ns < nsims; ns++) {
782  sim[ns].solver.rom_diff_norms[0]
783  = sim[ns].solver.rom_diff_norms[1]
784  = sim[ns].solver.rom_diff_norms[2]
785  = -1;
786  strcpy(sim[ns].solver.ConservationCheck,"no");
787  }
788 
789  ((libROMInterface*)context.rom_interface)->loadROM();
790  ((libROMInterface*)context.rom_interface)->projectInitialSolution(sim);
791 
792  {
793  int start_iter = sim[0].solver.restart_iter;
794  int n_iter = sim[0].solver.n_iter;
795  double dt = sim[0].solver.dt;
796 
797  double cur_time = start_iter * dt;
798  context.op_times_arr.push_back(cur_time);
799 
800  for (int iter = start_iter; iter < n_iter; iter++) {
801  cur_time += dt;
802  if ( ( (iter+1)%sim[0].solver.file_op_iter == 0)
803  && ( (iter+1) < n_iter) ) {
804  context.op_times_arr.push_back(cur_time);
805  }
806  }
807 
808  double t_final = n_iter*dt;
809  context.op_times_arr.push_back(t_final);
810  }
811 
812  double total_rom_predict_time = 0;
813  for (int iter = 0; iter < context.op_times_arr.size(); iter++) {
814 
815  double waqt = context.op_times_arr[iter];
816 
817  ((libROMInterface*)context.rom_interface)->predict(sim, waqt);
818  if (!rank) printf( "libROM: Predicted solution at time %1.4e using ROM, wallclock time: %f.\n",
819  waqt, ((libROMInterface*)context.rom_interface)->predictWallclockTime() );
820  total_rom_predict_time += ((libROMInterface*)context.rom_interface)->predictWallclockTime();
821 
822  /* write the solution to file */
823  for (int ns = 0; ns < nsims; ns++) {
824  if (sim[ns].solver.PhysicsOutput) {
825  sim[ns].solver.PhysicsOutput( &(sim[ns].solver),
826  &(sim[ns].mpi),
827  waqt );
828  }
829  }
830  OutputSolution(sim, nsims, waqt);
831 
832  }
833 
834  /* calculate error if exact solution has been provided */
835  for (int ns = 0; ns < nsims; ns++) {
836  CalculateError(&(sim[ns].solver),
837  &(sim[ns].mpi) );
838  }
839 
840  if (!rank) {
841  printf( "libROM: total prediction/query wallclock time: %f (seconds).\n",
842  total_rom_predict_time );
843  }
844 
845  }
846 
847  delete ((libROMInterface*)context.rom_interface);
848 #endif
849 
850  PetscFunctionReturn(0);
851 }
852 
853 #endif
CalculateError
int CalculateError(void *s, void *m)
Definition: CalculateError.c:25
PetscTimeError
PetscErrorCode PetscTimeError(TS)
Definition: PetscError.cpp:20
PetscIJacobianIMEX
PetscErrorCode PetscIJacobianIMEX(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *)
Definition: PetscIJacobianIMEX.cpp:61
PetscJacobianMatNonzeroEntriesImpl
int PetscJacobianMatNonzeroEntriesImpl(Mat, int, void *)
Definition: PetscJacobianMatNonzeroEntriesImpl.cpp:38
PetscJacobianFunction_JFNK
PetscErrorCode PetscJacobianFunction_JFNK(Mat, Vec, Vec)
Definition: PetscIJacobian.cpp:111
io_cpp.h
Function declarations for file I/O functions.
TransferVecFromPETSc
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
Definition: TransferFromPETSc.cpp:23
PetscJacobianFunctionIMEX_JFNK
PetscErrorCode PetscJacobianFunctionIMEX_JFNK(Mat, Vec, Vec)
Definition: PetscIJacobianIMEX.cpp:128
_EXPLICIT_
#define _EXPLICIT_
Definition: petscinterface_struct.h:27
PetscCleanup
int PetscCleanup(void *)
Definition: PetscCleanup.cpp:12
PetscGlobalDOF
int PetscGlobalDOF(void *)
Definition: PetscGlobalDOF.cpp:69
_ROM_MODE_INITIAL_GUESS_
#define _ROM_MODE_INITIAL_GUESS_
Definition: librom_interface.h:19
PetscJacobianFunction_Linear
PetscErrorCode PetscJacobianFunction_Linear(Mat, Vec, Vec)
Definition: PetscIJacobian.cpp:236
HyPar::PhysicsOutput
int(* PhysicsOutput)(void *, void *, double)
Definition: hypar.h:347
PETScContext::points
std::vector< int * > points
Definition: petscinterface_struct.h:89
PetscIJacobian
PetscErrorCode PetscIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *)
Definition: PetscIJacobian.cpp:52
librom_interface.h
libROM interface class
OutputROMSolution
int OutputROMSolution(void *s, int nsims, double a_time)
Definition: OutputROMSolution.cpp:24
HyPar::dt
double dt
Definition: hypar.h:67
ResetFilenameIndex
void ResetFilenameIndex(char *f, int len)
Definition: IncrementFilename.c:64
PETScContext::stage_times
std::vector< double > stage_times
Definition: petscinterface_struct.h:69
PETScContext::op_times_arr
std::vector< double > op_times_arr
Definition: petscinterface_struct.h:141
WriteArray
int WriteArray(int, int, int *, int *, int, double *, double *, void *, void *, char *)
Definition: WriteArray.c:27
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
_IMPLICIT_
#define _IMPLICIT_
Definition: petscinterface_struct.h:29
PETScContext::precon_matrix_type
std::string precon_matrix_type
Definition: petscinterface_struct.h:110
SolvePETSc
int SolvePETSc(void *, int, int, int)
Integrate in time with PETSc.
Definition: SolvePETSc.cpp:50
libROMInterface
Class implementing interface with libROM.
Definition: librom_interface.h:45
CalculateROMDiff
int CalculateROMDiff(void *s, void *m)
Definition: CalculateROMDiff.c:24
TransferVecToPETSc
int TransferVecToPETSc(const double *const, Vec, void *, const int, const int)
Definition: TransferToPETSc.cpp:24
PETScContext::globalDOF
std::vector< double * > globalDOF
Definition: petscinterface_struct.h:124
simulation_object.h
Simulation object.
OutputSolution
int OutputSolution(void *, int, double)
Definition: OutputSolution.cpp:27
PETScContext
Structure containing the variables for time-integration with PETSc.
Definition: petscinterface_struct.h:42
_ROM_MODE_PREDICT_
#define _ROM_MODE_PREDICT_
Definition: librom_interface.h:16
PetscCreatePointList
int PetscCreatePointList(void *)
Definition: PetscCreatePointList.cpp:25
HyPar::n_iter
int n_iter
Definition: hypar.h:55
PetscPreTimeStep
PetscErrorCode PetscPreTimeStep(TS)
Definition: PetscPreTimeStep.cpp:24
PetscRHSFunctionIMEX
PetscErrorCode PetscRHSFunctionIMEX(TS, PetscReal, Vec, Vec, void *)
Definition: PetscRHSFunctionIMEX.cpp:49
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
basic.h
Some basic definitions and macros.
PetscRHSFunctionExpl
PetscErrorCode PetscRHSFunctionExpl(TS, PetscReal, Vec, Vec, void *)
Definition: PetscRHSFunctionExpl.cpp:36
PetscPreStage
PetscErrorCode PetscPreStage(TS, PetscReal)
Definition: PetscPreStage.cpp:14
PetscPostTimeStep
PetscErrorCode PetscPostTimeStep(TS)
Definition: PetscPostTimeStep.cpp:21
PetscJacobianFunctionIMEX_Linear
PetscErrorCode PetscJacobianFunctionIMEX_Linear(Mat, Vec, Vec)
Definition: PetscIJacobianIMEX.cpp:290
PetscSetInitialGuessROM
PetscErrorCode PetscSetInitialGuessROM(SNES, Vec, void *)
Definition: PetscSetInitialGuessROM.cpp:18
_ROM_MODE_NONE_
#define _ROM_MODE_NONE_
Definition: librom_interface.h:12
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
PetscIFunctionImpl
PetscErrorCode PetscIFunctionImpl(TS, PetscReal, Vec, Vec, Vec, void *)
Definition: PetscIFunctionImpl.cpp:37
_ROM_MODE_TRAIN_
#define _ROM_MODE_TRAIN_
Definition: librom_interface.h:14
HyPar::rom_diff_norms
double rom_diff_norms[3]
Definition: hypar.h:405
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
PetscPostStage
PetscErrorCode PetscPostStage(TS, PetscReal, PetscInt, Vec *)
Definition: PetscPostStage.cpp:16
PetscIFunctionIMEX
PetscErrorCode PetscIFunctionIMEX(TS, PetscReal, Vec, Vec, Vec, void *)
Definition: PetscIFunctionIMEX.cpp:53
HyPar::restart_iter
int restart_iter
Definition: hypar.h:58
PetscRegisterTIMethods
int PetscRegisterTIMethods(int)
Definition: PetscRegisterTIMethods.cpp:46