25 #include "grins_config.h"
39 #include "libmesh/parallel.h"
40 #include "libmesh/exact_solution.h"
44 initial_values(
const libMesh::Point& p,
const libMesh::Parameters ¶ms,
45 const std::string& system_name,
const std::string& unknown_name );
47 int run(
int argc,
char* argv[],
const GetPot& input );
49 int main(
int argc,
char* argv[])
51 #ifdef USE_GRVY_TIMERS
52 GRVY::GRVY_Timer_Class grvy_timer;
53 grvy_timer.Init(
"GRINS Timer");
60 std::cerr <<
"Error: Must specify libMesh input file and exact solution file." << std::endl;
65 std::string libMesh_input_filename = argv[1];
68 GetPot libMesh_inputfile( libMesh_input_filename );
72 std::string chem_lib = libMesh_inputfile(
"Physics/ReactingLowMachNavierStokes/thermochemistry_library",
"DIE!");
74 if( chem_lib == std::string(
"cantera") )
76 #ifdef GRINS_HAVE_CANTERA
77 return_flag =
run(argc,argv,libMesh_inputfile);
82 else if( chem_lib == std::string(
"antioch") )
84 #ifdef GRINS_HAVE_ANTIOCH
85 return_flag =
run(argc,argv,libMesh_inputfile);
98 int run(
int argc,
char* argv[],
const GetPot& input )
100 #ifdef USE_GRVY_TIMERS
101 grvy_timer.BeginTimer(
"Initialize Solver");
105 libMesh::LibMeshInit libmesh_init(argc, argv);
111 libmesh_init.comm() );
114 std::string restart_file = input(
"restart-options/restart_file",
"none" );
116 if( restart_file ==
"none" )
119 std::string system_name = input(
"screen-options/system_name",
"GRINS" );
120 std::tr1::shared_ptr<libMesh::EquationSystems> es = grins.get_equation_system();
121 const libMesh::System& system = es->get_system(system_name);
123 libMesh::Parameters ¶ms = es->parameters;
125 libMesh::Real& w_N2 = params.set<libMesh::Real>(
"w_N2" );
126 w_N2 = input(
"Physics/ReactingLowMachNavierStokes/bound_species_3", 0.0, 0 );
128 libMesh::Real& w_N = params.set<libMesh::Real>(
"w_N" );
129 w_N = input(
"Physics/ReactingLowMachNavierStokes/bound_species_3", 0.0, 1 );
134 #ifdef USE_GRVY_TIMERS
135 grvy_timer.EndTimer(
"Initialize Solver");
138 grins.attach_grvy_timer( &grvy_timer );
143 #ifdef USE_GRVY_TIMERS
144 grvy_timer.Finalize();
146 if( Parallel::Communicator_World.rank() == 0 ) grvy_timer.Summarize();
150 std::tr1::shared_ptr<libMesh::EquationSystems>
151 es = grins.get_equation_system();
156 libMesh::ExactSolution exact_sol(*es);
158 libMesh::EquationSystems es_ref( es->get_mesh() );
161 std::string solution_file = std::string(argv[2]);
162 es_ref.read( solution_file );
164 exact_sol.attach_reference_solution( &es_ref );
167 exact_sol.compute_error(
"GRINS",
"u");
168 exact_sol.compute_error(
"GRINS",
"v");
170 if( (es->get_mesh()).mesh_dimension() == 3 )
171 exact_sol.compute_error(
"GRINS",
"w");
173 exact_sol.compute_error(
"GRINS",
"p");
174 exact_sol.compute_error(
"GRINS",
"T");
175 exact_sol.compute_error(
"GRINS",
"w_N2");
176 exact_sol.compute_error(
"GRINS",
"w_N");
178 double u_l2error = exact_sol.l2_error(
"GRINS",
"u");
179 double u_h1error = exact_sol.h1_error(
"GRINS",
"u");
181 double v_l2error = exact_sol.l2_error(
"GRINS",
"v");
182 double v_h1error = exact_sol.h1_error(
"GRINS",
"v");
184 double p_l2error = exact_sol.l2_error(
"GRINS",
"p");
185 double p_h1error = exact_sol.h1_error(
"GRINS",
"p");
187 double T_l2error = exact_sol.l2_error(
"GRINS",
"T");
188 double T_h1error = exact_sol.h1_error(
"GRINS",
"T");
190 double wN_l2error = exact_sol.l2_error(
"GRINS",
"w_N");
191 double wN_h1error = exact_sol.h1_error(
"GRINS",
"w_N");
193 double wN2_l2error = exact_sol.l2_error(
"GRINS",
"w_N2");
194 double wN2_h1error = exact_sol.h1_error(
"GRINS",
"w_N2");
196 double w_l2error = 0.0,
199 if( (es->get_mesh()).mesh_dimension() == 3 )
201 w_l2error = exact_sol.l2_error(
"GRINS",
"w");
202 w_h1error = exact_sol.h1_error(
"GRINS",
"w");
209 if( u_l2error > tol || u_h1error > tol ||
210 v_l2error > tol || v_h1error > tol ||
211 w_l2error > tol || w_h1error > tol ||
212 p_l2error > tol || p_h1error > tol ||
213 T_l2error > tol || T_h1error > tol ||
214 wN_l2error > tol || wN_h1error > tol ||
215 wN2_l2error > tol || wN2_h1error > tol )
219 std::cout <<
"Tolerance exceeded for thermally driven flow test." << std::endl
220 <<
"tolerance = " << tol << std::endl
221 <<
"u l2 error = " << u_l2error << std::endl
222 <<
"u h1 error = " << u_h1error << std::endl
223 <<
"v l2 error = " << v_l2error << std::endl
224 <<
"v h1 error = " << v_h1error << std::endl
225 <<
"w l2 error = " << w_l2error << std::endl
226 <<
"w h1 error = " << w_h1error << std::endl
227 <<
"p l2 error = " << p_l2error << std::endl
228 <<
"p h1 error = " << p_h1error << std::endl
229 <<
"T l2 error = " << T_l2error << std::endl
230 <<
"T h1 error = " << T_h1error << std::endl
231 <<
"w_N l2 error = " << wN_l2error << std::endl
232 <<
"w_N h1 error = " << wN_h1error << std::endl
233 <<
"w_N2 l2 error = " << wN2_l2error << std::endl
234 <<
"w_N2 h1 error = " << wN2_h1error << std::endl;
242 const std::string& ,
const std::string& unknown_name )
244 libMesh::Real value = 0.0;
246 if( unknown_name ==
"w_N2" )
247 value = params.get<libMesh::Real>(
"w_N2");
249 else if( unknown_name ==
"w_N" )
250 value = params.get<libMesh::Real>(
"w_N");
252 else if( unknown_name ==
"T" )
255 else if( unknown_name ==
"u" )
257 const libMesh::Real y = p(1);
258 value = 1.0*(-y*y+1);
int main(int argc, char *argv[])
int run(int argc, char *argv[], const GetPot &input)
libMesh::Real initial_values(const libMesh::Point &p, const libMesh::Parameters ¶ms, const std::string &system_name, const std::string &unknown_name)