text: // ---------------------------------------------------------------------
  //
  // Copyright (c) 2019 - 2019 by the IBAMR developers
  // All rights reserved.
  //
  // This file is part of IBAMR.
  //
  // IBAMR is free software and is distributed under the 3-clause BSD
  // license. The full text of the license can be found in the file
  // COPYRIGHT at the top level directory of IBAMR.
  //
  // ---------------------------------------------------------------------
  
  // Config files
  #include <SAMRAI_config.h>
  
  // Headers for basic PETSc functions
  #include <petscsys.h>
  
  // Headers for basic SAMRAI objects
  #include <BergerRigoutsos.h>
  #include <CartesianGridGeometry.h>
  #include <LoadBalancer.h>
  #include <StandardTagAndInitialize.h>
  
  // Headers for application-specific algorithm/data structure objects
  #include "ibamr/CFINSForcing.h"
  #include <ibamr/AdvDiffSemiImplicitHierarchyIntegrator.h>
  #include <ibamr/INSCollocatedHierarchyIntegrator.h>
  #include <ibamr/INSStaggeredHierarchyIntegrator.h>
  #include <ibamr/app_namespaces.h>
  
  #include <ibtk/AppInitializer.h>
  #include <ibtk/CartGridFunctionSet.h>
  #include <ibtk/IBTKInit.h>
  #include <ibtk/IBTK_MPI.h>
  #include <ibtk/muParserCartGridFunction.h>
  #include <ibtk/muParserRobinBcCoefs.h>
  
  #include <fstream>
  #include <iostream>
  // Function prototypes
  void output_data(Pointer<PatchHierarchy<NDIM> > patch_hierarchy,
                   Pointer<INSHierarchyIntegrator> ins_integrator,
                   const int iteration_num,
                   const double loop_time,
                   const string& data_dump_dirname);
  
  /*******************************************************************************
   * For each run, the input filename and restart information (if needed) must   *
   * be given on the command line.  For non-restarted case, command line is:     *
   *                                                                             *
   *    executable <input file name>                                             *
   *                                                                             *
   * For restarted run, command line is:                                         *
   *                                                                             *
   *    executable <input file name> <restart directory> <restart number>        *
   *                                                                             *
   *******************************************************************************/
  int
  main(int argc, char* argv[])
  {
      // Initialize IBAMR and libraries. Deinitialization is handled by this object as well.
      IBTKInit ibtk_init(argc, argv, MPI_COMM_WORLD);
  
      { // cleanup dynamically allocated objects prior to shutdown
  
          // Parse command line options, set some standard options from the input
          // file, initialize the restart database (if this is a restarted run),
          // and enable file logging.
          Pointer<AppInitializer> app_initializer = new AppInitializer(argc, argv, "INS.log");
          Pointer<Database> input_db = app_initializer->getInputDatabase();
  
          // Get various standard options set in the input file.
          const bool dump_viz_data = app_initializer->dumpVizData();
          const int viz_dump_interval = app_initializer->getVizDumpInterval();
          const bool uses_visit = dump_viz_data && app_initializer->getVisItDataWriter();
  
          const bool dump_restart_data = app_initializer->dumpRestartData();
          const int restart_dump_interval = app_initializer->getRestartDumpInterval();
          const string restart_dump_dirname = app_initializer->getRestartDumpDirectory();
  
          const bool dump_postproc_data = app_initializer->dumpPostProcessingData();
          const int postproc_data_dump_interval = app_initializer->getPostProcessingDataDumpInterval();
          const string postproc_data_dump_dirname = app_initializer->getPostProcessingDataDumpDirectory();
  
          const bool dump_timer_data = app_initializer->dumpTimerData();
          const int timer_dump_interval = app_initializer->getTimerDumpInterval();
  
          // Create major algorithm and data objects that comprise the
          // application.  These objects are configured from the input database
          // and, if this is a restarted run, from the restart database.
          Pointer<INSHierarchyIntegrator> time_integrator;
          const string solver_type =
              app_initializer->getComponentDatabase("Main")->getStringWithDefault("solver_type", "STAGGERED");
          if (solver_type == "STAGGERED")
          {
              time_integrator = new INSStaggeredHierarchyIntegrator(
                  "INSStaggeredHierarchyIntegrator",
                  app_initializer->getComponentDatabase("INSStaggeredHierarchyIntegrator"));
          }
          else if (solver_type == "COLLOCATED")
          {
              time_integrator = new INSCollocatedHierarchyIntegrator(
                  "INSCollocatedHierarchyIntegrator",
                  app_initializer->getComponentDatabase("INSCollocatedHierarchyIntegrator"));
          }
          else
          {
              TBOX_ERROR("Unsupported solver type: " << solver_type << "\n"
                                                     << "Valid options are: COLLOCATED, STAGGERED");
          }
          Pointer<AdvDiffSemiImplicitHierarchyIntegrator> adv_diff_integrator;
          adv_diff_integrator = new AdvDiffSemiImplicitHierarchyIntegrator(
              "AdvDiffSemiImplicitHierarchyIntegrator",
              app_initializer->getComponentDatabase("AdvDiffSemiImplicitHierarchyIntegrator"));
          Pointer<CartesianGridGeometry<NDIM> > grid_geometry = new CartesianGridGeometry<NDIM>(
              "CartesianGeometry", app_initializer->getComponentDatabase("CartesianGeometry"));
          Pointer<PatchHierarchy<NDIM> > patch_hierarchy = new PatchHierarchy<NDIM>("PatchHierarchy", grid_geometry);
          Pointer<StandardTagAndInitialize<NDIM> > error_detector =
              new StandardTagAndInitialize<NDIM>("StandardTagAndInitialize",
                                                 time_integrator,
                                                 app_initializer->getComponentDatabase("StandardTagAndInitialize"));
          Pointer<BergerRigoutsos<NDIM> > box_generator = new BergerRigoutsos<NDIM>();
          Pointer<LoadBalancer<NDIM> > load_balancer =
              new LoadBalancer<NDIM>("LoadBalancer", app_initializer->getComponentDatabase("LoadBalancer"));
          Pointer<GriddingAlgorithm<NDIM> > gridding_algorithm =
              new GriddingAlgorithm<NDIM>("GriddingAlgorithm",
                                          app_initializer->getComponentDatabase("GriddingAlgorithm"),
                                          error_detector,
                                          box_generator,
                                          load_balancer);
          time_integrator->registerAdvDiffHierarchyIntegrator(adv_diff_integrator);
          // Create initial condition specification objects.
          Pointer<CartGridFunction> u_init = new muParserCartGridFunction(
              "u_init", app_initializer->getComponentDatabase("VelocityInitialConditions"), grid_geometry);
          time_integrator->registerVelocityInitialConditions(u_init);
          Pointer<CartGridFunction> p_init = new muParserCartGridFunction(
              "p_init", app_initializer->getComponentDatabase("PressureInitialConditions"), grid_geometry);
          time_integrator->registerPressureInitialConditions(p_init);
  
          // Create boundary condition specification objects (when necessary).
          const IntVector<NDIM>& periodic_shift = grid_geometry->getPeriodicShift();
          vector<RobinBcCoefStrategy<NDIM>*> u_bc_coefs(NDIM);
          if (periodic_shift.min() > 0)
          {
              for (unsigned int d = 0; d < NDIM; ++d)
              {
                  u_bc_coefs[d] = NULL;
              }
          }
          else
          {
              for (unsigned int d = 0; d < NDIM; ++d)
              {
                  ostringstream bc_coefs_name_stream;
                  bc_coefs_name_stream << "u_bc_coefs_" << d;
                  const string bc_coefs_name = bc_coefs_name_stream.str();
  
                  ostringstream bc_coefs_db_name_stream;
                  bc_coefs_db_name_stream << "VelocityBcCoefs_" << d;
                  const string bc_coefs_db_name = bc_coefs_db_name_stream.str();
  
                  u_bc_coefs[d] = new muParserRobinBcCoefs(
                      bc_coefs_name, app_initializer->getComponentDatabase(bc_coefs_db_name), grid_geometry);
              }
              time_integrator->registerPhysicalBoundaryConditions(u_bc_coefs);
          }
  
          // Set up visualization plot file writers.
          Pointer<VisItDataWriter<NDIM> > visit_data_writer = app_initializer->getVisItDataWriter();
          if (uses_visit)
          {
              time_integrator->registerVisItDataWriter(visit_data_writer);
          }
          // Create body force function specification objects (when necessary).
          Pointer<CFINSForcing> polymericStressForcing;
          if (input_db->keyExists("ComplexFluid"))
          {
              polymericStressForcing = new CFINSForcing("PolymericStressForcing",
                                                        app_initializer->getComponentDatabase("ComplexFluid"),
                                                        time_integrator,
                                                        grid_geometry,
                                                        adv_diff_integrator,
                                                        visit_data_writer);
  
              time_integrator->registerBodyForceFunction(polymericStressForcing);
          }
          // Initialize hierarchy configuration and data on all patches.
          time_integrator->initializePatchHierarchy(patch_hierarchy, gridding_algorithm);
  
          // Deallocate initialization objects.
          // app_initializer.setNull();
  
          // Print the input database contents to the log file.
          plog << "Input database:\n";
          input_db->printClassData(plog);
  
          // Write out initial visualization data.
          int iteration_num = time_integrator->getIntegratorStep();
          double loop_time = time_integrator->getIntegratorTime();
          if (dump_viz_data && uses_visit)
          {
              pout << "\n\nWriting visualization files...\n\n";
              time_integrator->setupPlotData();
              visit_data_writer->writePlotData(patch_hierarchy, iteration_num, loop_time);
          }
  
          // Main time step loop.
          double loop_time_end = time_integrator->getEndTime();
          double dt = 0.0;
          while (!MathUtilities<double>::equalEps(loop_time, loop_time_end) && time_integrator->stepsRemaining())
          {
              iteration_num = time_integrator->getIntegratorStep();
              loop_time = time_integrator->getIntegratorTime();
  
              pout << "\n";
              pout << "+++++++++++++++++++++++++++++++++++++++++++++++++++\n";
              pout << "At beginning of timestep # " << iteration_num << "\n";
              pout << "Simulation time is " << loop_time << "\n";
  
              dt = time_integrator->getMaximumTimeStepSize();
              time_integrator->advanceHierarchy(dt);
              loop_time += dt;
  
              pout << "\n";
              pout << "At end       of timestep # " << iteration_num << "\n";
              pout << "Simulation time is " << loop_time << "\n";
              pout << "+++++++++++++++++++++++++++++++++++++++++++++++++++\n";
              pout << "\n";
  
              // At specified intervals, write visualization and restart files,
              // print out timer data, and store hierarchy data for post
              // processing.
              iteration_num += 1;
              const bool last_step = !time_integrator->stepsRemaining();
              if (dump_viz_data && uses_visit && (iteration_num % viz_dump_interval == 0 || last_step))
              {
                  pout << "\nWriting visualization files...\n\n";
                  time_integrator->setupPlotData();
                  visit_data_writer->writePlotData(patch_hierarchy, iteration_num, loop_time);
              }
              if (dump_restart_data && (iteration_num % restart_dump_interval == 0 || last_step))
              {
                  pout << "\nWriting restart files...\n\n";
                  RestartManager::getManager()->writeRestartFile(restart_dump_dirname, iteration_num);
              }
              if (dump_timer_data && (iteration_num % timer_dump_interval == 0 || last_step))
              {
                  pout << "\nWriting timer data...\n\n";
                  TimerManager::getManager()->print(plog);
              }
              if (dump_postproc_data && (iteration_num % postproc_data_dump_interval == 0 || last_step))
              {
                  output_data(patch_hierarchy, time_integrator, iteration_num, loop_time, postproc_data_dump_dirname);
              }
          }
  
      } // cleanup dynamically allocated objects prior to shutdown
  } // main
  
  void
  output_data(Pointer<PatchHierarchy<NDIM> > patch_hierarchy,
              Pointer<INSHierarchyIntegrator> ins_integrator,
              const int iteration_num,
              const double loop_time,
              const string& data_dump_dirname)
  {
      plog << "writing hierarchy data at iteration " << iteration_num << " to disk" << endl;
      plog << "simulation time is " << loop_time << endl;
      string file_name = data_dump_dirname + "/" + "hier_data.";
      char temp_buf[128];
      sprintf(temp_buf, "%05d.samrai.%05d", iteration_num, IBTK_MPI::getRank());
      file_name += temp_buf;
      Pointer<HDFDatabase> hier_db = new HDFDatabase("hier_db");
      hier_db->create(file_name);
      VariableDatabase<NDIM>* var_db = VariableDatabase<NDIM>::getDatabase();
      ComponentSelector hier_data;
      hier_data.setFlag(var_db->mapVariableAndContextToIndex(ins_integrator->getVelocityVariable(),
                                                             ins_integrator->getCurrentContext()));
      hier_data.setFlag(var_db->mapVariableAndContextToIndex(ins_integrator->getPressureVariable(),
                                                             ins_integrator->getCurrentContext()));
      patch_hierarchy->putToDatabase(hier_db->putDatabase("PatchHierarchy"), hier_data);
      hier_db->putDouble("loop_time", loop_time);
      hier_db->putInteger("iteration_num", iteration_num);
      hier_db->close();
      return;
  } // output_data
  