diff --git a/src/qoi/src/qoi_output.C b/src/qoi/src/qoi_output.C index 15c5791b8..997bf1784 100644 --- a/src/qoi/src/qoi_output.C +++ b/src/qoi/src/qoi_output.C @@ -71,7 +71,7 @@ namespace GRINS if( comm.rank() == 0 ) { std::ofstream output; - output.open( _file_prefix+".dat" ); + output.open( _file_prefix+".dat",std::ofstream::app ); qois.output_qoi(output); diff --git a/src/qoi/src/rayfire_mesh.C b/src/qoi/src/rayfire_mesh.C index f6946da54..3e344f80e 100644 --- a/src/qoi/src/rayfire_mesh.C +++ b/src/qoi/src/rayfire_mesh.C @@ -887,7 +887,8 @@ namespace GRINS J_inv.vector_mult(delta,F); // check for convergence - if ( delta.l2_norm() < libMesh::TOLERANCE ) + libMesh::Real tol = std::min( libMesh::TOLERANCE, libMesh::TOLERANCE*edge_elem->hmax() ); + if ( delta.l2_norm() < tol ) { libMesh::Point intersect(X,Y); @@ -1136,7 +1137,8 @@ namespace GRINS J.lu_solve(F,delta); // check for convergence - if ( delta.l2_norm() < libMesh::TOLERANCE ) + libMesh::Real tol = std::min( libMesh::TOLERANCE, libMesh::TOLERANCE*side_elem->hmax() ); + if ( delta.l2_norm() < tol ) { libMesh::Point intersect(X,Y,Z); diff --git a/src/solver/include/grins/simulation.h b/src/solver/include/grins/simulation.h index c232c662d..926de8ada 100644 --- a/src/solver/include/grins/simulation.h +++ b/src/solver/include/grins/simulation.h @@ -140,6 +140,7 @@ namespace GRINS // Visualization options bool _output_vis; + bool _output_every_amr; bool _output_adjoint; bool _output_residual; bool _output_residual_sensitivities; diff --git a/src/solver/include/grins/solver_context.h b/src/solver/include/grins/solver_context.h index 1bcfb8598..dd8980fbf 100644 --- a/src/solver/include/grins/solver_context.h +++ b/src/solver/include/grins/solver_context.h @@ -61,6 +61,7 @@ namespace GRINS unsigned int timesteps_per_vis; unsigned int timesteps_per_perflog; bool output_vis; + bool output_every_amr; bool output_adjoint; bool output_residual; bool output_residual_sensitivities; diff --git a/src/solver/include/grins/steady_mesh_adaptive_solver.h b/src/solver/include/grins/steady_mesh_adaptive_solver.h index daa4195a1..3640551f6 100644 --- a/src/solver/include/grins/steady_mesh_adaptive_solver.h +++ b/src/solver/include/grins/steady_mesh_adaptive_solver.h @@ -60,6 +60,8 @@ namespace GRINS protected: + bool amr_helper(SolverContext & context, unsigned int r_step, bool do_amr); + virtual void init_time_solver( MultiphysicsSystem* system ); void check_qoi_error_option_consistency( SolverContext& context ); diff --git a/src/solver/src/simulation.C b/src/solver/src/simulation.C index a950b2fba..264ea8c39 100644 --- a/src/solver/src/simulation.C +++ b/src/solver/src/simulation.C @@ -65,6 +65,7 @@ namespace GRINS _print_scalars( input("screen-options/print_scalars", false ) ), _qoi_output( new QoIOutput(input) ), _output_vis( input("vis-options/output_vis", false ) ), + _output_every_amr( input("vis-options/output_every_amr", false ) ), _output_adjoint( input("vis-options/output_adjoint", false ) ), _output_residual( input( "vis-options/output_residual", false ) ), _output_residual_sensitivities( input( "vis-options/output_residual_sensitivities", false ) ), @@ -230,6 +231,7 @@ namespace GRINS context.timesteps_per_vis = _timesteps_per_vis; context.timesteps_per_perflog = _timesteps_per_perflog; context.output_vis = _output_vis; + context.output_every_amr = _output_every_amr; context.output_residual = _output_residual; context.output_residual_sensitivities = _output_residual_sensitivities; context.output_solution_sensitivities = _output_solution_sensitivities; diff --git a/src/solver/src/solver_context.C b/src/solver/src/solver_context.C index fe89b08c1..7966663cf 100644 --- a/src/solver/src/solver_context.C +++ b/src/solver/src/solver_context.C @@ -38,6 +38,7 @@ namespace GRINS timesteps_per_vis( 1 ), timesteps_per_perflog( 1 ), output_vis( false ), + output_every_amr( false ), output_adjoint(false), output_residual( false ), output_residual_sensitivities( false ), diff --git a/src/solver/src/steady_mesh_adaptive_solver.C b/src/solver/src/steady_mesh_adaptive_solver.C index 6a44a3d14..d9cb0c203 100644 --- a/src/solver/src/steady_mesh_adaptive_solver.C +++ b/src/solver/src/steady_mesh_adaptive_solver.C @@ -86,6 +86,8 @@ namespace GRINS libMesh::MeshBase& mesh = context.equation_system->get_mesh(); this->build_mesh_refinement( mesh ); + bool converged = false; + /*! \todo This output cannot be toggled in the input file, but it should be able to be. */ std::cout << "==========================================================" << std::endl << "Performing " << _mesh_adaptivity_options.max_refinement_steps() @@ -94,71 +96,92 @@ namespace GRINS for ( unsigned int r_step = 0; r_step < _mesh_adaptivity_options.max_refinement_steps(); r_step++ ) { - std::cout << "==========================================================" << std::endl + std::cout << std::endl + << "==========================================================" << std::endl << "Adaptive Refinement Step " << r_step << std::endl << "==========================================================" << std::endl; - // Solve the forward problem - context.system->solve(); + converged = this->amr_helper(context,r_step,true); + if (converged) + break; - if( context.output_vis ) - { - context.postprocessing->update_quantities( *(context.equation_system) ); - context.vis->output( context.equation_system ); - } + } // r_step for-loop - // Solve adjoint system - if(context.do_adjoint_solve) - this->steady_adjoint_solve(context); + // this will print output the error estimates and/or meshes after the final AMR run + if (! converged) + { + std::cout << std::endl + << "==========================================================" << std::endl + << "Post-AMR Output " << std::endl + << "==========================================================" << std::endl; + this->amr_helper(context,_mesh_adaptivity_options.max_refinement_steps(),false); + } - if(context.output_adjoint) - context.vis->output_adjoint(context.equation_system, context.system); + return; + } - if( context.output_residual ) - { - context.vis->output_residual( context.equation_system, context.system ); - } + bool SteadyMeshAdaptiveSolver::amr_helper(SolverContext & context, unsigned int r_step, bool do_amr) + { + // Solve the forward problem + context.system->solve(); - // Now we construct the data structures for the mesh refinement process - libMesh::ErrorVector error; - this->estimate_error_for_amr( context, error ); + if( context.output_vis ) + { + context.postprocessing->update_quantities( *(context.equation_system) ); - // Get the global error estimate if you can and are asked to - if( _error_estimator_options.compute_qoi_error_estimate() ) - for(unsigned int i = 0; i != context.system->qoi.size(); i++) - { - libMesh::AdjointRefinementEstimator* adjoint_ref_error_estimator = libMesh::cast_ptr( context.error_estimator.get() ); - std::cout<<"The error estimate for QoI("<get_global_QoI_error_estimate(i)<output_amr( context.equation_system,r_step ); + else + context.vis->output( context.equation_system ); + } - // Check for convergence of error - bool converged = this->check_for_convergence( context, error ); + // Solve adjoint system + if(context.do_adjoint_solve) + this->steady_adjoint_solve(context); - if( converged ) - { - // Break out of adaptive loop - std::cout << "==========================================================" << std::endl - << "Convergence detected!" << std::endl - << "==========================================================" << std::endl; - break; - } - else + if(context.output_adjoint) + context.vis->output_adjoint(context.equation_system, context.system); + + if( context.output_residual ) + { + context.vis->output_residual( context.equation_system, context.system ); + } + + // Now we construct the data structures for the mesh refinement process + libMesh::ErrorVector error; + this->estimate_error_for_amr( context, error ); + + // Get the global error estimate if you can and are asked to + if( _error_estimator_options.compute_qoi_error_estimate() ) + for(unsigned int i = 0; i != context.system->qoi.size(); i++) + { + libMesh::AdjointRefinementEstimator* adjoint_ref_error_estimator = libMesh::cast_ptr( context.error_estimator.get() ); + std::cout<<"The error estimate for QoI("<get_global_QoI_error_estimate(i)<check_for_convergence(context,error); + + if( converged ) + { + // Break out of adaptive loop + std::cout << "==========================================================" << std::endl + << "Convergence detected!" << std::endl + << "==========================================================" << std::endl; + } + else + { + if (do_amr) { - // Only bother refining if we're on the last step. - if( r_step < _mesh_adaptivity_options.max_refinement_steps() -1 ) - { - this->perform_amr(context, error); - - // It's helpful to print the qoi along the way, but only do it if the user - // asks for it - if( context.qoi_output->output_qoi_set() ) - this->print_qoi(context); - } + // print the QoI before AMR + if( context.qoi_output->output_qoi_set() ) + this->print_qoi(context); + + this->perform_amr(context, error); } + } - } // r_step for-loop - - return; + return converged; } void SteadyMeshAdaptiveSolver::adjoint_qoi_parameter_sensitivity diff --git a/src/solver/src/unsteady_mesh_adaptive_solver.C b/src/solver/src/unsteady_mesh_adaptive_solver.C index 0f93ceecc..3b47e3f9f 100644 --- a/src/solver/src/unsteady_mesh_adaptive_solver.C +++ b/src/solver/src/unsteady_mesh_adaptive_solver.C @@ -44,6 +44,11 @@ namespace GRINS void UnsteadyMeshAdaptiveSolver::solve( SolverContext& context ) { + if (context.output_every_amr) + std::cout << std::endl + << "WARNING: 'output_every_amr' flag is not yet implemented" + << "in UnsteadyMeshAdaptiveSolver" <deltat = this->_deltat; libMesh::Real sim_time; diff --git a/src/visualization/include/grins/visualization.h b/src/visualization/include/grins/visualization.h index a1bc62302..941b6b308 100644 --- a/src/visualization/include/grins/visualization.h +++ b/src/visualization/include/grins/visualization.h @@ -58,9 +58,13 @@ namespace GRINS virtual ~Visualization(); void output( std::shared_ptr equation_system ); + void output( std::shared_ptr equation_system, const unsigned int time_step, const libMesh::Real time ); + void output_amr( std::shared_ptr equation_system, + const unsigned int amr_step ); + void output_residual( std::shared_ptr equation_system, GRINS::MultiphysicsSystem* system ); diff --git a/src/visualization/src/visualization.C b/src/visualization/src/visualization.C index 3fb1a9952..0a17574af 100644 --- a/src/visualization/src/visualization.C +++ b/src/visualization/src/visualization.C @@ -87,8 +87,6 @@ namespace GRINS { _output_format.push_back( input("vis-options/output_format", "DIE", i ) ); } - - return; } Visualization::~Visualization() @@ -99,8 +97,6 @@ namespace GRINS void Visualization::output( std::shared_ptr equation_system ) { this->dump_visualization( equation_system, _vis_output_file_prefix, 0.0 ); - - return; } void Visualization::output @@ -116,15 +112,26 @@ namespace GRINS filename+="."+suffix.str(); this->dump_visualization( equation_system, filename, time ); + } - return; + void Visualization::output_amr + ( std::shared_ptr equation_system, + const unsigned int amr_step ) + { + std::stringstream suffix; + + suffix << "amr" << amr_step; + + std::string filename = this->_vis_output_file_prefix; + filename+="."+suffix.str(); + + this->dump_visualization( equation_system, filename, 0.0 ); } void Visualization::output_residual( std::shared_ptr equation_system, MultiphysicsSystem* system ) { this->output_residual( equation_system, system, 0, 0.0 ); - return; } void Visualization::output_residual_sensitivities @@ -249,8 +256,6 @@ namespace GRINS libmesh_error(); } } // End loop over formats - - return; } } // namespace GRINS diff --git a/test/unit/rayfire_test_3d.C b/test/unit/rayfire_test_3d.C index 80b20e239..6abef86e4 100644 --- a/test/unit/rayfire_test_3d.C +++ b/test/unit/rayfire_test_3d.C @@ -72,6 +72,7 @@ namespace GRINSTesting CPPUNIT_TEST( hex_27elem_3x3x3 ); CPPUNIT_TEST( fire_through_vertex ); CPPUNIT_TEST( origin_between_elems ); + CPPUNIT_TEST( hex27_from_larger_mesh ); CPPUNIT_TEST_SUITE_END(); @@ -222,6 +223,67 @@ namespace GRINSTesting this->run_test_with_mesh_from_file(origin,pts,exit_ids,"mesh_hex27_27elem.in"); } + //! This test replicates a failure from a larger mesh, hence the seemingly arbitrary coordinates + void hex27_from_larger_mesh() + { + libMesh::Point origin = libMesh::Point(0.0030875001992899573, 0.069253246753246747, -6.7698557104360308e-11); + libMesh::Real theta = 1.57079632679; + libMesh::Real phi = 1.57079632679; + + std::shared_ptr mesh( new libMesh::SerialMesh(*TestCommWorld) ); + + mesh->set_mesh_dimension(3); + + mesh->add_point( libMesh::Point(0.001122532015133644, 0.070792207792207781, 0.0011225320151336397),0 ); + mesh->add_point( libMesh::Point(0.0015874999999999999, 0.070792207792207795, -4.1403598545004364e-18),1 ); + mesh->add_point( libMesh::Point(0.0031562499999999993, 0.070792207792207781, -3.9482433878841953e-18),2 ); + mesh->add_point( libMesh::Point(0.0029363425824496629, 0.070792207792207795, 0.0013735659049402692),3 ); + + mesh->add_point( libMesh::Point(0.001122532015133644, 0.069253246753246747, 0.0011225320151336397),4 ); + mesh->add_point( libMesh::Point(0.0015874999999999999, 0.069253246753246747, -4.0461256689816297e-18),5 ); + mesh->add_point( libMesh::Point(0.0031562499999999985, 0.069253246753246733, -3.8540092023653886e-18),6 ); + mesh->add_point( libMesh::Point(0.0029363425824496621, 0.069253246753246719, 0.0013735659049402692),7 ); + + mesh->add_point( libMesh::Point(0.0014666587578616678, 0.070792207792207795, 0.00060750994887957582),8 ); + mesh->add_point( libMesh::Point(0.0023718749999999999, 0.070792207792207795, -4.0443016211923158e-18),9 ); + mesh->add_point( libMesh::Point(0.0030462962912248311, 0.070792207792207795, 0.00068678295247013264),10 ); + mesh->add_point( libMesh::Point(0.0020294372987916536, 0.070792207792207795, 0.0012480489600369543),11 ); + + mesh->add_point( libMesh::Point(0.0011225320151336442, 0.070022727272727264, 0.0011225320151336399),12 ); + mesh->add_point( libMesh::Point(0.0015874999999999999, 0.070022727272727264, -4.0932427617410323e-18),13 ); + mesh->add_point( libMesh::Point(0.0031562499999999989, 0.070022727272727264, -3.9011262951247919e-18),14 ); + mesh->add_point( libMesh::Point(0.0029363425824496625, 0.070022727272727264, 0.0013735659049402692),15 ); + + mesh->add_point( libMesh::Point(0.0014666587578616678, 0.069253246753246747, 0.00060750994887957582),16 ); + mesh->add_point( libMesh::Point(0.002371874999999999, 0.069253246753246733, -3.9500674356735084e-18),17 ); + mesh->add_point( libMesh::Point(0.0030462962912248303, 0.069253246753246733, 0.00068678295247013264),18 ); + mesh->add_point( libMesh::Point(0.0020294372987916531, 0.069253246753246733, 0.0012480489600369543),19 ); + + mesh->add_point( libMesh::Point(0.0022564775245432493, 0.070792207792207795, 0.00064714645067485423),20 ); + mesh->add_point( libMesh::Point(0.0014666587578616675, 0.070022727272727264, 0.00060750994887957593),21 ); + mesh->add_point( libMesh::Point(0.0023718749999999994, 0.070022727272727264, -3.9971845284329117e-18),22 ); + mesh->add_point( libMesh::Point(0.0030462962912248307, 0.070022727272727264, 0.00068678295247013264),23 ); + mesh->add_point( libMesh::Point(0.0020294372987916531, 0.070022727272727264, 0.0012480489600369543),24 ); + mesh->add_point( libMesh::Point(0.0022564775245432489, 0.069253246753246733, 0.00064714645067485434),25 ); + + mesh->add_point( libMesh::Point(0.0022564775245432489, 0.070022727272727264, 0.00064714645067485413),26 ); + + libMesh::Elem* elem = mesh->add_elem( new libMesh::Hex27 ); + for (unsigned int n=0; n<27; n++) + elem->set_node(n) = mesh->node_ptr(n); + + mesh->prepare_for_use(); + + std::shared_ptr rayfire( new GRINS::RayfireMesh(origin,theta,phi) ); + + rayfire->init(*mesh); + + // make sure we get a rayfire elem + const libMesh::Elem * rayfire_elem = rayfire->map_to_rayfire_elem(0); + + CPPUNIT_ASSERT(rayfire_elem); + } + }; CPPUNIT_TEST_SUITE_REGISTRATION( RayfireTest3D );