HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
Macros | Functions
PetscRegisterTIMethods.cpp File Reference

Register a custom time-integration method. More...

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <petscinterface_struct.h>

Go to the source code of this file.

Macros

#define __FUNCT__   "PetscRegisterTIMethods"
 

Functions

int PetscRegisterTIMethods (int rank)
 

Detailed Description

Register a custom time-integration method.

Author
Debojyoti Ghosh

Definition in file PetscRegisterTIMethods.cpp.

Macro Definition Documentation

◆ __FUNCT__

#define __FUNCT__   "PetscRegisterTIMethods"

Definition at line 14 of file PetscRegisterTIMethods.cpp.

Function Documentation

◆ PetscRegisterTIMethods()

int PetscRegisterTIMethods ( int  rank)

This function allows the registering of custom multi-stage time integration methods and using it with PETSc. More than one method may be provided, each with a unique name (the name must not conflict with the names of existing methods in PETSc). The methods should be provided in a file "time_method.inp", and the following must be provided:

  • Name
  • Class (i.e. RK or ARK)
  • Number of stages
  • Theoretical order
  • Butcher tableaux

The method can then be used by invoking it by its name. For example, if a custom ARKIMEX method "foo" is provided through "time_method.inp", it can be used by specifying "-ts_type arkimex -ts_arkimex_type foo" in the command line or ".petscrc" file.

See /Examples/PETScInputs for examples of the input files required to provide a time integration method.

Currently supported classes of time integrators for which custom methods can be provided:

To do:

  • Add support for TSGLEE when it gets merged to PETSc's master.
Parameters
rankMPI rank

Definition at line 46 of file PetscRegisterTIMethods.cpp.

