26 #include "grins_config.h"
39 #include "libmesh/dirichlet_boundaries.h"
40 #include "libmesh/dof_map.h"
41 #include "libmesh/fe_base.h"
42 #include "libmesh/fe_interface.h"
43 #include "libmesh/mesh_function.h"
44 #include "libmesh/serial_mesh.h"
45 #include "libmesh/edge_edge2.h"
46 #include "libmesh/mesh_generation.h"
47 #include "libmesh/linear_implicit_system.h"
48 #include "libmesh/gmv_io.h"
49 #include "libmesh/exact_solution.h"
50 #include "libmesh/zero_function.h"
51 #include "libmesh/composite_function.h"
52 #include "libmesh/zero_function.h"
55 #ifdef GRINS_HAVE_GRVY
61 class MultiphysicsSystem;
65 const std::string& system_name,
66 const std::string& var,
67 const std::string& norm,
79 { this->_initialized =
true; }
81 virtual libMesh::Number
operator() (
const libMesh::Point&,
const libMesh::Real = 0)
82 { libmesh_not_implemented(); }
84 virtual libMesh::Number
compute_final_value(
const libMesh::DenseVector<libMesh::Number>& u_nu_values ) =0;
87 const libMesh::Real t,
88 libMesh::DenseVector<libMesh::Number>& output)
94 libMesh::Point p_copy(p);
96 p_copy(0) = p_copy(1);
102 p_copy(0) = 1 - p_copy(0);
104 p_copy(0) = 2*p_copy(0);
106 libMesh::DenseVector<libMesh::Number> u_nu_values;
127 {
return u_nu_values(0)/21.995539; }
129 virtual libMesh::UniquePtr<libMesh::FunctionBase<libMesh::Number> >
clone()
const
142 {
return u_nu_values(1)/(2.0*21.995539); }
144 virtual libMesh::UniquePtr<libMesh::FunctionBase<libMesh::Number> >
clone()
const
176 virtual libMesh::UniquePtr<libMesh::FunctionBase<libMesh::Number> >
179 std::vector<std::string>& var_names,
182 libmesh_assert_equal_to(var_names.size(), 2);
183 libmesh_assert_equal_to(var_names[0], std::string(
"u"));
184 libmesh_assert_equal_to(var_names[1], std::string(
"v"));
186 libMesh::UniquePtr<libMesh::CompositeFunction<libMesh::Number> >
189 libMesh::UniquePtr<libMesh::FunctionBase<libMesh::Number> >
192 std::vector<GRINS::VariableIndex> vars(1);
193 vars[0] = system.variable_number(var_names[0]);
194 composite_func->attach_subfunction(*bound_func, vars);
196 vars[0] = system.variable_number(var_names[1]);
197 composite_func->attach_subfunction(libMesh::ZeroFunction<libMesh::Number>(), vars);
199 return libMesh::UniquePtr<libMesh::FunctionBase<libMesh::Number> >( composite_func.release() );
213 virtual libMesh::UniquePtr<libMesh::FunctionBase<libMesh::Number> >
216 std::vector<std::string>& var_names,
219 libmesh_assert_equal_to(var_names.size(), 1);
220 libmesh_assert_equal_to(var_names[0], std::string(
"nu"));
222 libMesh::UniquePtr<libMesh::CompositeFunction<libMesh::Number> >
225 libMesh::UniquePtr<libMesh::FunctionBase<libMesh::Number> >
228 std::vector<GRINS::VariableIndex> vars(1);
229 vars[0] = system.variable_number(var_names[0]);
230 composite_func->attach_subfunction(*bound_func, vars);
232 return libMesh::UniquePtr<libMesh::FunctionBase<libMesh::Number> >( composite_func.release() );
241 int main(
int argc,
char* argv[])
244 #ifdef GRINS_USE_GRVY_TIMERS
245 GRVY::GRVY_Timer_Class grvy_timer;
246 grvy_timer.Init(
"GRINS Timer");
253 std::cerr <<
"Error: Must specify libMesh input file." << std::endl;
258 std::string libMesh_input_filename = argv[1];
261 GetPot libMesh_inputfile( libMesh_input_filename );
263 #ifdef GRINS_USE_GRVY_TIMERS
264 grvy_timer.BeginTimer(
"Initialize Solver");
268 libMesh::LibMeshInit libmesh_init(argc, argv);
271 libMesh::SerialMesh mesh(libmesh_init.comm());
273 GetPot command_line(argc,argv);
275 std::string oned_mesh = command_line(
"mesh-1d",
"DIE!");
276 mesh.read(oned_mesh);
279 libMesh::EquationSystems equation_systems (mesh);
281 std::string oned_data = command_line(
"data-1d",
"DIE!");
282 equation_systems.read(oned_data, libMesh::XdrMODE::READ,
283 libMesh::EquationSystems::READ_HEADER |
284 libMesh::EquationSystems::READ_DATA |
285 libMesh::EquationSystems::READ_ADDITIONAL_DATA);
288 libMesh::System & turbulent_bc_system = equation_systems.get_system<libMesh::System>(
"flow");
292 turbulent_bc_system.update();
295 equation_systems.print_info();
298 libMesh::UniquePtr<libMesh::MeshFunction> turbulent_bc_values;
300 libMesh::UniquePtr<libMesh::NumericVector<libMesh::Number> > turbulent_bc_soln = libMesh::NumericVector<libMesh::Number>::build(equation_systems.comm());
302 std::vector<libMesh::Number> flow_soln;
304 turbulent_bc_system.update_global_solution(flow_soln);
306 turbulent_bc_soln->init(turbulent_bc_system.solution->size(),
true, libMesh::SERIAL);
308 (*turbulent_bc_soln) = flow_soln;
310 std::vector<unsigned int>turbulent_bc_system_variables;
311 turbulent_bc_system_variables.push_back(0);
312 turbulent_bc_system_variables.push_back(1);
314 turbulent_bc_values = libMesh::UniquePtr<libMesh::MeshFunction>
315 (
new libMesh::MeshFunction(equation_systems,
317 turbulent_bc_system.get_dof_map(),
318 turbulent_bc_system_variables ));
320 turbulent_bc_values->init();
328 libmesh_init.comm() );
330 #ifdef GRINS_USE_GRVY_TIMERS
331 grvy_timer.EndTimer(
"Initialize Solver");
334 grins.attach_grvy_timer( &grvy_timer );
341 GRINS::SharedPtr<libMesh::EquationSystems> es = grins.get_equation_system();
347 libMesh::ExactSolution exact_sol(*es);
349 libMesh::EquationSystems es_ref( es->get_mesh() );
352 std::string solution_file = command_line(
"soln-data",
"DIE!");
353 es_ref.read( solution_file, libMesh::XdrMODE::DECODE,
354 libMesh::EquationSystems::READ_HEADER |
355 libMesh::EquationSystems::READ_DATA |
356 libMesh::EquationSystems::READ_ADDITIONAL_DATA);
358 exact_sol.attach_reference_solution( &es_ref );
361 unsigned int n_vars = command_line.vector_variable_size(
"vars");
362 std::vector<std::string> vars(n_vars);
363 for(
unsigned int v = 0; v < n_vars; v++ )
365 vars[v] = command_line(
"vars",
"DIE!", v);
369 unsigned int n_norms = command_line.vector_variable_size(
"norms");
370 std::vector<std::string> norms(n_norms);
371 for(
unsigned int n = 0; n < n_norms; n++ )
373 norms[n] = command_line(
"norms",
"DIE!", n);
374 if( norms[n] != std::string(
"L2") &&
375 norms[n] != std::string(
"H1") )
377 std::cerr <<
"ERROR: Invalid norm input " << norms[n] << std::endl
378 <<
" Valid values are: L2" << std::endl
379 <<
" H1" << std::endl;
383 const std::string& system_name = grins.get_multiphysics_system_name();
386 for(
unsigned int v = 0; v < n_vars; v++ )
388 exact_sol.compute_error(system_name, vars[v]);
393 double tol = command_line(
"tol", 1.0e-10);
396 for(
unsigned int v = 0; v < n_vars; v++ )
398 for(
unsigned int n = 0; n < n_norms; n++ )
400 test_error_norm( exact_sol, system_name, vars[v], norms[n], tol, return_flag );
408 const std::string& system_name,
409 const std::string& var,
410 const std::string& norm,
421 std::cout <<
"==========================================================" << std::endl
422 <<
"Checking variable " << var <<
" using error norm " << norm <<
" with tol " << tol <<
"...";
424 if( norm == std::string(
"L2") )
426 error = exact_sol.l2_error(system_name, var);
428 else if( norm == std::string(
"H1") )
430 error = exact_sol.h1_error(system_name, var);
434 std::cerr <<
"ERROR: Invalid norm " << norm << std::endl;
442 std::cerr <<
"Tolerance exceeded for generic regression test!" << std::endl
443 <<
"tolerance = " << tol << std::endl
444 <<
"norm of error = " << error << std::endl
445 <<
"norm type = " << norm << std::endl
446 <<
"var = " << var << std::endl;
450 std::cout <<
"PASSED!" << std::endl
451 <<
"==========================================================" << std::endl;
SATurbBCFactoryBase(const std::string &bc_type_name)
SATurbUBCFactory grins_factory_testing_turb_u_bc("testing_turb_u")
TurbBoundFuncBase(libMesh::MeshFunction *turbulent_bc_values)
virtual libMesh::UniquePtr< libMesh::FunctionBase< libMesh::Number > > build_func(const GetPot &, GRINS::MultiphysicsSystem &system, std::vector< std::string > &var_names, const std::string &)
Builds the FunctionBase object for boundary condition.
SATurbNuBCFactory grins_factory_testing_turb_nu_bc("testing_turb_nu")
virtual libMesh::UniquePtr< libMesh::FunctionBase< libMesh::Number > > clone() const
libMesh::MeshFunction * _turbulent_bc_values
virtual libMesh::UniquePtr< libMesh::FunctionBase< libMesh::Number > > build_func(const GetPot &, GRINS::MultiphysicsSystem &system, std::vector< std::string > &var_names, const std::string &)
Builds the FunctionBase object for boundary condition.
virtual libMesh::Number operator()(const libMesh::Point &, const libMesh::Real=0)
virtual libMesh::UniquePtr< libMesh::FunctionBase< libMesh::Number > > clone() const
virtual libMesh::Number compute_final_value(const libMesh::DenseVector< libMesh::Number > &u_nu_values)
static libMesh::MeshFunction * _turbulent_bc_values
SATurbNuBCFactory(const std::string &bc_type_name)
TurbulentBdyFunctionU(libMesh::MeshFunction *turbulent_bc_values)
DirichletBCFactoryFunctionBase(const std::string &bc_type_name)
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)
virtual libMesh::Number compute_final_value(const libMesh::DenseVector< libMesh::Number > &u_nu_values)=0
Interface with libMesh for solving Multiphysics problems.
static void set_turb_bc_values(libMesh::MeshFunction *turbulent_bc_values)
TurbulentBdyFunctionNu(libMesh::MeshFunction *turbulent_bc_values)
SATurbUBCFactory(const std::string &bc_type_name)
int main(int argc, char *argv[])
virtual libMesh::Number compute_final_value(const libMesh::DenseVector< libMesh::Number > &u_nu_values)