HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
Solve.cpp File Reference
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <vector>
#include <common_cpp.h>
#include <io_cpp.h>
#include <timeintegration_cpp.h>
#include <mpivars_cpp.h>
#include <simulation_object.h>
#include <librom_interface.h>

Go to the source code of this file.

Functions

int ComputeRHSOperators (void *, void *, double)
 
int CalculateError (void *, void *)
 
int OutputSolution (void *, int, double)
 
void ResetFilenameIndex (char *, int)
 
int CalculateROMDiff (void *, void *)
 
int OutputROMSolution (void *, int, double)
 
int Solve (void *s, int nsims, int rank, int nproc)
 

Function Documentation

int ComputeRHSOperators ( void *  s,
void *  m,
double  t 
)

Computes the matrices representing the Jacobians of the hyperbolic (HyPar::HyperbolicFunction), parabolic (HyPar::ParabolicFunction), and the source (HyPar::SourceFunction) terms.

Each element is compute through finite-differences, by perturbing one DOF of the solution at a time. The matrices are not stored in the memory, instead the non-zero are written to file as soon as they are computed. The filenames for the matrices are "Mat_*Function_nnnnn.dat", where "nnnnn" is a time-dependent index.

  • The format (ASCII text) of the file is as follows:
    n
    i j val
    i j val
    ...
    i j val
    where n is the size of the matrix, followed by all the non-zero elements (each line contains the indices i,j and the value val of one non-zero element).
  • This function is called after HyPar::file_op_iter iterations (i.e. the same frequency at which solution files are written).
  • If a splitting for the hyperbolic flux is defined, then the Jacobians of the complete hyperbolic term, as well as the split terms are computed.
  • The eigenvalues of these matrices can be computed and plotted in MATLAB using the scripts Examples/Matlab/ComputeEvals.m and Examples/Matlab/PlotEvals.m respectively.
  • To use this function, the code must be compiled with the flag -Dcompute_rhs_operators, and run on a single processor. This is very slow for larger domains, and should be used only for the numerical analysis of test cases.
  • The current solution stored in HyPar::u is used as the reference state at which the Jacobians are computed (this is relevant to know for non-linear right-hand-sides functions).

Note: Parabolic terms have not yet been included.

Parameters
sSolver object of type HyPar
mMPI object of type MPIVariables
tCurrent simulation time

Definition at line 50 of file ComputeRHSOperators.c.

