Output after every AMR step

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);
 

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);
 

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;

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;

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 );

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;

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 ),

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<libMesh::AdjointRefinementEstimator*>( context.error_estimator.get() );
-              std::cout<<"The error estimate for QoI("<<i<<") is: "<<adjoint_ref_error_estimator->get_global_QoI_error_estimate(i)<<std::endl;
-            }
+        if (context.output_every_amr)
+          context.vis->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<libMesh::AdjointRefinementEstimator*>( context.error_estimator.get() );
+          std::cout<<"The error estimate for QoI("<<i<<") is: "<<adjoint_ref_error_estimator->get_global_QoI_error_estimate(i)<<std::endl;
+        }
+
+    // Check for convergence of error
+    bool converged = this->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

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" <<std::endl;
+
     context.system->deltat = this->_deltat;
 
     libMesh::Real sim_time;

src/visualization/include/grins/visualization.h

@@ -58,9 +58,13 @@ namespace GRINS
     virtual ~Visualization();
 
     void output( std::shared_ptr<libMesh::EquationSystems> equation_system );
+
     void output( std::shared_ptr<libMesh::EquationSystems> equation_system,
                  const unsigned int time_step, const libMesh::Real time );
 
+    void output_amr( std::shared_ptr<libMesh::EquationSystems> equation_system,
+                 const unsigned int amr_step );
+
     void output_residual( std::shared_ptr<libMesh::EquationSystems> equation_system,
                           GRINS::MultiphysicsSystem* system );
 

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<libMesh::EquationSystems> 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<libMesh::EquationSystems> 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<libMesh::EquationSystems> 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

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<libMesh::UnstructuredMesh> 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<GRINS::RayfireMesh> 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 );