26 #include "grins_config.h"
38 #include "libmesh/dirichlet_boundaries.h"
39 #include "libmesh/dof_map.h"
40 #include "libmesh/fe_base.h"
41 #include "libmesh/fe_interface.h"
42 #include "libmesh/mesh_function.h"
43 #include "libmesh/serial_mesh.h"
44 #include "libmesh/edge_edge2.h"
45 #include "libmesh/mesh_generation.h"
46 #include "libmesh/linear_implicit_system.h"
47 #include "libmesh/gmv_io.h"
48 #include "libmesh/exact_solution.h"
51 #ifdef GRINS_HAVE_GRVY
56 const std::string& system_name,
57 const std::string& var,
58 const std::string& norm,
73 std::multimap< GRINS::PhysicsName, GRINS::DBCContainer >
build_dirichlet( );
86 { this->_initialized =
true; }
88 virtual libMesh::Number
operator() (
const libMesh::Point&,
const libMesh::Real = 0)
89 { libmesh_not_implemented(); }
92 const libMesh::Real t,
93 libMesh::DenseVector<libMesh::Number>& output)
99 libMesh::Point p_copy(p);
101 p_copy(0) = p_copy(1);
107 p_copy(0) = 1 - p_copy(0);
109 p_copy(0) = 2*p_copy(0);
111 libMesh::DenseVector<libMesh::Number> u_nu_values;
114 output(0) = u_nu_values(0)/21.995539;
117 virtual libMesh::AutoPtr<libMesh::FunctionBase<libMesh::Number> >
clone()
const
130 { this->_initialized =
true; }
132 virtual libMesh::Number
operator() (
const libMesh::Point&,
const libMesh::Real = 0)
133 { libmesh_not_implemented(); }
136 const libMesh::Real t,
137 libMesh::DenseVector<libMesh::Number>& output)
143 libMesh::Point p_copy(p);
145 p_copy(0) = p_copy(1);
151 p_copy(0) = 1 - p_copy(0);
153 p_copy(0) = 2*p_copy(0);
155 libMesh::DenseVector<libMesh::Number> u_nu_values;
158 output(0) = u_nu_values(1)/(2.0*21.995539);
161 virtual libMesh::AutoPtr<libMesh::FunctionBase<libMesh::Number> >
clone()
const
169 int main(
int argc,
char* argv[])
172 #ifdef GRINS_USE_GRVY_TIMERS
173 GRVY::GRVY_Timer_Class grvy_timer;
174 grvy_timer.Init(
"GRINS Timer");
181 std::cerr <<
"Error: Must specify libMesh input file." << std::endl;
186 std::string libMesh_input_filename = argv[1];
189 GetPot libMesh_inputfile( libMesh_input_filename );
191 #ifdef GRINS_USE_GRVY_TIMERS
192 grvy_timer.BeginTimer(
"Initialize Solver");
196 libMesh::LibMeshInit libmesh_init(argc, argv);
199 libMesh::SerialMesh mesh(libmesh_init.comm());
201 GetPot command_line(argc,argv);
203 std::string oned_mesh = command_line(
"mesh-1d",
"DIE!");
204 mesh.read(oned_mesh);
207 libMesh::EquationSystems equation_systems (mesh);
209 std::string oned_data = command_line(
"data-1d",
"DIE!");
210 equation_systems.read(oned_data, libMesh::XdrMODE::READ,
211 libMesh::EquationSystems::READ_HEADER |
212 libMesh::EquationSystems::READ_DATA |
213 libMesh::EquationSystems::READ_ADDITIONAL_DATA);
216 libMesh::System & turbulent_bc_system = equation_systems.get_system<libMesh::System>(
"flow");
220 turbulent_bc_system.update();
223 equation_systems.print_info();
226 libMesh::AutoPtr<libMesh::MeshFunction> turbulent_bc_values;
228 libMesh::AutoPtr<libMesh::NumericVector<libMesh::Number> > turbulent_bc_soln = libMesh::NumericVector<libMesh::Number>::build(equation_systems.comm());
230 std::vector<libMesh::Number> flow_soln;
232 turbulent_bc_system.update_global_solution(flow_soln);
234 turbulent_bc_soln->init(turbulent_bc_system.solution->size(),
true, libMesh::SERIAL);
236 (*turbulent_bc_soln) = flow_soln;
238 std::vector<unsigned int>turbulent_bc_system_variables;
239 turbulent_bc_system_variables.push_back(0);
240 turbulent_bc_system_variables.push_back(1);
242 turbulent_bc_values = libMesh::AutoPtr<libMesh::MeshFunction>
243 (
new libMesh::MeshFunction(equation_systems,
245 turbulent_bc_system.get_dof_map(),
246 turbulent_bc_system_variables ));
248 turbulent_bc_values->init();
252 std::tr1::shared_ptr<TurbulentBCFactory> bc_factory(
new TurbulentBCFactory(turbulent_bc_values.get()) );
258 libmesh_init.comm() );
260 #ifdef GRINS_USE_GRVY_TIMERS
261 grvy_timer.EndTimer(
"Initialize Solver");
264 grins.attach_grvy_timer( &grvy_timer );
271 std::tr1::shared_ptr<libMesh::EquationSystems> es = grins.get_equation_system();
277 libMesh::ExactSolution exact_sol(*es);
279 libMesh::EquationSystems es_ref( es->get_mesh() );
282 std::string solution_file = command_line(
"soln-data",
"DIE!");
283 es_ref.read( solution_file );
285 exact_sol.attach_reference_solution( &es_ref );
288 unsigned int n_vars = command_line.vector_variable_size(
"vars");
289 std::vector<std::string> vars(n_vars);
290 for(
unsigned int v = 0; v < n_vars; v++ )
292 vars[v] = command_line(
"vars",
"DIE!", v);
296 unsigned int n_norms = command_line.vector_variable_size(
"norms");
297 std::vector<std::string> norms(n_norms);
298 for(
unsigned int n = 0; n < n_norms; n++ )
300 norms[n] = command_line(
"norms",
"DIE!", n);
301 if( norms[n] != std::string(
"L2") &&
302 norms[n] != std::string(
"H1") )
304 std::cerr <<
"ERROR: Invalid norm input " << norms[n] << std::endl
305 <<
" Valid values are: L2" << std::endl
306 <<
" H1" << std::endl;
310 const std::string& system_name = grins.get_multiphysics_system_name();
313 for(
unsigned int v = 0; v < n_vars; v++ )
315 exact_sol.compute_error(system_name, vars[v]);
320 double tol = command_line(
"tol", 1.0e-10);
323 for(
unsigned int v = 0; v < n_vars; v++ )
325 for(
unsigned int n = 0; n < n_norms; n++ )
327 test_error_norm( exact_sol, system_name, vars[v], norms[n], tol, return_flag );
335 const std::string& system_name,
336 const std::string& var,
337 const std::string& norm,
348 std::cout <<
"==========================================================" << std::endl
349 <<
"Checking variable " << var <<
" using error norm " << norm <<
" with tol " << tol <<
"...";
351 if( norm == std::string(
"L2") )
353 error = exact_sol.l2_error(system_name, var);
355 else if( norm == std::string(
"H1") )
357 error = exact_sol.h1_error(system_name, var);
361 std::cerr <<
"ERROR: Invalid norm " << norm << std::endl;
369 std::cerr <<
"Tolerance exceeded for generic regression test!" << std::endl
370 <<
"tolerance = " << tol << std::endl
371 <<
"norm of error = " << error << std::endl
372 <<
"norm type = " << norm << std::endl
373 <<
"var = " << var << std::endl;
377 std::cout <<
"PASSED!" << std::endl
378 <<
"==========================================================" << std::endl;
395 cont_u.
set_func( turbulent_inlet_u );
401 cont_nu.
set_func( turbulent_inlet_nu );
403 std::multimap< GRINS::PhysicsName, GRINS::DBCContainer > mymap;
const PhysicsName incompressible_navier_stokes
Object for constructing boundary condition function objects.
Simple helper class to setup general Dirichlet boundary conditions.
libMesh::MeshFunction * turbulent_bc_values
libMesh::MeshFunction * turbulent_bc_values
const PhysicsName spalart_allmaras
std::multimap< GRINS::PhysicsName, GRINS::DBCContainer > build_dirichlet()
Builds all required libMesh::DirichletBoundary objects and adds them to DofMap.
BoundaryConditionsFactory()
void add_var_name(const GRINS::VariableName &var)
Add variables that are constrained by the Dirichlet bc.
TurbulentBdyFunctionU(libMesh::MeshFunction *_turbulent_bc_values)
libMesh::MeshFunction * turbulent_bc_values
virtual libMesh::Number operator()(const libMesh::Point &, const libMesh::Real=0)
void test_error_norm(libMesh::ExactSolution &exact_sol, const std::string &system_name, const std::string &var, const std::string &norm, const double tol, int &return_flag)
TurbulentBdyFunctionNu(libMesh::MeshFunction *_turbulent_bc_values)
void add_bc_id(const GRINS::BoundaryID bc_id)
Add boundary id's for which this Dirichlet bc is to be applied.
virtual libMesh::AutoPtr< libMesh::FunctionBase< libMesh::Number > > clone() const
void set_func(std::tr1::shared_ptr< libMesh::FunctionBase< libMesh::Number > > func)
Add the Dirichlet bc functor.
void attach_bc_factory(std::tr1::shared_ptr< BoundaryConditionsFactory > bc_factory)
virtual libMesh::AutoPtr< libMesh::FunctionBase< libMesh::Number > > clone() const
int main(int argc, char *argv[])
TurbulentBCFactory(libMesh::MeshFunction *_turbulent_bc_values)
virtual libMesh::Number operator()(const libMesh::Point &, const libMesh::Real=0)