55 {
56  HyPar *solver = (HyPar*) s;
57  MPIVariables *mpi = (MPIVariables*) m;
58  int i, j, size, ndof;
59  double *f, *f0, *u, *u0, *rhs, *drhs;
60  FILE *fout;
61  char filename[_MAX_STRING_SIZE_];
62 
63  if (mpi->nproc > 1) return(0);
64 
65  int ndims = solver->ndims;
66  int nvars = solver->nvars;
67  int ghosts = solver->ghosts;
68  int *dim = solver->dim_local;
69  int index[ndims];
70 
71  double epsilon = 1e-6;
72  double tolerance = 1e-15;
73 
74  /* allocate arrays */
75  size = solver->npoints_local_wghosts * nvars;
76  u = (double*) calloc (size,sizeof(double));
77  u0 = (double*) calloc (size,sizeof(double));
78  rhs = (double*) calloc (size,sizeof(double));
79  drhs = (double*) calloc (size,sizeof(double));
80  ndof = solver->npoints_local * nvars;
81  f = (double*) calloc (ndof,sizeof(double));
82  f0 = (double*) calloc (ndof,sizeof(double));
83 
84  /* copy the current solution to u0 */
85  _ArrayCopy1D_(solver->u,u0,size);
86  /* apply boundary conditions to the solution u0 */
87  IERR solver->ApplyBoundaryConditions(solver,mpi,u0,NULL,t);CHECKERR(ierr);
88  IERR solver->ApplyIBConditions(solver,mpi,u0,t);CHECKERR(ierr);
89  IERR MPIExchangeBoundariesnD(ndims,nvars,dim,ghosts,mpi,u0); CHECKERR(ierr);
90 
91  /* compute linearized matrix for the hyperbolic FFunction */
92  if (solver->FFunction) {
93  strcpy(filename,"Mat_FFunction_");
94  strcat(filename,solver->filename_index);
95  strcat(filename,".dat");
96  printf("ComputeRHSOperators(): Computing linearized matrix operator for FFunction. ndof=%d.\n",ndof);
97  printf("ComputeRHSOperators(): Writing to sparse matrix file %s.\n",filename);
98  fout = fopen(filename,"w");
99  fprintf(fout,"%d\n",ndof);
100  /* compute the FFunction of u0 */
101  _ArraySetValue_(f0,ndof,0.0);
102  IERR solver->HyperbolicFunction(rhs,u0,solver,mpi,t,1,solver->FFunction,
103  solver->Upwind); CHECKERR(ierr);
104  _ArrayScale1D_(rhs,-1.0,size);
105  /* transfer to an array f0 which as no ghost points */
106  IERR ArrayCopynD(ndims,rhs,f0,dim,ghosts,0,index,nvars); CHECKERR(ierr);
107  for (i=0; i<ndof; i++) {
108  /* copy the solution u0 to u */
109  _ArrayCopy1D_(u0,u,size);
110  /* find the 1D index p in an array with ghosts points (u),
111  * corresponding to index i which assumes no ghost points */
112  int ii, p, v;
113  ii = i / nvars;
114  v = i - ii*nvars;
115  _ArrayIndexnD_(ndims,ii,dim,index,0);
116  _ArrayIndex1D_(ndims,dim,index,ghosts,p);
117  /* add a perturbation */
118  u[nvars*p+v] += epsilon;
119  /* apply boundary conditions to the perturbed solution u */
120  IERR solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);CHECKERR(ierr);
121  IERR solver->ApplyIBConditions(solver,mpi,u,t);CHECKERR(ierr);
122  IERR MPIExchangeBoundariesnD(ndims,nvars,dim,ghosts,mpi,u); CHECKERR(ierr);
123  /* compute the FFunction of u */
124  IERR solver->HyperbolicFunction(rhs,u,solver,mpi,t,1,solver->FFunction,
125  solver->Upwind); CHECKERR(ierr);
126  _ArrayScale1D_(rhs,-1.0,size);
127  /* transfer to an array f which as no ghost points */
128  _ArraySetValue_(f,ndof,0.0);
129  IERR ArrayCopynD(ndims,rhs,f,dim,ghosts,0,index,nvars); CHECKERR(ierr);
130  /* subtract f0 from f */
131  _ArrayAXPY_(f0,-1.0,f,ndof);
132  /* f is now espilon*column of the matrix */
133  for (j=0; j<ndof; j++) {
134  double mat_elem = f[j] / epsilon;
135  /* write to file if element is non-zero */
136  if (absolute(mat_elem) > tolerance) fprintf(fout,"%5d %5d %+1.16e\n",j+1,i+1,mat_elem);
137  }
138  }
139  fclose(fout);
140  }
141  /* compute linearized matrix for the split hyperbolic (FFunction-dFFunction) and dFFunction
142  functions, if dFFunction is available */
143  if (solver->dFFunction) {
144  strcpy(filename,"Mat_FdFFunction_");
145  strcat(filename,solver->filename_index);
146  strcat(filename,".dat");
147  printf("ComputeRHSOperators(): Computing linearized matrix operator for (FFunction-dFFunction). ndof=%d.\n",ndof);
148  printf("ComputeRHSOperators(): Writing to sparse matrix file %s.\n",filename);
149  fout = fopen(filename,"w");
150  fprintf(fout,"%d\n",ndof);
151  /* compute the FFunction of u0 */
152  _ArraySetValue_(f0,ndof,0.0);
153  if (solver->flag_fdf_specified) {
154  IERR solver->HyperbolicFunction(rhs,u0,solver,mpi,t,1,solver->FdFFunction,
155  solver->UpwindFdF); CHECKERR(ierr);
156  _ArrayScale1D_(rhs,-1.0,size);
157  } else {
158  IERR solver->HyperbolicFunction(rhs,u0,solver,mpi,t,1,solver->FFunction,
159  solver->Upwind); CHECKERR(ierr);
160  IERR solver->HyperbolicFunction(drhs,u0,solver,mpi,t,0,solver->dFFunction,
161  solver->UpwinddF); CHECKERR(ierr);
162  _ArrayScale1D_(rhs,-1.0,size);
163  _ArrayAXPY_(drhs,1.0,rhs,size);
164  }
165  /* transfer to an array f0 which as no ghost points */
166  IERR ArrayCopynD(ndims,rhs,f0,dim,ghosts,0,index,nvars); CHECKERR(ierr);
167  for (i=0; i<ndof; i++) {
168  /* copy the solution u0 to u */
169  _ArrayCopy1D_(u0,u,size);
170  /* find the 1D index p in an array with ghosts points (u),
171  * corresponding to index i which assumes no ghost points */
172  int ii, p, v;
173  ii = i / nvars;
174  v = i - ii*nvars;
175  _ArrayIndexnD_(ndims,ii,dim,index,0);
176  _ArrayIndex1D_(ndims,dim,index,ghosts,p);
177  /* add a perturbation */
178  u[nvars*p+v] += epsilon;
179  /* apply boundary conditions to the perturbed solution u */
180  IERR solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);CHECKERR(ierr);
181  IERR solver->ApplyIBConditions(solver,mpi,u,t);CHECKERR(ierr);
182  IERR MPIExchangeBoundariesnD(ndims,nvars,dim,ghosts,mpi,u); CHECKERR(ierr);
183  /* compute the FFunction of u */
184  if (solver->flag_fdf_specified) {
185  IERR solver->HyperbolicFunction(rhs,u,solver,mpi,t,1,solver->FdFFunction,
186  solver->UpwindFdF); CHECKERR(ierr);
187  _ArrayScale1D_(rhs,-1.0,size);
188  } else {
189  IERR solver->HyperbolicFunction(rhs,u,solver,mpi,t,1,solver->FFunction,
190  solver->Upwind); CHECKERR(ierr);
191  IERR solver->HyperbolicFunction(drhs,u,solver,mpi,t,0,solver->dFFunction,
192  solver->UpwinddF); CHECKERR(ierr);
193  _ArrayScale1D_(rhs,-1.0,size);
194  _ArrayAXPY_(drhs,1.0,rhs,size);
195  }
196  /* transfer to an array f which as no ghost points */
197  _ArraySetValue_(f,ndof,0.0);
198  IERR ArrayCopynD(ndims,rhs,f,dim,ghosts,0,index,nvars); CHECKERR(ierr);
199  /* subtract f0 from f */
200  _ArrayAXPY_(f0,-1.0,f,ndof);
201  /* f is now espilon*column of the matrix */
202  for (j=0; j<ndof; j++) {
203  double mat_elem = f[j] / epsilon;
204  /* write to file if element is non-zero */
205  if (absolute(mat_elem) > tolerance) fprintf(fout,"%5d %5d %+1.16e\n",j+1,i+1,mat_elem);
206  }
207  }
208  fclose(fout);
209 
210  strcpy(filename,"Mat_dFFunction_");
211  strcat(filename,solver->filename_index);
212  strcat(filename,".dat");
213  printf("ComputeRHSOperators(): Computing linearized matrix operator for dFFunction. ndof=%d.\n",ndof);
214  printf("ComputeRHSOperators(): Writing to sparse matrix file %s.\n",filename);
215  fout = fopen(filename,"w");
216  fprintf(fout,"%d\n",ndof);
217  /* compute the FFunction of u0 */
218  _ArraySetValue_(f0,ndof,0.0);
219  IERR solver->HyperbolicFunction(rhs,u0,solver,mpi,t,0,solver->dFFunction,
220  solver->UpwinddF); CHECKERR(ierr);
221  _ArrayScale1D_(rhs,-1.0,size);
222  /* transfer to an array f0 which as no ghost points */
223  IERR ArrayCopynD(ndims,rhs,f0,dim,ghosts,0,index,nvars); CHECKERR(ierr);
224  for (i=0; i<ndof; i++) {
225  /* copy the solution u0 to u */
226  _ArrayCopy1D_(u0,u,size);
227  /* find the 1D index p in an array with ghosts points (u),
228  * corresponding to index i which assumes no ghost points */
229  int ii, p, v;
230  ii = i / nvars;
231  v = i - ii*nvars;
232  _ArrayIndexnD_(ndims,ii,dim,index,0);
233  _ArrayIndex1D_(ndims,dim,index,ghosts,p);
234  /* add a perturbation */
235  u[nvars*p+v] += epsilon;
236  /* apply boundary conditions to the perturbed solution u */
237  IERR solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);CHECKERR(ierr);
238  IERR solver->ApplyIBConditions(solver,mpi,u,t);CHECKERR(ierr);
239  IERR MPIExchangeBoundariesnD(ndims,nvars,dim,ghosts,mpi,u); CHECKERR(ierr);
240  /* compute the FFunction of u */
241  IERR solver->HyperbolicFunction(rhs,u,solver,mpi,t,0,solver->dFFunction,
242  solver->UpwinddF); CHECKERR(ierr);
243  _ArrayScale1D_(rhs,-1.0,size);
244  /* transfer to an array f which as no ghost points */
245  _ArraySetValue_(f,ndof,0.0);
246  IERR ArrayCopynD(ndims,rhs,f,dim,ghosts,0,index,nvars); CHECKERR(ierr);
247  /* subtract f0 from f */
248  _ArrayAXPY_(f0,-1.0,f,ndof);
249  /* f is now espilon*column of the matrix */
250  for (j=0; j<ndof; j++) {
251  double mat_elem = f[j] / epsilon;
252  /* write to file if element is non-zero */
253  if (absolute(mat_elem) > tolerance) fprintf(fout,"%5d %5d %+1.16e\n",j+1,i+1,mat_elem);
254  }
255  }
256  fclose(fout);
257  }
258 
259  /* compute linearized matrix for the source SFunction */
260  if (solver->SFunction) {
261  strcpy(filename,"Mat_SFunction_");
262  strcat(filename,solver->filename_index);
263  strcat(filename,".dat");
264  printf("ComputeRHSOperators(): Computing linearized matrix operator for SFunction. ndof=%d.\n",ndof);
265  printf("ComputeRHSOperators(): Writing to sparse matrix file %s.\n",filename);
266  fout = fopen(filename,"w");
267  fprintf(fout,"%d\n",ndof);
268  /* compute the FFunction of u0 */
269  _ArraySetValue_(f0,ndof,0.0);
270  IERR solver->SourceFunction(rhs,u0,solver,mpi,t); CHECKERR(ierr);
271  /* transfer to an array f0 which as no ghost points */
272  IERR ArrayCopynD(ndims,rhs,f0,dim,ghosts,0,index,nvars); CHECKERR(ierr);
273  for (i=0; i<ndof; i++) {
274  /* copy the solution u0 to u */
275  _ArrayCopy1D_(u0,u,size);
276  /* find the 1D index p in an array with ghosts points (u),
277  * corresponding to index i which assumes no ghost points */
278  int ii, p, v;
279  ii = i / nvars;
280  v = i - ii*nvars;
281  _ArrayIndexnD_(ndims,ii,dim,index,0);
282  _ArrayIndex1D_(ndims,dim,index,ghosts,p);
283  /* add a perturbation */
284  u[nvars*p+v] += epsilon;
285  /* apply boundary conditions to the perturbed solution u */
286  IERR solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);CHECKERR(ierr);
287  IERR solver->ApplyIBConditions(solver,mpi,u,t);CHECKERR(ierr);
288  IERR MPIExchangeBoundariesnD(ndims,nvars,dim,ghosts,mpi,u); CHECKERR(ierr);
289  /* compute the SFunction of u */
290  IERR solver->SourceFunction(rhs,u,solver,mpi,t); CHECKERR(ierr);
291  /* transfer to an array f which as no ghost points */
292  _ArraySetValue_(f,ndof,0.0);
293  IERR ArrayCopynD(ndims,rhs,f,dim,ghosts,0,index,nvars); CHECKERR(ierr);
294  /* subtract f0 from f */
295  _ArrayAXPY_(f0,-1.0,f,ndof);
296  /* f is now espilon*column of the matrix */
297  for (j=0; j<ndof; j++) {
298  double mat_elem = f[j] / epsilon;
299  /* write to file if element is non-zero */
300  if (absolute(mat_elem) > tolerance) fprintf(fout,"%5d %5d %+1.16e\n",j+1,i+1,mat_elem);
301  }
302  }
303  fclose(fout);
304  }
305 
306  /* clean up */
307  free(u);
308  free(u0);
309  free(rhs);
310  free(drhs);
311  free(f);
312  free(f0);
313  return(0);
314 }
HyPar::npoints_local_wghosts
int npoints_local_wghosts
Definition: hypar.h:42
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::dFFunction
int(* dFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:280
HyPar::npoints_local
int npoints_local
Definition: hypar.h:42
HyPar::u
double * u
Definition: hypar.h:116
_ArrayIndexnD_
#define _ArrayIndexnD_(N, index, imax, i, ghost)
Definition: arrayfunctions.h:21
ArrayCopynD
INLINE int ArrayCopynD(int, const double *, double *, int *, int, int, int *, int)
Definition: arrayfunctions.h:410
HyPar::Upwind
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
HyPar::flag_fdf_specified
int flag_fdf_specified
Definition: hypar.h:290
HyPar::UpwinddF
int(* UpwinddF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:300
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
HyPar::UpwindFdF
int(* UpwindFdF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:307
HyPar::FdFFunction
int(* FdFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:286
HyPar::FFunction
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
_ArrayScale1D_
#define _ArrayScale1D_(x, a, size)
Definition: arrayfunctions.h:190
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
HyPar::filename_index
char * filename_index
Definition: hypar.h:197
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
HyPar::SourceFunction
int(* SourceFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:259
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::HyperbolicFunction
int(* HyperbolicFunction)(double *, double *, void *, void *, double, int, int(*)(double *, double *, int, void *, double), int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
Definition: hypar.h:250
HyPar::SFunction
int(* SFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:317
HyPar::nvars
int nvars
Definition: hypar.h:29
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
HyPar::ndims
int ndims
Definition: hypar.h:26
absolute
#define absolute(a)
Definition: math_ops.h:32
HyPar::ApplyBoundaryConditions
int(* ApplyBoundaryConditions)(void *, void *, double *, double *, double)
Definition: hypar.h:214
IERR
#define IERR
Definition: basic.h:16
MPIExchangeBoundariesnD
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Definition: MPIExchangeBoundariesnD.c:42
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
_ArrayAXPY_
#define _ArrayAXPY_(x, a, y, size)
Definition: arrayfunctions.h:271
HyPar::ApplyIBConditions
int(* ApplyIBConditions)(void *, void *, double *, double)
Definition: hypar.h:217
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int CalculateError ( void *  s,
void *  m 
)

Calculate the error in the final solution

Calculates the error in the solution if the exact solution is available. If the exact solution is not available, the errors are reported as zero. The exact solution should be provided in the file "exact.inp" in the same format as the initial solution.

Parameters
sSolver object of type HyPar
mMPI object of type MPIVariables

Definition at line 25 of file CalculateError.c.

29 {
30  HyPar *solver = (HyPar*) s;
31  MPIVariables *mpi = (MPIVariables*) m;
32  int exact_flag = 0, i, size;
33  double sum = 0, global_sum = 0;
34  double *uex = NULL;
35  _DECLARE_IERR_;
36 
37  size = solver->nvars;
38  for (i = 0; i < solver->ndims; i++)
39  size *= (solver->dim_local[i]+2*solver->ghosts);
40  uex = (double*) calloc (size, sizeof(double));
41 
42  char fname_root[_MAX_STRING_SIZE_] = "exact";
43  if (solver->nsims > 1) {
44  char index[_MAX_STRING_SIZE_];
45  GetStringFromInteger(solver->my_idx, index, (int)log10(solver->nsims)+1);
46  strcat(fname_root, "_");
47  strcat(fname_root, index);
48  }
49 
50  static const double tolerance = 1e-15;
51  IERR ExactSolution( solver,
52  mpi,
53  uex,
54  fname_root,
55  &exact_flag ); CHECKERR(ierr);
56 
57  if (!exact_flag) {
58 
59  /* No exact solution */
60  IERR TimeError(solver,mpi,NULL); CHECKERR(ierr);
61  solver->error[0] = solver->error[1] = solver->error[2] = -1;
62 
63  } else {
64 
65  IERR TimeError(solver,mpi,uex); CHECKERR(ierr);
66 
67  /* calculate solution norms (for relative error) */
68  double solution_norm[3] = {0.0,0.0,0.0};
69  /* L1 */
70  sum = ArraySumAbsnD (solver->nvars,solver->ndims,solver->dim_local,
71  solver->ghosts,solver->index,uex);
72  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
73  solution_norm[0] = global_sum/((double)solver->npoints_global);
74  /* L2 */
75  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
76  solver->ghosts,solver->index,uex);
77  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
78  solution_norm[1] = sqrt(global_sum/((double)solver->npoints_global));
79  /* Linf */
80  sum = ArrayMaxnD (solver->nvars,solver->ndims,solver->dim_local,
81  solver->ghosts,solver->index,uex);
82  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
83  solution_norm[2] = global_sum;
84 
85  /* compute error = difference between exact and numerical solution */
86  _ArrayAXPY_(solver->u,-1.0,uex,size);
87 
88  /* calculate L1 norm of error */
89  sum = ArraySumAbsnD (solver->nvars,solver->ndims,solver->dim_local,
90  solver->ghosts,solver->index,uex);
91  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
92  solver->error[0] = global_sum/((double)solver->npoints_global);
93 
94  /* calculate L2 norm of error */
95  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
96  solver->ghosts,solver->index,uex);
97  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
98  solver->error[1] = sqrt(global_sum/((double)solver->npoints_global));
99 
100  /* calculate Linf norm of error */
101  sum = ArrayMaxnD (solver->nvars,solver->ndims,solver->dim_local,
102  solver->ghosts,solver->index,uex);
103  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
104  solver->error[2] = global_sum;
105 
106  /*
107  decide whether to normalize and report relative errors,
108  or report absolute errors.
109  */
110  if ( (solution_norm[0] > tolerance)
111  && (solution_norm[1] > tolerance)
112  && (solution_norm[2] > tolerance) ) {
113  solver->error[0] /= solution_norm[0];
114  solver->error[1] /= solution_norm[1];
115  solver->error[2] /= solution_norm[2];
116  }
117  }
118 
119  free(uex);
120  return(0);
121 }
HyPar::u
double * u
Definition: hypar.h:116
GetStringFromInteger
void GetStringFromInteger(int, char *, int)
Definition: CommonFunctions.c:15
TimeError
int TimeError(void *, void *, double *)
Definition: TimeError.c:27
HyPar::npoints_global
int npoints_global
Definition: hypar.h:40
ArraySumAbsnD
INLINE double ArraySumAbsnD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:510
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
MPISum_double
int MPISum_double(double *, double *, int, void *)
Definition: MPISum.c:39
MPIMax_double
int MPIMax_double(double *, double *, int, void *)
Definition: MPIMax.c:38
HyPar::ghosts
int ghosts
Definition: hypar.h:52
HyPar::nsims
int nsims
Definition: hypar.h:64
MPIVariables::world
MPI_Comm world
Definition: mpivars_struct.h:37
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
ArraySumSquarenD
INLINE double ArraySumSquarenD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:529
HyPar::nvars
int nvars
Definition: hypar.h:29
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
StdVecOps::sum
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18
ExactSolution
int ExactSolution(void *, void *, double *, char *, int *)
Definition: ExactSolution.c:16
HyPar::ndims
int ndims
Definition: hypar.h:26
HyPar::index
int * index
Definition: hypar.h:102
HyPar::my_idx
int my_idx
Definition: hypar.h:61
IERR
#define IERR
Definition: basic.h:16
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
ArrayMaxnD
INLINE double ArrayMaxnD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:486
_ArrayAXPY_
#define _ArrayAXPY_(x, a, y, size)
Definition: arrayfunctions.h:271
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
HyPar::error
double error[3]
Definition: hypar.h:371
int OutputSolution ( void *  s,
int  nsims,
double  a_time 
)

Write solutions to file

Write out the solution to file

Parameters
sArray of simulation objects of type SimulationObject
nsimsNumber of simulation objects
a_timeCurrent simulation time

Definition at line 27 of file OutputSolution.cpp.

30 {
31  SimulationObject* simobj = (SimulationObject*) s;
32  int ns;
33  _DECLARE_IERR_;
34 
35  for (ns = 0; ns < nsims; ns++) {
36 
37  HyPar* solver = &(simobj[ns].solver);
38  MPIVariables* mpi = &(simobj[ns].mpi);
39 
40  if ((!solver->WriteOutput) && (strcmp(solver->plot_solution,"yes"))) continue;
41 
42  /* time integration module may have auxiliary arrays to write out, so get them */
43  int NSolutions = 0;
44  IERR TimeGetAuxSolutions(&NSolutions,NULL,solver,-1,ns); CHECKERR(ierr);
45  if (NSolutions > 10) NSolutions = 10;
46 
47  int nu;
48  char fname_root[_MAX_STRING_SIZE_];
49  char aux_fname_root[_MAX_STRING_SIZE_];
50  strcpy(fname_root, solver->op_fname_root);
51  strcpy(aux_fname_root, solver->aux_op_fname_root);
52 
53  if (nsims > 1) {
54  char index[_MAX_STRING_SIZE_];
55  GetStringFromInteger(ns, index, (int)log10(nsims)+1);
56  strcat(fname_root, "_");
57  strcat(fname_root, index);
58  strcat(aux_fname_root, "_");
59  strcat(aux_fname_root, index);
60  }
61 
62  for (nu=0; nu<NSolutions; nu++) {
63 
64  double *uaux = NULL;
65  IERR TimeGetAuxSolutions(&NSolutions,&uaux,solver,nu,ns); CHECKERR(ierr);
66 
67  IERR WriteArray( solver->ndims,
68  solver->nvars,
69  solver->dim_global,
70  solver->dim_local,
71  solver->ghosts,
72  solver->x,
73  uaux,
74  solver,
75  mpi,
76  aux_fname_root ); CHECKERR(ierr);
77 
78  aux_fname_root[2]++;
79  }
80 
81 #if defined(HAVE_CUDA)
82  if (solver->use_gpu) {
83  /* Copy values from GPU memory to CPU memory for writing. */
84  gpuMemcpy(solver->x, solver->gpu_x, sizeof(double)*solver->size_x, gpuMemcpyDeviceToHost);
85 
86 #ifdef CUDA_VAR_ORDERDING_AOS
87  gpuMemcpy( solver->u,
88  solver->gpu_u,
89  sizeof(double)*solver->ndof_cells_wghosts,
90  gpuMemcpyDeviceToHost );
91 #else
92  double *h_u = (double *) malloc(sizeof(double)*solver->ndof_cells_wghosts);
93  gpuMemcpy(h_u, solver->gpu_u, sizeof(double)*solver->ndof_cells_wghosts, gpuMemcpyDeviceToHost);
94  for (int i=0; i<solver->npoints_local_wghosts; i++) {
95  for (int v=0; v<solver->nvars; v++) {
96  solver->u[i*solver->nvars+v] = h_u[i+v*solver->npoints_local_wghosts];
97  }
98  }
99  free(h_u);
100 #endif
101  }
102 #endif
103 
104  WriteArray( solver->ndims,
105  solver->nvars,
106  solver->dim_global,
107  solver->dim_local,
108  solver->ghosts,
109  solver->x,
110  solver->u,
111  solver,
112  mpi,
113  fname_root );
114 
115  if (!strcmp(solver->plot_solution, "yes")) {
116  PlotArray( solver->ndims,
117  solver->nvars,
118  solver->dim_global,
119  solver->dim_local,
120  solver->ghosts,
121  solver->x,
122  solver->u,
123  a_time,
124  solver,
125  mpi,
126  fname_root );
127  }
128 
129  /* increment the index string, if required */
130  if ((!strcmp(solver->output_mode,"serial")) && (!strcmp(solver->op_overwrite,"no"))) {
131  IncrementFilenameIndex(solver->filename_index,solver->index_length);
132  }
133 
134  }
135 
136  return(0);
137 }
HyPar::npoints_local_wghosts
int npoints_local_wghosts
Definition: hypar.h:42
HyPar::dim_global
int * dim_global
Definition: hypar.h:33
HyPar::u
double * u
Definition: hypar.h:116
HyPar::aux_op_fname_root
char aux_op_fname_root[_MAX_STRING_SIZE_]
Definition: hypar.h:208
GetStringFromInteger
void GetStringFromInteger(int, char *, int)
Definition: CommonFunctions.c:15
HyPar::plot_solution
char plot_solution[_MAX_STRING_SIZE_]
Definition: hypar.h:194
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
gpuMemcpy
void gpuMemcpy(void *, const void *, size_t, enum gpuMemcpyKind)
Definition: ArrayImplementations_GPU.cu:94
TimeGetAuxSolutions
int TimeGetAuxSolutions(int *, double **, void *, int, int)
Definition: TimeGetAuxSolutions.c:29
HyPar::ghosts
int ghosts
Definition: hypar.h:52
HyPar::filename_index
char * filename_index
Definition: hypar.h:197
HyPar::op_fname_root
char op_fname_root[_MAX_STRING_SIZE_]
Definition: hypar.h:206
HyPar::size_x
int size_x
Definition: hypar.h:48
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
HyPar::gpu_u
double * gpu_u
Definition: hypar.h:459
HyPar::op_overwrite
char op_overwrite[_MAX_STRING_SIZE_]
Definition: hypar.h:191
HyPar::ndof_cells_wghosts
int ndof_cells_wghosts
Definition: hypar.h:47
HyPar::gpu_x
double * gpu_x
Definition: hypar.h:457
HyPar::nvars
int nvars
Definition: hypar.h:29
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
PlotArray
int PlotArray(int, int, int *, int *, int, double *, double *, double, void *, void *, char *)
Definition: PlotArray.c:31
HyPar::use_gpu
int use_gpu
Definition: hypar.h:449
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
HyPar::ndims
int ndims
Definition: hypar.h:26
IERR
#define IERR
Definition: basic.h:16
HyPar::index_length
int index_length
Definition: hypar.h:199
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar::x
double * x
Definition: hypar.h:107
HyPar::output_mode
char output_mode[_MAX_STRING_SIZE_]
Definition: hypar.h:183
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
IncrementFilenameIndex
void IncrementFilenameIndex(char *f, int len)
Definition: IncrementFilename.c:46
HyPar::WriteOutput
int(* WriteOutput)(int, int, int *, double *, double *, char *, int *)
Definition: hypar.h:211
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
void ResetFilenameIndex ( char *  f,
int  len 
)

Reset filename index

Resets the index to "0000..." of a desired length.

Parameters
fCharacter string representing the integer
lenLength of the string

Definition at line 64 of file IncrementFilename.c.

66 {
67  if (!f) return;
68  int i;
69  for (i = 0; i < len; i++) {
70  f[i] = '0';
71  }
72  return;
73 }
int CalculateROMDiff ( void *  s,
void *  m 
)

Calculate the diff of PDE and ROM solutions

Calculates the L1, L2, & Linf norms of the diff between the PDE solution and the predicted solution by libROM. error in the solution if the exact solution is

Parameters
sSolver object of type HyPar
mMPI object of type MPIVariables

Definition at line 24 of file CalculateROMDiff.c.

26 {
27  HyPar* solver = (HyPar*) s;
28  MPIVariables* mpi = (MPIVariables*) m;
29  double sum = 0, global_sum = 0;
30 
31  static const double tolerance = 1e-15;
32 
33  int size = solver->npoints_local_wghosts * solver->nvars;
34  double* u_diff = (double*) calloc (size, sizeof(double));
35 
36  /* calculate solution norms (for relative error) */
37  double solution_norm[3] = {0.0,0.0,0.0};
38  /* L1 */
39  sum = ArraySumAbsnD ( solver->nvars,
40  solver->ndims,
41  solver->dim_local,
42  solver->ghosts,
43  solver->index,
44  solver->u );
45  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
46  solution_norm[0] = global_sum/((double)solver->npoints_global);
47  /* L2 */
48  sum = ArraySumSquarenD ( solver->nvars,
49  solver->ndims,
50  solver->dim_local,
51  solver->ghosts,
52  solver->index,
53  solver->u );
54  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
55  solution_norm[1] = sqrt(global_sum/((double)solver->npoints_global));
56  /* Linf */
57  sum = ArrayMaxnD ( solver->nvars,
58  solver->ndims,
59  solver->dim_local,
60  solver->ghosts,
61  solver->index
62  ,solver->u );
63  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
64  solution_norm[2] = global_sum;
65 
66  /* set u_diff to PDE solution */
67  _ArrayCopy1D_(solver->u, u_diff, size);
68  /* subtract the ROM solutions */
69  _ArrayAXPY_(solver->u_rom_predicted, -1.0, u_diff, size);
70 
71  /* calculate L1 norm of error */
72  sum = ArraySumAbsnD ( solver->nvars,
73  solver->ndims,
74  solver->dim_local,
75  solver->ghosts,
76  solver->index,
77  u_diff );
78  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
79  solver->rom_diff_norms[0] = global_sum/((double)solver->npoints_global);
80 
81  /* calculate L2 norm of error */
82  sum = ArraySumSquarenD ( solver->nvars,
83  solver->ndims,
84  solver->dim_local,
85  solver->ghosts,
86  solver->index,
87  u_diff );
88  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
89  solver->rom_diff_norms[1] = sqrt(global_sum/((double)solver->npoints_global));
90 
91  /* calculate Linf norm of error */
92  sum = ArrayMaxnD ( solver->nvars,
93  solver->ndims,
94  solver->dim_local,
95  solver->ghosts,
96  solver->index,
97  u_diff );
98  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
99  solver->rom_diff_norms[2] = global_sum;
100 
101  /* decide whether to normalize and report relative diff norms,
102  or report absolute diff norms. */
103  if ( (solution_norm[0] > tolerance)
104  && (solution_norm[1] > tolerance)
105  && (solution_norm[2] > tolerance) ) {
106  solver->rom_diff_norms[0] /= solution_norm[0];
107  solver->rom_diff_norms[1] /= solution_norm[1];
108  solver->rom_diff_norms[2] /= solution_norm[2];
109  }
110 
111  free(u_diff);
112  return 0;
113 }
HyPar::npoints_local_wghosts
int npoints_local_wghosts
Definition: hypar.h:42
HyPar::u
double * u
Definition: hypar.h:116
HyPar::npoints_global
int npoints_global
Definition: hypar.h:40
ArraySumAbsnD
INLINE double ArraySumAbsnD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:510
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
MPISum_double
int MPISum_double(double *, double *, int, void *)
Definition: MPISum.c:39
MPIMax_double
int MPIMax_double(double *, double *, int, void *)
Definition: MPIMax.c:38
HyPar::ghosts
int ghosts
Definition: hypar.h:52
MPIVariables::world
MPI_Comm world
Definition: mpivars_struct.h:37
ArraySumSquarenD
INLINE double ArraySumSquarenD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:529
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::nvars
int nvars
Definition: hypar.h:29
StdVecOps::sum
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18
HyPar::ndims
int ndims
Definition: hypar.h:26
HyPar::index
int * index
Definition: hypar.h:102
HyPar::u_rom_predicted
double * u_rom_predicted
Definition: hypar.h:403
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
ArrayMaxnD
INLINE double ArrayMaxnD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:486
_ArrayAXPY_
#define _ArrayAXPY_(x, a, y, size)
Definition: arrayfunctions.h:271
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
int OutputROMSolution ( void *  s,
int  nsims,
double  a_time 
)

Write ROM solutions to file

Write out the ROM solution to file

Parameters
sArray of simulation objects of type SimulationObject
nsimsNumber of simulation objects
a_timeCurrent simulation time

Definition at line 24 of file OutputROMSolution.cpp.

27 {
28  SimulationObject* simobj = (SimulationObject*) s;
29  int ns;
30 
31  for (ns = 0; ns < nsims; ns++) {
32 
33  HyPar* solver = &(simobj[ns].solver);
34  MPIVariables* mpi = &(simobj[ns].mpi);
35 
36  if ((!solver->WriteOutput) && (strcmp(solver->plot_solution,"yes"))) continue;
37 
38  char fname_root[_MAX_STRING_SIZE_];
39  strcpy(fname_root, solver->op_rom_fname_root);
40 
41  if (nsims > 1) {
42  char index[_MAX_STRING_SIZE_];
43  GetStringFromInteger(ns, index, (int)log10(nsims)+1);
44  strcat(fname_root, "_");
45  strcat(fname_root, index);
46  }
47 
48  WriteArray( solver->ndims,
49  solver->nvars,
50  solver->dim_global,
51  solver->dim_local,
52  solver->ghosts,
53  solver->x,
54  solver->u_rom_predicted,
55  solver,
56  mpi,
57  fname_root );
58 
59  if (!strcmp(solver->plot_solution, "yes")) {
60  PlotArray( solver->ndims,
61  solver->nvars,
62  solver->dim_global,
63  solver->dim_local,
64  solver->ghosts,
65  solver->x,
66  solver->u,
67  a_time,
68  solver,
69  mpi,
70  fname_root );
71  }
72 
73  /* increment the index string, if required */
74  if ((!strcmp(solver->output_mode,"serial")) && (!strcmp(solver->op_overwrite,"no"))) {
75  IncrementFilenameIndex(solver->filename_index,solver->index_length);
76  }
77 
78  }
79 
80  return 0;
81 }
HyPar::dim_global
int * dim_global
Definition: hypar.h:33
HyPar::u
double * u
Definition: hypar.h:116
GetStringFromInteger
void GetStringFromInteger(int, char *, int)
Definition: CommonFunctions.c:15
HyPar::plot_solution
char plot_solution[_MAX_STRING_SIZE_]
Definition: hypar.h:194
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
HyPar::ghosts
int ghosts
Definition: hypar.h:52
HyPar::filename_index
char * filename_index
Definition: hypar.h:197
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
HyPar::op_overwrite
char op_overwrite[_MAX_STRING_SIZE_]
Definition: hypar.h:191
HyPar::op_rom_fname_root
char op_rom_fname_root[_MAX_STRING_SIZE_]
Definition: hypar.h:407
HyPar::nvars
int nvars
Definition: hypar.h:29
PlotArray
int PlotArray(int, int, int *, int *, int, double *, double *, double, void *, void *, char *)
Definition: PlotArray.c:31
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
HyPar::ndims
int ndims
Definition: hypar.h:26
HyPar::u_rom_predicted
double * u_rom_predicted
Definition: hypar.h:403
HyPar::index_length
int index_length
Definition: hypar.h:199
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar::x
double * x
Definition: hypar.h:107
HyPar::output_mode
char output_mode[_MAX_STRING_SIZE_]
Definition: hypar.h:183
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
IncrementFilenameIndex
void IncrementFilenameIndex(char *f, int len)
Definition: IncrementFilename.c:46
HyPar::WriteOutput
int(* WriteOutput)(int, int, int *, double *, double *, char *, int *)
Definition: hypar.h:211
int Solve ( void *  s,
int  nsims,
int  rank,
int  nproc 
)

This function integrates the semi-discrete ODE (obtained from discretizing the PDE in space) using natively implemented time integration methods. It initializes the time integration object, iterates the simulation for the required number of time steps, and calculates the errors. After the specified number of iterations, it writes out some information to the screen and the solution to a file.

Parameters
sArray of simulation objects of type SimulationObject
nsimsnumber of simulation objects
rankMPI rank of this process
nprocNumber of MPI processes

Definition at line 37 of file Solve.cpp.

42 {
43  SimulationObject* sim = (SimulationObject*) s;
44 
45  /* make sure none of the simulation objects sent in the array
46  * are "barebones" type */
47  for (int ns = 0; ns < nsims; ns++) {
48  if (sim[ns].is_barebones == 1) {
49  fprintf(stderr, "Error in Solve(): simulation object %d on rank %d is barebones!\n",
50  ns, rank );
51  return 1;
52  }
53  }
54 
55  /* write out iblank to file for visualization */
56  for (int ns = 0; ns < nsims; ns++) {
57  if (sim[ns].solver.flag_ib) {
58 
59  char fname_root[_MAX_STRING_SIZE_] = "iblank";
60  if (nsims > 1) {
61  char index[_MAX_STRING_SIZE_];
62  GetStringFromInteger(ns, index, (int)log10((nsims)+1));
63  strcat(fname_root, "_");
64  strcat(fname_root, index);
65  }
66 
67  WriteArray( sim[ns].solver.ndims,
68  1,
69  sim[ns].solver.dim_global,
70  sim[ns].solver.dim_local,
71  sim[ns].solver.ghosts,
72  sim[ns].solver.x,
73  sim[ns].solver.iblank,
74  &(sim[ns].solver),
75  &(sim[ns].mpi),
76  fname_root );
77  }
78  }
79 
80 #ifdef with_librom
81  if (!rank) printf("Setting up libROM interface.\n");
82  libROMInterface rom_interface( sim, nsims, rank, nproc, sim[0].solver.dt );
83  const std::string& rom_mode( rom_interface.mode() );
84  std::vector<double> op_times_arr(0);
85 #endif
86 
87 #ifdef with_librom
88  if ((rom_mode == _ROM_MODE_TRAIN_) || (rom_mode == _ROM_MODE_NONE_)) {
89 #endif
90  /* Define and initialize the time-integration object */
91  TimeIntegration TS;
92  if (!rank) printf("Setting up time integration.\n");
93  TimeInitialize(sim, nsims, rank, nproc, &TS);
94  double ti_runtime = 0.0;
95 
96  if (!rank) printf("Solving in time (from %d to %d iterations)\n",TS.restart_iter,TS.n_iter);
97  for (TS.iter = TS.restart_iter; TS.iter < TS.n_iter; TS.iter++) {
98 
99  /* Write initial solution to file if this is the first iteration */
100  if (!TS.iter) {
101  for (int ns = 0; ns < nsims; ns++) {
102  if (sim[ns].solver.PhysicsOutput) {
103  sim[ns].solver.PhysicsOutput( &(sim[ns].solver),
104  &(sim[ns].mpi),
105  TS.waqt );
106  }
107  }
108  OutputSolution(sim, nsims, TS.waqt);
109 #ifdef with_librom
110  op_times_arr.push_back(TS.waqt);
111 #endif
112  }
113 
114 #ifdef with_librom
115  if ((rom_mode == _ROM_MODE_TRAIN_) && (TS.iter%rom_interface.samplingFrequency() == 0)) {
116  rom_interface.takeSample( sim, TS.waqt );
117  }
118 #endif
119 
120  /* Call pre-step function */
121  TimePreStep (&TS);
122 #ifdef compute_rhs_operators
123  /* compute and write (to file) matrix operators representing the right-hand side */
124 // if (((TS.iter+1)%solver->file_op_iter == 0) || (!TS.iter))
125 // { ComputeRHSOperators(solver,mpi,TS.waqt);
126 #endif
127 
128  /* Step in time */
129  TimeStep (&TS);
130 
131  /* Call post-step function */
132  TimePostStep (&TS);
133 
134  ti_runtime += TS.iter_wctime;
135 
136  /* Print information to screen */
137  TimePrintStep(&TS);
138 
139  /* Write intermediate solution to file */
140  if ( ((TS.iter+1)%sim[0].solver.file_op_iter == 0)
141  && ((TS.iter+1) < TS.n_iter) ) {
142  for (int ns = 0; ns < nsims; ns++) {
143  if (sim[ns].solver.PhysicsOutput) {
144  sim[ns].solver.PhysicsOutput( &(sim[ns].solver),
145  &(sim[ns].mpi),
146  TS.waqt );
147  }
148  }
149  OutputSolution(sim, nsims, TS.waqt);
150 #ifdef with_librom
151  op_times_arr.push_back(TS.waqt);
152 #endif
153  }
154 
155  }
156 
157  double t_final = TS.waqt;
158  TimeCleanup(&TS);
159 
160  if (!rank) {
161  printf( "Completed time integration (Final time: %f), total wctime: %f (seconds).\n",
162  t_final, ti_runtime );
163  if (nsims > 1) printf("\n");
164  }
165 
166  /* calculate error if exact solution has been provided */
167  for (int ns = 0; ns < nsims; ns++) {
168  CalculateError(&(sim[ns].solver),
169  &(sim[ns].mpi) );
170  }
171 
172  /* write a final solution file */
173  for (int ns = 0; ns < nsims; ns++) {
174  if (sim[ns].solver.PhysicsOutput) {
175  sim[ns].solver.PhysicsOutput( &(sim[ns].solver),
176  &(sim[ns].mpi),
177  t_final );
178  }
179  }
180  OutputSolution(sim, nsims, t_final);
181 
182 #ifdef with_librom
183  op_times_arr.push_back(TS.waqt);
184 
185  for (int ns = 0; ns < nsims; ns++) {
186  ResetFilenameIndex( sim[ns].solver.filename_index,
187  sim[ns].solver.index_length );
188  }
189 
190  if (rom_interface.mode() == _ROM_MODE_TRAIN_) {
191 
192  rom_interface.train();
193  if (!rank) printf("libROM: total training wallclock time: %f (seconds).\n",
194  rom_interface.trainWallclockTime() );
195 
196  double total_rom_predict_time = 0;
197  for (int iter = 0; iter < op_times_arr.size(); iter++) {
198 
199  double waqt = op_times_arr[iter];
200 
201  rom_interface.predict(sim, waqt);
202  if (!rank) printf( "libROM: Predicted solution at time %1.4e using ROM, wallclock time: %f.\n",
203  waqt, rom_interface.predictWallclockTime() );
204  total_rom_predict_time += rom_interface.predictWallclockTime();
205 
206  /* calculate diff between ROM and PDE solutions */
207  if (iter == (op_times_arr.size()-1)) {
208  if (!rank) printf("libROM: Calculating diff between PDE and ROM solutions.\n");
209  for (int ns = 0; ns < nsims; ns++) {
210  CalculateROMDiff( &(sim[ns].solver),
211  &(sim[ns].mpi) );
212  }
213  }
214  /* write the ROM solution to file */
215  OutputROMSolution(sim, nsims,waqt);
216 
217  }
218 
219  if (!rank) {
220  printf( "libROM: total prediction/query wallclock time: %f (seconds).\n",
221  total_rom_predict_time );
222  }
223 
224  rom_interface.saveROM();
225 
226  } else {
227 
228  for (int ns = 0; ns < nsims; ns++) {
229  sim[ns].solver.rom_diff_norms[0]
230  = sim[ns].solver.rom_diff_norms[1]
231  = sim[ns].solver.rom_diff_norms[2]
232  = -1;
233  }
234 
235  }
236 
237  } else if (rom_mode == _ROM_MODE_PREDICT_) {
238 
239  for (int ns = 0; ns < nsims; ns++) {
240  sim[ns].solver.rom_diff_norms[0]
241  = sim[ns].solver.rom_diff_norms[1]
242  = sim[ns].solver.rom_diff_norms[2]
243  = -1;
244  strcpy(sim[ns].solver.ConservationCheck,"no");
245  }
246 
247  rom_interface.loadROM();
248  rom_interface.projectInitialSolution(sim);
249 
250  {
251  int start_iter = sim[0].solver.restart_iter;
252  int n_iter = sim[0].solver.n_iter;
253  double dt = sim[0].solver.dt;
254 
255  double cur_time = start_iter * dt;
256  op_times_arr.push_back(cur_time);
257 
258  for (int iter = start_iter; iter < n_iter; iter++) {
259  cur_time += dt;
260  if ( ( (iter+1)%sim[0].solver.file_op_iter == 0)
261  && ( (iter+1) < n_iter) ) {
262  op_times_arr.push_back(cur_time);
263  }
264  }
265 
266  double t_final = n_iter*dt;
267  op_times_arr.push_back(t_final);
268  }
269 
270  double total_rom_predict_time = 0;
271  for (int iter = 0; iter < op_times_arr.size(); iter++) {
272 
273  double waqt = op_times_arr[iter];
274 
275  rom_interface.predict(sim, waqt);
276  if (!rank) printf( "libROM: Predicted solution at time %1.4e using ROM, wallclock time: %f.\n",
277  waqt, rom_interface.predictWallclockTime() );
278  total_rom_predict_time += rom_interface.predictWallclockTime();
279 
280  /* write the solution to file */
281  for (int ns = 0; ns < nsims; ns++) {
282  if (sim[ns].solver.PhysicsOutput) {
283  sim[ns].solver.PhysicsOutput( &(sim[ns].solver),
284  &(sim[ns].mpi),
285  waqt );
286  }
287  }
288  OutputSolution(sim, nsims, waqt);
289 
290  }
291 
292  /* calculate error if exact solution has been provided */
293  for (int ns = 0; ns < nsims; ns++) {
294  CalculateError(&(sim[ns].solver),
295  &(sim[ns].mpi) );
296  }
297 
298  if (!rank) {
299  printf( "libROM: total prediction/query wallclock time: %f (seconds).\n",
300  total_rom_predict_time );
301  }
302 
303  }
304 #endif
305 
306  return 0;
307 }
TimePostStep
int TimePostStep(void *)
Definition: TimePostStep.c:28
TimePreStep
int TimePreStep(void *)
Definition: TimePreStep.c:22
CalculateError
int CalculateError(void *s, void *m)
Definition: CalculateError.c:25
HyPar::dim_global
int * dim_global
Definition: hypar.h:33
HyPar::iblank
double * iblank
Definition: hypar.h:436
TimeStep
int TimeStep(void *)
Definition: TimeStep.c:13
GetStringFromInteger
void GetStringFromInteger(int, char *, int)
Definition: CommonFunctions.c:15
TimeIntegration
Structure of variables/parameters and function pointers for time integration.
Definition: timeintegration_struct.h:29
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
HyPar::PhysicsOutput
int(* PhysicsOutput)(void *, void *, double)
Definition: hypar.h:347
TimeCleanup
int TimeCleanup(void *)
Definition: TimeCleanup.c:18
OutputROMSolution
int OutputROMSolution(void *s, int nsims, double a_time)
Definition: OutputROMSolution.cpp:24
HyPar::ghosts
int ghosts
Definition: hypar.h:52
HyPar::dt
double dt
Definition: hypar.h:67
ResetFilenameIndex
void ResetFilenameIndex(char *f, int len)
Definition: IncrementFilename.c:64
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
libROMInterface
Class implementing interface with libROM.
Definition: librom_interface.h:45
CalculateROMDiff
int CalculateROMDiff(void *s, void *m)
Definition: CalculateROMDiff.c:24
HyPar::file_op_iter
int file_op_iter
Definition: hypar.h:171
OutputSolution
int OutputSolution(void *, int, double)
Definition: OutputSolution.cpp:27
_ROM_MODE_PREDICT_
#define _ROM_MODE_PREDICT_
Definition: librom_interface.h:16
HyPar::n_iter
int n_iter
Definition: hypar.h:55
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
TimePrintStep
int TimePrintStep(void *)
Definition: TimePrintStep.c:16
_ROM_MODE_NONE_
#define _ROM_MODE_NONE_
Definition: librom_interface.h:12
HyPar::index_length
int index_length
Definition: hypar.h:199
HyPar::x
double * x
Definition: hypar.h:107
_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::restart_iter
int restart_iter
Definition: hypar.h:58
TimeInitialize
int TimeInitialize(void *, int, int, int, void *)
Definition: TimeInitialize.c:28