47 {
48  PetscErrorCode ierr;
49  int ierr2;
50 
51  PetscFunctionBegin;
52 
53  /* Note: all processors read and register the custom methods */
54  /* instead of root doing it and broadcasting. */
55  FILE *in;
56  in = fopen("time_method.inp","r");
57  if (in) {
58  int n, N; /* N = total number of methods specified in the file */
59  ierr2 = fscanf(in,"%d",&N); if (ierr2 != 1) return(1);
60  for (n = 0; n < N; n++) {
61  char name[_MAX_STRING_SIZE_]; /* name of the scheme */
62  char type[_MAX_STRING_SIZE_]; /* type of scheme - ARKIMEX or RK */
63  PetscInt s, order; /* number of stages and order */
64  PetscReal *A , *b, *c; /* Butcher tableaux entries for non-stiff terms */
65  PetscReal *At, *bt, *ct; /* Butcher tableaux entries for the stiff terms */
66  PetscReal *bemb, *bembt; /* Embedded method coefficients */
67  PetscInt pinterp; /* order of interpolation scheme */
68  PetscReal *bint,*bintt; /* Dense output interpolation coefficients */
69 
70  /* Initializations */
71  strcpy(name,"");
72  s = order = pinterp = 0;
73  A = b = c = NULL;
74  At = bt = ct = NULL;
75  bemb = bembt = NULL;
76  bint = bintt = NULL;
77 
78  /* Read the method */
79  char word[_MAX_STRING_SIZE_];
80  ierr2 = fscanf(in,"%s",word); if (ierr2 != 1) return(1);
81  if (!strcmp(word,"begin")) {
82  while (strcmp(word, "end")) {
83  ierr2 = fscanf(in,"%s",word); if (ierr2 != 1) return(1);
84  if (!strcmp(word,"name")) { ierr2 = fscanf(in,"%s",name); if (ierr2 != 1) return(1); }
85  else if (!strcmp(word,"class")) { ierr2 = fscanf(in,"%s",type); if (ierr2 != 1) return(1); }
86  else if (!strcmp(word,"nstages")) { ierr2 = fscanf(in,"%d",&s); if (ierr2 != 1) return(1); }
87  else if (!strcmp(word,"order")) { ierr2 = fscanf(in,"%d",&order); if (ierr2 != 1) return(1); }
88  else if (!strcmp(word,"pinterp")) { ierr2 = fscanf(in,"%d",&pinterp);if (ierr2 != 1) return(1); }
89  else if (!strcmp(word,"At")) {
90  if (s == 0) {
91  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
92  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
93  return(1);
94  } else {
95  At = (PetscReal*) calloc (s*s, sizeof(PetscReal));
96  int i, j;
97  for (i = 0; i < s; i++) {
98  for (j = 0; j < s; j++) {
99  ierr2 = fscanf(in,"%lf",&At[i*s+j]); if (ierr2 != 1) return(1);
100  }
101  }
102  }
103  } else if (!strcmp(word,"A")) {
104  if (s == 0) {
105  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
106  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
107  return(1);
108  } else {
109  A = (PetscReal*) calloc (s*s, sizeof(PetscReal));
110  int i, j;
111  for (i = 0; i < s; i++) {
112  for (j = 0; j < s; j++) {
113  ierr2 = fscanf(in,"%lf",&A[i*s+j]); if (ierr2 != 1) return(1);
114  }
115  }
116  }
117  } else if (!strcmp(word,"bt")) {
118  if (s == 0) {
119  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
120  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
121  return(1);
122  } else {
123  bt = (PetscReal*) calloc (s, sizeof(PetscReal));
124  int i;
125  for (i = 0; i < s; i++) ierr2 = fscanf(in,"%lf",&bt[i]); if (ierr2 != 1) return(1);
126  }
127  } else if (!strcmp(word,"b")) {
128  if (s == 0) {
129  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
130  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
131  return(1);
132  } else {
133  b = (PetscReal*) calloc (s, sizeof(PetscReal));
134  int i;
135  for (i = 0; i < s; i++) ierr2 = fscanf(in,"%lf",&b[i]); if (ierr2 != 1) return(1);
136  }
137  } else if (!strcmp(word,"ct")) {
138  if (s == 0) {
139  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
140  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
141  return(1);
142  } else {
143  ct = (PetscReal*) calloc (s, sizeof(PetscReal));
144  int i;
145  for (i = 0; i < s; i++) ierr2 = fscanf(in,"%lf",&ct[i]); if (ierr2 != 1) return(1);
146  }
147  } else if (!strcmp(word,"c")) {
148  if (s == 0) {
149  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
150  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
151  return(1);
152  } else {
153  c = (PetscReal*) calloc (s, sizeof(PetscReal));
154  int i;
155  for (i = 0; i < s; i++) ierr2 = fscanf(in,"%lf",&c[i]); if (ierr2 != 1) return(1);
156  }
157  } else if (!strcmp(word,"bembt")) {
158  if (s == 0) {
159  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
160  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
161  return(1);
162  } else {
163  bembt = (PetscReal*) calloc (s, sizeof(PetscReal));
164  int i;
165  for (i = 0; i < s; i++) ierr2 = fscanf(in,"%lf",&bembt[i]); if (ierr2 != 1) return(1);
166  }
167  } else if (!strcmp(word,"bemb")) {
168  if (s == 0) {
169  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
170  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
171  return(1);
172  } else {
173  bemb = (PetscReal*) calloc (s, sizeof(PetscReal));
174  int i;
175  for (i = 0; i < s; i++) ierr2 = fscanf(in,"%lf",&bemb[i]); if (ierr2 != 1) return(1);
176  }
177  } else if (!strcmp(word,"bintt")) {
178  if (s == 0 || pinterp == 0) {
179  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages and pinterp must be " );
180  if (!rank) fprintf(stderr,"defined as positive values before specifying interpolation coeffs.\n");
181  return(1);
182  } else {
183  bintt = (PetscReal*) calloc (s*pinterp, sizeof(PetscReal));
184  int i, j;
185  for (i = 0; i < s; i++) {
186  for (j = 0; j < pinterp; j++) {
187  ierr2 = fscanf(in,"%lf",&bintt[i*s+j]); if (ierr2 != 1) return(1);
188  }
189  }
190  }
191  } else if (!strcmp(word,"bint")) {
192  if (s == 0 || pinterp == 0) {
193  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages and pinterp must be " );
194  if (!rank) fprintf(stderr,"defined as positive values before specifying interpolation coeffs.\n");
195  return(1);
196  } else {
197  bint = (PetscReal*) calloc (s*pinterp, sizeof(PetscReal));
198  int i, j;
199  for (i = 0; i < s; i++) {
200  for (j = 0; j < pinterp; j++) {
201  ierr2 = fscanf(in,"%lf",&bint[i*s+j]); if (ierr2 != 1) return(1);
202  }
203  }
204  }
205  }
206  }
207  } else {
208  if (!rank) fprintf(stderr,"Error: Illegal format in file \"time_method.inp\" (expected keyword \"begin\").\n");
209  return(1);
210  }
211 
212  /* Register the method */
213  if (!strcmp(type,"arkimex")) {
214  if (A && At) {
215  ierr = TSARKIMEXRegister(name,order,s,At,bt,ct,A,b,c,bembt,bemb,pinterp,bintt,bint); CHKERRQ(ierr);
216  if (!rank) {
217  printf("\nRegistered custom ARKIMEX scheme \"%s\" with the following Butcher tableaux:-\n",name);
218  int i,j;
219  for (i = 0; i < s; i++) {
220  if (c) printf(" %+1.5lf |",c[i]);
221  else printf(" |");
222  for (j = 0; j < s; j++) printf (" %+1.5lf :",A[i*s+j]);
223  printf(" ");
224  if (ct) printf("%+1.5lf |",ct[i]);
225  else printf(" |");
226  for (j = 0; j < s; j++) printf (" %+1.5lf :",At[i*s+j]);
227  printf("\n");
228  }
229  printf(" ---------|");
230  for (j = 0; j < s; j++) printf("-----------");
231  printf(" ");
232  printf("---------|");
233  for (j = 0; j < s; j++) printf("-----------");
234  printf("\n");
235  printf(" |");
236  if (b) for (j = 0; j < s; j++) printf(" %+1.5lf :",b[j]);
237  else for (j = 0; j < s; j++) printf(" :");
238  printf(" ");
239  printf(" |");
240  if (bt) for (j = 0; j < s; j++) printf(" %+1.5lf :",bt[j]);
241  else for (j = 0; j < s; j++) printf(" :");
242  printf("\n\n");
243  }
244  } else {
245  if (!rank) {
246  fprintf(stderr,"Warning in PetscRegisterTIMethods(): Failed to register method ");
247  fprintf(stderr,"(A or At not defined).\n");
248  }
249  }
250  } else if (!strcmp(type,"rk")) {
251  if (A) {
252  ierr = TSRKRegister(name,order,s,A,b,c,bemb,pinterp,bint); CHKERRQ(ierr);
253  if (!rank) {
254  printf("\nRegistered custom RK scheme \"%s\" with the following Butcher tableaux:-\n",name);
255  int i,j;
256  for (i = 0; i < s; i++) {
257  if (c) printf(" %+1.5lf |",c[i]);
258  else printf(" |");
259  for (j = 0; j < s; j++) printf (" %+1.5lf :",A[i*s+j]);
260  printf("\n");
261  }
262  printf(" ---------|");
263  for (j = 0; j < s; j++) printf("-----------");
264  printf("\n");
265  printf(" |");
266  if (b) for (j = 0; j < s; j++) printf(" %+1.5lf :",b[j]);
267  else for (j = 0; j < s; j++) printf(" :");
268  printf("\n\n");
269  }
270  } else {
271  if (!rank) {
272  fprintf(stderr,"Warning in PetscRegisterTIMethods(): Failed to register method ");
273  fprintf(stderr,"(A not defined).\n");
274  }
275  }
276  } else {
277  if (!rank){
278  fprintf(stderr,"Error in PetscRegisterTIMethods(): %s class of time-integration schemes ",type);
279  fprintf(stderr,"does not support custom method registration and usage.\n");
280  }
281  }
282 
283  /* Free the arrays */
284  if (At) free(At);
285  if (bt) free(bt);
286  if (ct) free(ct);
287  if (A) free(A);
288  if (b) free(b);
289  if (c) free(c);
290  if (bembt) free(bembt);
291  if (bemb) free(bemb);
292  if (bintt) free(bintt);
293  if (bint) free(bint);
294 
295  }
296  }
297  PetscFunctionReturn(0);
298 }
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14