HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
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

◆ ComputeRHSOperators()

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 }
absolute
#define absolute(a)
Definition: math_ops.h:32
HyPar::nvars
int nvars
Definition: hypar.h:29
IERR
#define IERR
Definition: basic.h:16
HyPar::filename_index
char * filename_index
Definition: hypar.h:197
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
MPIExchangeBoundariesnD
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Definition: MPIExchangeBoundariesnD.c:42
HyPar::ApplyBoundaryConditions
int(* ApplyBoundaryConditions)(void *, void *, double *, double *, double)
Definition: hypar.h:214
HyPar::u
double * u
Definition: hypar.h:116
HyPar::npoints_local
int npoints_local
Definition: hypar.h:42
HyPar::FdFFunction
int(* FdFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:286
HyPar::ApplyIBConditions
int(* ApplyIBConditions)(void *, void *, double *, double)
Definition: hypar.h:217
HyPar::SFunction
int(* SFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:317
_ArrayIndexnD_
#define _ArrayIndexnD_(N, index, imax, i, ghost)
Definition: arrayfunctions.h:21
HyPar::Upwind
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
HyPar::SourceFunction
int(* SourceFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:259
HyPar::ndims
int ndims
Definition: hypar.h:26
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::FFunction
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
_ArrayAXPY_
#define _ArrayAXPY_(x, a, y, size)
Definition: arrayfunctions.h:271
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
ArrayCopynD
INLINE int ArrayCopynD(int, const double *, double *, int *, int, int, int *, int)
Definition: arrayfunctions.h:410
HyPar::UpwindFdF
int(* UpwindFdF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:307
HyPar::flag_fdf_specified
int flag_fdf_specified
Definition: hypar.h:290
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
HyPar::dFFunction
int(* dFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:280
HyPar::UpwinddF
int(* UpwinddF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:300
HyPar::ghosts
int ghosts
Definition: hypar.h:52
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
_ArrayScale1D_
#define _ArrayScale1D_(x, a, size)
Definition: arrayfunctions.h:190
HyPar::npoints_local_wghosts
int npoints_local_wghosts
Definition: hypar.h:42
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319

◆ CalculateError()

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::nvars
int nvars
Definition: hypar.h:29
IERR
#define IERR
Definition: basic.h:16
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
MPISum_double
int MPISum_double(double *, double *, int, void *)
Definition: MPISum.c:39
HyPar::nsims
int nsims
Definition: hypar.h:64
HyPar::error
double error[3]
Definition: hypar.h:371
StdVecOps::sum
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18
HyPar::u
double * u
Definition: hypar.h:116
HyPar::ndims
int ndims
Definition: hypar.h:26
ExactSolution
int ExactSolution(void *, void *, double *, char *, int *)
Definition: ExactSolution.c:16
ArraySumAbsnD
INLINE double ArraySumAbsnD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:510
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
MPIMax_double
int MPIMax_double(double *, double *, int, void *)
Definition: MPIMax.c:38
_ArrayAXPY_
#define _ArrayAXPY_(x, a, y, size)
Definition: arrayfunctions.h:271
HyPar::npoints_global
int npoints_global
Definition: hypar.h:40
MPIVariables::world
MPI_Comm world
Definition: mpivars_struct.h:37
HyPar::index
int * index
Definition: hypar.h:102
GetStringFromInteger
void GetStringFromInteger(int, char *, int)
Definition: CommonFunctions.c:15
HyPar::my_idx
int my_idx
Definition: hypar.h:61
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
ArrayMaxnD
INLINE double ArrayMaxnD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:486
HyPar::ghosts
int ghosts
Definition: hypar.h:52
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
ArraySumSquarenD
INLINE double ArraySumSquarenD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:529
TimeError
int TimeError(void *, void *, double *)
Definition: TimeError.c:27

◆ OutputSolution()

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::plot_solution
char plot_solution[_MAX_STRING_SIZE_]
Definition: hypar.h:194
HyPar::nvars
int nvars
Definition: hypar.h:29
IERR
#define IERR
Definition: basic.h:16
HyPar::filename_index
char * filename_index
Definition: hypar.h:197
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
PlotArray
int PlotArray(int, int, int *, int *, int, double *, double *, double, void *, void *, char *)
Definition: PlotArray.c:31
TimeGetAuxSolutions
int TimeGetAuxSolutions(int *, double **, void *, int, int)
Definition: TimeGetAuxSolutions.c:29
HyPar::u
double * u
Definition: hypar.h:116
HyPar::output_mode
char output_mode[_MAX_STRING_SIZE_]
Definition: hypar.h:183
HyPar::x
double * x
Definition: hypar.h:107
HyPar::size_x
int size_x
Definition: hypar.h:48
HyPar::gpu_x
double * gpu_x
Definition: hypar.h:457
HyPar::aux_op_fname_root
char aux_op_fname_root[_MAX_STRING_SIZE_]
Definition: hypar.h:208
HyPar::ndims
int ndims
Definition: hypar.h:26
HyPar::ndof_cells_wghosts
int ndof_cells_wghosts
Definition: hypar.h:47
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_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
WriteArray
int WriteArray(int, int, int *, int *, int, double *, double *, void *, void *, char *)
Definition: WriteArray.c:27
GetStringFromInteger
void GetStringFromInteger(int, char *, int)
Definition: CommonFunctions.c:15
HyPar::index_length
int index_length
Definition: hypar.h:199
HyPar::op_fname_root
char op_fname_root[_MAX_STRING_SIZE_]
Definition: hypar.h:206
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
gpuMemcpy
void gpuMemcpy(void *, const void *, size_t, enum gpuMemcpyKind)
IncrementFilenameIndex
void IncrementFilenameIndex(char *, int)
Definition: IncrementFilename.c:46
HyPar::WriteOutput
int(* WriteOutput)(int, int, int *, double *, double *, char *, int *)
Definition: hypar.h:211
HyPar::ghosts
int ghosts
Definition: hypar.h:52
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar::use_gpu
int use_gpu
Definition: hypar.h:449
HyPar::npoints_local_wghosts
int npoints_local_wghosts
Definition: hypar.h:42
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
HyPar::dim_global
int * dim_global
Definition: hypar.h:33

◆ ResetFilenameIndex()

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 }

◆ CalculateROMDiff()

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::nvars
int nvars
Definition: hypar.h:29
MPISum_double
int MPISum_double(double *, double *, int, void *)
Definition: MPISum.c:39
StdVecOps::sum
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18
HyPar::u
double * u
Definition: hypar.h:116
HyPar::ndims
int ndims
Definition: hypar.h:26
ArraySumAbsnD
INLINE double ArraySumAbsnD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:510
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
MPIMax_double
int MPIMax_double(double *, double *, int, void *)
Definition: MPIMax.c:38
_ArrayAXPY_
#define _ArrayAXPY_(x, a, y, size)
Definition: arrayfunctions.h:271
HyPar::npoints_global
int npoints_global
Definition: hypar.h:40
MPIVariables::world
MPI_Comm world
Definition: mpivars_struct.h:37
HyPar::index
int * index
Definition: hypar.h:102
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
HyPar::rom_diff_norms
double rom_diff_norms[3]
Definition: hypar.h:405
ArrayMaxnD
INLINE double ArrayMaxnD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:486
HyPar::ghosts
int ghosts
Definition: hypar.h:52
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar::npoints_local_wghosts
int npoints_local_wghosts
Definition: hypar.h:42
HyPar::u_rom_predicted
double * u_rom_predicted
Definition: hypar.h:403
ArraySumSquarenD
INLINE double ArraySumSquarenD(int, int, int *, int, int *, double *)
Definition: arrayfunctions.h:529
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319

◆ OutputROMSolution()

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::plot_solution
char plot_solution[_MAX_STRING_SIZE_]
Definition: hypar.h:194
HyPar::nvars
int nvars
Definition: hypar.h:29
HyPar::filename_index
char * filename_index
Definition: hypar.h:197
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
PlotArray
int PlotArray(int, int, int *, int *, int, double *, double *, double, void *, void *, char *)
Definition: PlotArray.c:31
HyPar::op_rom_fname_root
char op_rom_fname_root[_MAX_STRING_SIZE_]
Definition: hypar.h:407
HyPar::u
double * u
Definition: hypar.h:116
HyPar::output_mode
char output_mode[_MAX_STRING_SIZE_]
Definition: hypar.h:183
HyPar::x
double * x
Definition: hypar.h:107
IncrementFilenameIndex
void IncrementFilenameIndex(char *, int)
Definition: IncrementFilename.c:46
HyPar::ndims
int ndims
Definition: hypar.h:26
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
HyPar::op_overwrite
char op_overwrite[_MAX_STRING_SIZE_]
Definition: hypar.h:191
WriteArray
int WriteArray(int, int, int *, int *, int, double *, double *, void *, void *, char *)
Definition: WriteArray.c:27
GetStringFromInteger
void GetStringFromInteger(int, char *, int)
Definition: CommonFunctions.c:15
HyPar::index_length
int index_length
Definition: hypar.h:199
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
HyPar::WriteOutput
int(* WriteOutput)(int, int, int *, double *, double *, char *, int *)
Definition: hypar.h:211
HyPar::ghosts
int ghosts
Definition: hypar.h:52
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar::u_rom_predicted
double * u_rom_predicted
Definition: hypar.h:403
HyPar::dim_global
int * dim_global
Definition: hypar.h:33

◆ Solve()

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 }
CalculateError
int CalculateError(void *, void *)
Definition: CalculateError.c:25
TimeIntegration
Structure of variables/parameters and function pointers for time integration.
Definition: timeintegration_struct.h:29
TimeInitialize
int TimeInitialize(void *, int, int, int, void *)
Definition: TimeInitialize.c:28
HyPar::PhysicsOutput
int(* PhysicsOutput)(void *, void *, double)
Definition: hypar.h:347
TimePreStep
int TimePreStep(void *)
Definition: TimePreStep.c:22
SimulationObject
Structure defining a simulation.
Definition: simulation_object.h:23
OutputSolution
int OutputSolution(void *, int, double)
Definition: OutputSolution.cpp:27
HyPar::n_iter
int n_iter
Definition: hypar.h:55
_ROM_MODE_PREDICT_
#define _ROM_MODE_PREDICT_
Definition: librom_interface.h:16
TimePrintStep
int TimePrintStep(void *)
Definition: TimePrintStep.c:16
HyPar::restart_iter
int restart_iter
Definition: hypar.h:58
OutputROMSolution
int OutputROMSolution(void *, int, double)
Definition: OutputROMSolution.cpp:24
TimePostStep
int TimePostStep(void *)
Definition: TimePostStep.c:28
_ROM_MODE_NONE_
#define _ROM_MODE_NONE_
Definition: librom_interface.h:12
_MAX_STRING_SIZE_
#define _MAX_STRING_SIZE_
Definition: basic.h:14
_ROM_MODE_TRAIN_
#define _ROM_MODE_TRAIN_
Definition: librom_interface.h:14
WriteArray
int WriteArray(int, int, int *, int *, int, double *, double *, void *, void *, char *)
Definition: WriteArray.c:27
GetStringFromInteger
void GetStringFromInteger(int, char *, int)
Definition: CommonFunctions.c:15
TimeCleanup
int TimeCleanup(void *)
Definition: TimeCleanup.c:18
HyPar::dt
double dt
Definition: hypar.h:67
HyPar::rom_diff_norms
double rom_diff_norms[3]
Definition: hypar.h:405
TimeStep
int TimeStep(void *)
Definition: TimeStep.c:13
libROMInterface
Class implementing interface with libROM.
Definition: librom_interface.h:45
ResetFilenameIndex
void ResetFilenameIndex(char *, int)
Definition: IncrementFilename.c:64
CalculateROMDiff
int CalculateROMDiff(void *, void *)
Definition: CalculateROMDiff.c:24