<commit_before>// Copyright (c) 2002-2010, <NAME>
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
//    * Redistributions of source code must retain the above copyright notice,
//      this list of conditions and the following disclaimer.
//
//    * Redistributions in binary form must reproduce the above copyright
//      notice, this list of conditions and the following disclaimer in the
//      documentation and/or other materials provided with the distribution.
//
//    * Neither the name of New York University nor the names of its
//      contributors may be used to endorse or promote products derived from
//      this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.

// Config files
#include <IBAMR_config.h>
#include <SAMRAI_config.h>

// Headers for basic PETSc functions
#include <petsc.h>

// Headers for basic libMesh objects
#include <exodusII_io.h>
#include <mesh.h>
#include <mesh_generation.h>
#include <periodic_boundaries.h>
#include <quadrature.h>
#include <string_to_enum.h>

// Headers for application-specific algorithm/data structure objects
#include <ibamr/IBFEHierarchyIntegrator.h>
#include <ibtk/muParserCartGridFunction.h>
#include <ibtk/muParserRobinBcCoefs.h>
#include <BergerRigoutsos.h>
#include <StandardTagAndInitialize.h>

// C++ namespace delcarations
#include <ibamr/namespaces.h>
using namespace IBAMR;
using namespace IBTK;
using namespace libMesh;
using namespace std;

// Elasticity model data.
namespace ModelData
{
// Problem parameters.
static const double R = 0.25;
static const double w = 0.0625;
static const double gamma = 0.0;
static const double mu = 1.0;

// Coordinate mapping function.
void
coordinate_mapping_function(
    Point& X,
    const Point& s,
    void* /*ctx*/)
{
    const double X0 = (R+      s(1))*cos(s(0)/R);
    const double X1 = (R+gamma+s(1))*sin(s(0)/R);

    X(0) = cos(22.5*M_PI/180.0)*X0 - sin(22.5*M_PI/180.0)*X1 + 0.5;
    X(1) = sin(22.5*M_PI/180.0)*X0 + cos(22.5*M_PI/180.0)*X1 + 0.5;
    return;
}// coordinate_mapping_function

// Stress tensor function.
bool smooth_case = false;
void
PK1_stress_function(
    TensorValue<double>& PP,
    const TensorValue<double>& dX_ds,
    const Point& /*X*/,
    const Point& /*s*/,
    Elem* const /*elem*/,
    NumericVector<double>& /*X_vec*/,
    const std::vector<NumericVector<double>*>& /*system_data*/,
    const double& /*time*/,
    void* /*ctx*/)
{
    PP = (mu/w)*dX_ds;
    if (smooth_case)
    {
        PP(0,1) = 0.0;
        PP(1,1) = 0.0;
    }
    return;
}// PK1_stress_function
}
using namespace ModelData;

namespace Postprocessing
{
void
dump_hier_data(
    const int iteration_num,
    const double loop_time,
    const double dt,
    const string hier_dump_dirname,
    const int hier_dump_interval,
    Pointer<PatchHierarchy<NDIM> > patch_hierarchy,
    const ComponentSelector& hier_data,
    Mesh& mesh,
    EquationSystems& equation_systems)
{
    string file_name;

    // Write Cartesian data.
    file_name = hier_dump_dirname + "/" + "hier_data.";
    char temp_buf[128];
    sprintf(temp_buf, "%05d.samrai.%05d", iteration_num, SAMRAI_MPI::getRank());
    file_name += temp_buf;

    Pointer<HDFDatabase> hier_db = new HDFDatabase("hier_db");
    hier_db->create(file_name);
    patch_hierarchy->putToDatabase(hier_db->putDatabase("PatchHierarchy"), hier_data);
    hier_db->putDouble("loop_time", loop_time);
    hier_db->putDouble("dt", dt);
    hier_db->putInteger("iteration_num", iteration_num);
    hier_db->putInteger("hier_dump_interval", hier_dump_interval);
    hier_db->close();

    // Write Lagrangian data.
    file_name = hier_dump_dirname + "/" + "fe_mesh.";
    sprintf(temp_buf, "%05d", iteration_num);
    file_name += temp_buf;
    file_name += ".xda";
    mesh.write(file_name);

    file_name = hier_dump_dirname + "/" + "fe_equation_systems.";
    sprintf(temp_buf, "%05d", iteration_num);
    file_name += temp_buf;
    equation_systems.write(file_name, (EquationSystems::WRITE_DATA | EquationSystems::WRITE_ADDITIONAL_DATA));
    return;
}// dump_hier_data
}
using namespace Postprocessing;

int
main(
    int argc,
    char* argv[])
{
    if (argc == 1)
    {
        pout << "USAGE:  " << argv[0] << " <input filename> <restart dir> <restore number> [options]\n"
             << "OPTIONS: PETSc command line options; use -help for more information\n";
        return -1;
    }

    // Initialize libMesh, PETSc, MPI, and SAMRAI.
    LibMeshInit init(argc, argv);
    SAMRAI_MPI::setCallAbortInSerialInsteadOfExit();
    SAMRAI_MPI::setCommunicator(PETSC_COMM_WORLD);
    SAMRAIManager::startup();

    // Process command line options and enable logging.
    const string input_filename = argv[1];
    string restart_read_dirname;
    int restore_num = 0;
    bool is_from_restart = false;
    if (argc >= 4)
    {
        // Check whether this appears to be a restarted run.
        FILE* fstream = (SAMRAI_MPI::getRank() == 0 ? fopen(argv[2], "r") : NULL);
        if (SAMRAI_MPI::bcast(fstream != NULL ? 1 : 0, 0) == 1)
        {
            restart_read_dirname = argv[2];
            restore_num = atoi(argv[3]);
            is_from_restart = true;
        }
        if (fstream != NULL)
        {
            fclose(fstream);
        }
    }

    // Create input database and parse all data in input file.
    Pointer<Database> input_db = new InputDatabase("input_db");
    InputManager::getManager()->parseInputFile(input_filename, input_db);

    // Create a simple FE mesh with periodic boundary conditions in the "x"
    // direction.
    //
    // Note that periodic boundary conditions must be registered with each FE
    // system below once they are created after the call to
    // IBFEHierarchyIntegrator::initializeHierarchyIntegrator().
    Mesh mesh(NDIM);
    const double R = 0.25;
    const double w = 0.0625;
    const int M = input_db->getIntegerWithDefault("M", 4);
    string elem_type = input_db->getStringWithDefault("elem_type", "QUAD4");
    MeshTools::Generation::build_square(mesh,
                                        28*M, M,
                                        0.0, 2.0*M_PI*R,
                                        0.0, w,
                                        Utility::string_to_enum<ElemType>(elem_type));
    VectorValue<double> boundary_translation(2.0*M_PI*R, 0.0, 0.0);
    PeriodicBoundary pbc(boundary_translation);
    pbc.myboundary = 3;
    pbc.pairedboundary = 1;

    // Configure stress tensor options.
    smooth_case = input_db->getBoolWithDefault("smooth_case", smooth_case);

    // Create the FE data manager that manages mappings between the FE mesh and
    // the Cartesian grid.
    const string quad_type = input_db->getStringWithDefault("quad_type", "QGAUSS");
    const string quad_order = input_db->getStringWithDefault("quad_order", "SIXTH");
    AutoPtr<QBase> qrule = QBase::build(Utility::string_to_enum<QuadratureType>(quad_type),NDIM,Utility::string_to_enum<Order>(quad_order));
    const string weighting_fcn = input_db->getStringWithDefault("weighting_fcn", "IB_4");
    const bool use_consistent_mass_matrix = input_db->getBoolWithDefault("use_consistent_mass_matrix", true);
    FEDataManager* const fe_data_manager = FEDataManager::getManager("IBFE Manager", weighting_fcn, weighting_fcn, qrule.get(), use_consistent_mass_matrix);
    const int mesh_level_number = input_db->getInteger("MAX_LEVELS")-1;
    EquationSystems equation_systems(mesh);
    fe_data_manager->setEquationSystems(&equation_systems, mesh_level_number);

    // Process "Main" section of the input database.
    Pointer<Database> main_db = input_db->getDatabase("Main");

    // Configure logging options.
    const string log_file_name = main_db->getStringWithDefault("log_file_name","IB.log");
    const bool log_all_nodes = main_db->getBoolWithDefault("log_all_nodes",false);
    if (log_all_nodes)
    {
        PIO::logAllNodes(log_file_name);
    }
    else
    {
        PIO::logOnlyNodeZero(log_file_name);
    }

    // Configure visualization options.
    const int viz_dump_interval = main_db->getIntegerWithDefault("viz_dump_interval",0);
    const bool viz_dump_data = viz_dump_interval > 0;
    string viz_dump_dirname;
    bool uses_visit = false;
    bool uses_exodus = false;
    int visit_number_procs_per_file = 1;
    string exodus_filename;
    if (viz_dump_data)
    {
        Array<string> viz_writer;
        if (main_db->keyExists("viz_writer"))
        {
            viz_writer = main_db->getStringArray("viz_writer");
        }
        for (int i = 0; i < viz_writer.getSize(); i++)
        {
            if (viz_writer[i] == "VisIt"   ) uses_visit  = true;
            if (viz_writer[i] == "ExodusII") uses_exodus = true;
        }

        if (main_db->keyExists("viz_dump_dirname"))
        {
            viz_dump_dirname = main_db->getString("viz_dump_dirname");
            if (viz_dump_dirname.empty())
            {
                TBOX_ERROR("viz_dump_interval > 0, but `viz_dump_dirname' is empty" << endl);
            }
        }
        else
        {
            TBOX_ERROR("viz_dump_interval > 0, but key `viz_dump_dirname' not specifed in input file" << endl);
        }

        if (uses_visit)
        {
            visit_number_procs_per_file = main_db->getIntegerWithDefault("visit_number_procs_per_file",visit_number_procs_per_file);
        }

        if (uses_exodus)
        {
            exodus_filename = main_db->getStringWithDefault("exodus_filename","output.ex2");
            ostringstream os;
            os << viz_dump_dirname << "/" << exodus_filename;
            exodus_filename = os.str();
        }
    }

    // Configure restart options.
    const int restart_interval = main_db->getIntegerWithDefault("restart_interval",0);
    const bool write_restart = restart_interval > 0;
    string restart_write_dirname;
    if (write_restart)
    {
        if (main_db->keyExists("restart_write_dirname"))
        {
            restart_write_dirname = main_db->getString("restart_write_dirname");
            if (restart_write_dirname.empty())
            {
                TBOX_ERROR("restart_interval > 0, but `restart_write_dirname' is empty" << endl);
            }
        }
        else
        {
            TBOX_ERROR("restart_interval > 0, but key `restart_write_dirname' not specifed in input file" << endl);
        }
    }

    // Configure timing options.
    const int timer_dump_interval = main_db->getIntegerWithDefault("timer_dump_interval",0);
    const bool write_timer_data = timer_dump_interval > 0;
    if (write_timer_data)
    {
        TimerManager::createManager(input_db->getDatabase("TimerManager"));
    }

    // Configure load balancing options.
    const bool use_nonuniform_load_balancer = main_db->getBoolWithDefault("use_nonuniform_load_balancer", false);

    // Configure postprocessing options.
    const int hier_dump_interval = main_db->getIntegerWithDefault("hier_dump_interval",0);
    const bool write_hier_data = hier_dump_interval > 0;
    string hier_dump_dirname;
    if (write_hier_data)
    {
        if (main_db->keyExists("hier_dump_dirname"))
        {
            hier_dump_dirname = main_db->getString("hier_dump_dirname");
            if (hier_dump_dirname.empty())
            {
                TBOX_ERROR("hier_dump_interval > 0, but `hier_dump_dirname' is empty" << endl);
            }
        }
        else
        {
            TBOX_ERROR("hier_dump_interval > 0, but key `hier_dump_dirname' not specifed in input file" << endl);
        }
    }

    // Process restart data if this is a restarted run.
    if (is_from_restart)
    {
        RestartManager::getManager()->openRestartFile(restart_read_dirname, restore_num, SAMRAI_MPI::getNodes());
    }

    // Create major algorithm and data objects which comprise the application.
    // These objects are configured from the input database and, if this is a
    // restarted run, from the restart database.
    Pointer<CartesianGridGeometry<NDIM> > grid_geometry = new CartesianGridGeometry<NDIM>(
        "CartesianGeometry", input_db->getDatabase("CartesianGeometry"));

    Pointer<PatchHierarchy<NDIM> > patch_hierarchy = new PatchHierarchy<NDIM>(
        "PatchHierarchy", grid_geometry);

    Pointer<INSStaggeredHierarchyIntegrator> navier_stokes_integrator = new INSStaggeredHierarchyIntegrator(
        "INSStaggeredHierarchyIntegrator", input_db->getDatabase("INSStaggeredHierarchyIntegrator"), patch_hierarchy);

    Pointer<IBFEHierarchyIntegrator> time_integrator = new IBFEHierarchyIntegrator(
        "IBFEHierarchyIntegrator", input_db->getDatabase("IBFEHierarchyIntegrator"),
        patch_hierarchy, navier_stokes_integrator, fe_data_manager);

    Pointer<StandardTagAndInitialize<NDIM> > error_detector = new StandardTagAndInitialize<NDIM>(
        "StandardTagAndInitialize", time_integrator, input_db->getDatabase("StandardTagAndInitialize"));

    Pointer<BergerRigoutsos<NDIM> > box_generator = new BergerRigoutsos<NDIM>();

    Pointer<LoadBalancer<NDIM> > load_balancer = new LoadBalancer<NDIM>(
        "LoadBalancer", input_db->getDatabase("LoadBalancer"));

    Pointer<GriddingAlgorithm<NDIM> > gridding_algorithm = new GriddingAlgorithm<NDIM>(
        "GriddingAlgorithm", input_db->getDatabase("GriddingAlgorithm"), error_detector, box_generator, load_balancer);

    // Configure the IBFE solver.
    time_integrator->registerInitialCoordinateMappingFunction(&coordinate_mapping_function);
    time_integrator->registerPK1StressTensorFunction(&PK1_stress_function);
    if (use_nonuniform_load_balancer)
    {
        time_integrator->registerLoadBalancer(load_balancer);
    }

    // Create initial conditions specification objects.  These objects also are
    // used to specify exact solution values for error analysis.
    Pointer<CartGridFunction> u_init = new muParserCartGridFunction("u_init", input_db->getDatabase("VelocityInitialConditions"), grid_geometry);
    navier_stokes_integrator->registerVelocityInitialConditions(u_init);

    Pointer<CartGridFunction> p_init = new muParserCartGridFunction("p_init", input_db->getDatabase("PressureInitialConditions"), grid_geometry);
    navier_stokes_integrator->registerPressureInitialConditions(p_init);

    // Create boundary condition specification objects (when necessary).
    const IntVector<NDIM>& periodic_shift = grid_geometry->getPeriodicShift();
    const bool has_physical_boundaries = periodic_shift.min() == 0;
    blitz::TinyVector<RobinBcCoefStrategy<NDIM>*,NDIM> u_bc_coefs;
    for (unsigned int d = 0; d < NDIM; ++d)
    {
        if (periodic_shift(d) == 0)
        {
            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();

            if (input_db->isDatabase(bc_coefs_db_name))
            {
                u_bc_coefs[d] = new muParserRobinBcCoefs(bc_coefs_name, input_db->getDatabase(bc_coefs_db_name), grid_geometry);
            }
            else
            {
                TBOX_ERROR("dimension " << d << "is nonperiodic.\n" <<
                           "boundary condition data must be provided in database named " << bc_coefs_db_name << endl);
            }
        }
    }
    if (has_physical_boundaries) time_integrator->registerVelocityPhysicalBcCoefs(u_bc_coefs);

    // Create body force function specification objects (when necessary).
    if (input_db->keyExists("ForcingFunction"))
    {
        tbox::Pointer<CartGridFunction> f_fcn = new muParserCartGridFunction(
            "f_fcn", input_db->getDatabase("ForcingFunction"), grid_geometry);
        time_integrator->registerEulBodyForceFunction(f_fcn);
    }

    // Setup visualization plot file writers.
    Pointer<VisItDataWriter<NDIM> > visit_data_writer;
    if (uses_visit)
    {
        visit_data_writer = new VisItDataWriter<NDIM>(
            "VisIt Writer",
            viz_dump_dirname, visit_number_procs_per_file);
        time_integrator->registerVisItDataWriter(visit_data_writer);
    }
    AutoPtr<ExodusII_IO> exodus_io(uses_exodus ? new ExodusII_IO(mesh) : NULL);

    // Initialize hierarchy configuration and data on all patches.
    time_integrator->initializeHierarchyIntegrator(gridding_algorithm);
    for (unsigned int k = 0; k < equation_systems.n_systems(); ++k)
    {
        System& system = equation_systems.get_system(k);
        system.get_dof_map().add_periodic_boundary(pbc);
    }
    double dt_now = time_integrator->initializeHierarchy();

    // Close the restart manager.
    RestartManager::getManager()->closeRestartFile();

    // Print the input database contents to the log file.
    plog << "Input database:" << endl;
    input_db->printClassData(plog);

    // Setup data used to determine the accuracy of the computed solution.
    VariableDatabase<NDIM>* var_db = VariableDatabase<NDIM>::getDatabase();

    const Pointer<SideVariable<NDIM,double> > u_var = navier_stokes_integrator->getVelocityVar();
    const Pointer<VariableContext> u_ctx = navier_stokes_integrator->getCurrentContext();

    const int u_idx = var_db->mapVariableAndContextToIndex(u_var, u_ctx);
    const int u_cloned_idx = var_db->registerClonedPatchDataIndex(u_var, u_idx);

    const Pointer<CellVariable<NDIM,double> > p_var = navier_stokes_integrator->getPressureVar();
    const Pointer<VariableContext> p_ctx = navier_stokes_integrator->getCurrentContext();

    const int p_idx = var_db->mapVariableAndContextToIndex(p_var, p_ctx);
    const int p_cloned_idx = var_db->registerClonedPatchDataIndex(p_var, p_idx);
    visit_data_writer->registerPlotQuantity("P error", "SCALAR", p_cloned_idx);

    // Indicate the Eulerian data components to be saved for postprocessing.
    const int U_data_idx = var_db->mapVariableAndContextToIndex(
        navier_stokes_integrator->getVelocityVar(), navier_stokes_integrator->getCurrentContext());
    const int P_data_idx = var_db->mapVariableAndContextToIndex(
        navier_stokes_integrator->getPressureVar(), navier_stokes_integrator->getCurrentContext());
    const int P_extrap_data_idx = var_db->mapVariableAndContextToIndex(
        navier_stokes_integrator->getExtrapolatedPressureVar(), navier_stokes_integrator->getCurrentContext());
    hier::ComponentSelector hier_data_comps;
    hier_data_comps.setFlag(U_data_idx);
    hier_data_comps.setFlag(P_data_idx);
    hier_data_comps.setFlag(P_extrap_data_idx);

    // Write out initial data.
    double loop_time = time_integrator->getIntegratorTime();
    int iteration_num = time_integrator->getIntegratorStep();
    if (viz_dump_data && iteration_num%viz_dump_interval == 0)
    {
        pout << "\nWriting visualization files...\n\n";
        if (uses_visit)
        {
            visit_data_writer->writePlotData(patch_hierarchy, iteration_num, loop_time);
        }
        if (uses_exodus)
        {
            exodus_io->write_timestep(exodus_filename, equation_systems, iteration_num/viz_dump_interval+1, loop_time);
        }
    }
    if (write_hier_data && iteration_num%hier_dump_interval == 0)
    {
        pout << "\nWriting data files for postprocessing...\n\n";
        dump_hier_data(iteration_num, loop_time, dt_now, hier_dump_dirname, hier_dump_interval,
                       patch_hierarchy, hier_data_comps, mesh, equation_systems);
    }

    // Main time step loop.
    const double loop_time_end = time_integrator->getEndTime();
    while (!MathUtilities<double>::equalEps(loop_time,loop_time_end) && time_integrator->stepsRemaining())
    {
        iteration_num = time_integrator->getIntegratorStep() + 1;

        pout <<                                                       endl;
        pout << "++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
        pout << "At beginning of timestep # " << iteration_num - 1 << endl;
        pout << "Simulation time is " << loop_time                 << endl;

        double dt_new = time_integrator->advanceHierarchy(dt_now);
        loop_time += dt_now;
        dt_now = dt_new;

        pout <<                                                       endl;
        pout << "At end       of timestep # " << iteration_num - 1 << endl;
        pout << "Simulation time is " << loop_time                 << endl;
        pout << "++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
        pout <<                                                       endl;

        // Compute velocity and pressure error norms.
        const int coarsest_ln = 0;
        const int finest_ln = patch_hierarchy->getFinestLevelNumber();
        for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
        {
            Pointer<PatchLevel<NDIM> > level = patch_hierarchy->getPatchLevel(ln);
            if (!level->checkAllocated(u_cloned_idx)) level->allocatePatchData(u_cloned_idx, 0.0);
            if (!level->checkAllocated(p_cloned_idx)) level->allocatePatchData(p_cloned_idx, 0.0);
        }

        u_init->setDataOnPatchHierarchy(u_cloned_idx, u_var, patch_hierarchy, loop_time);
        p_init->setDataOnPatchHierarchy(p_cloned_idx, p_var, patch_hierarchy, loop_time-0.5*dt_now);

        HierarchyMathOps hier_math_ops("HierarchyMathOps", patch_hierarchy);
        hier_math_ops.setPatchHierarchy(patch_hierarchy);
        hier_math_ops.resetLevels(coarsest_ln, finest_ln);
        const double volume  = hier_math_ops.getVolumeOfPhysicalDomain();
        const int wgt_cc_idx = hier_math_ops.getCellWeightPatchDescriptorIndex();
        const int wgt_sc_idx = hier_math_ops.getSideWeightPatchDescriptorIndex();

        HierarchySideDataOpsReal<NDIM,double> hier_sc_data_ops(patch_hierarchy, coarsest_ln, finest_ln);
        hier_sc_data_ops.subtract(u_cloned_idx, u_idx, u_cloned_idx);
        for (int ln = patch_hierarchy->getFinestLevelNumber(); ln > 0; --ln)
        {
            Pointer<CoarsenOperator<NDIM> > coarsen_operator = grid_geometry->lookupCoarsenOperator(u_var, "CONSERVATIVE_COARSEN");
            CoarsenAlgorithm<NDIM> calg;
            calg.registerCoarsen(u_cloned_idx, u_cloned_idx, coarsen_operator);
            Pointer<PatchLevel<NDIM> > level = patch_hierarchy->getPatchLevel(ln);
            Pointer<PatchLevel<NDIM> > coarser_level = patch_hierarchy->getPatchLevel(ln-1);
            calg.createSchedule(coarser_level, level)->coarsenData();
        }
        pout << "Error in u at time " << loop_time << ":\n"
             << "  L1-norm:  " << hier_sc_data_ops.L1Norm(u_cloned_idx,wgt_sc_idx)  << "\n"
             << "  L2-norm:  " << hier_sc_data_ops.L2Norm(u_cloned_idx,wgt_sc_idx)  << "\n"
             << "  max-norm: " << hier_sc_data_ops.maxNorm(u_cloned_idx,wgt_sc_idx) << "\n";

        HierarchyCellDataOpsReal<NDIM,double> hier_cc_data_ops(patch_hierarchy, coarsest_ln, finest_ln);
        const double p_mean = (1.0/volume)*hier_cc_data_ops.integral(p_idx, wgt_cc_idx);
        hier_cc_data_ops.addScalar(p_idx, p_idx, -p_mean);
        const double p_cloned_mean = (1.0/volume)*hier_cc_data_ops.integral(p_cloned_idx, wgt_cc_idx);
        hier_cc_data_ops.addScalar(p_cloned_idx, p_cloned_idx, -p_cloned_mean);
        hier_cc_data_ops.subtract(p_cloned_idx, p_idx, p_cloned_idx);
        for (int ln = patch_hierarchy->getFinestLevelNumber(); ln > 0; --ln)
        {
            Pointer<CoarsenOperator<NDIM> > coarsen_operator = grid_geometry->lookupCoarsenOperator(p_var, "CONSERVATIVE_COARSEN");
            CoarsenAlgorithm<NDIM> calg;
            calg.registerCoarsen(p_cloned_idx, p_cloned_idx, coarsen_operator);
            Pointer<PatchLevel<NDIM> > level = patch_hierarchy->getPatchLevel(ln);
            Pointer<PatchLevel<NDIM> > coarser_level = patch_hierarchy->getPatchLevel(ln-1);
            calg.createSchedule(coarser_level, level)->coarsenData();
        }
        pout << "Error in p at time " << loop_time-0.5*dt_now << ":\n"
             << "  L1-norm:  " << hier_cc_data_ops.L1Norm(p_cloned_idx,wgt_cc_idx)  << "\n"
             << "  L2-norm:  " << hier_cc_data_ops.L2Norm(p_cloned_idx,wgt_cc_idx)  << "\n"
             << "  max-norm: " << hier_cc_data_ops.maxNorm(p_cloned_idx,wgt_cc_idx) << "\n";

        // At specified intervals, write visualization, restart, and
        // postprocessing files, and print out timer data.
        if (viz_dump_data && iteration_num%viz_dump_interval == 0)
        {
            pout << "\nWriting visualization files...\n\n";
            if (uses_visit)
            {
                visit_data_writer->writePlotData(patch_hierarchy, iteration_num, loop_time);
            }
            if (uses_exodus)
            {
                exodus_io->write_timestep(exodus_filename, equation_systems, iteration_num/viz_dump_interval+1, loop_time);
            }
        }

        if (write_restart && iteration_num%restart_interval == 0)
        {
            pout << "\nWriting restart files...\n\n";
            RestartManager::getManager()->writeRestartFile(restart_write_dirname, iteration_num);
        }

        if (write_hier_data && iteration_num%hier_dump_interval == 0)
        {
            pout << "\nWriting data files for postprocessing...\n\n";
            dump_hier_data(iteration_num, loop_time, dt_now, hier_dump_dirname, hier_dump_interval,
                           patch_hierarchy, hier_data_comps, mesh, equation_systems);
        }

        if (write_timer_data && iteration_num%timer_dump_interval == 0)
        {
            pout << "\nWriting timer data...\n\n";
            TimerManager::getManager()->print(plog);
        }
    }

    for (unsigned int d = 0; d < NDIM; ++d)
    {
        if (periodic_shift(d) == 0) delete u_bc_coefs[d];
    }

    // Shutdown SAMRAI.
    SAMRAIManager::shutdown();
    return 0;
}// main
<commit_msg>branch boyceg: reverted IBFE/explicit/ex0/main.C to r2433<commit_after>// Copyright (c) 2002-2010, <NAME>
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
//    * Redistributions of source code must retain the above copyright notice,
//      this list of conditions and the following disclaimer.
//
//    * Redistributions in binary form must reproduce the above copyright
//      notice, this list of conditions and the following disclaimer in the
//      documentation and/or other materials provided with the distribution.
//
//    * Neither the name of New York University nor the names of its
//      contributors may be used to endorse or promote products derived from
//      this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.

// Config files
#include <IBAMR_config.h>
#include <SAMRAI_config.h>

// Headers for basic PETSc functions
#include <petsc.h>

// Headers for basic libMesh objects
#include <exodusII_io.h>
#include <mesh.h>
#include <mesh_generation.h>
#include <periodic_boundaries.h>
#include <quadrature.h>
#include <string_to_enum.h>

// Headers for application-specific algorithm/data structure objects
#include <ibamr/IBFEHierarchyIntegrator.h>
#include <ibtk/muParserCartGridFunction.h>
#include <ibtk/muParserRobinBcCoefs.h>
#include <BergerRigoutsos.h>
#include <StandardTagAndInitialize.h>

// C++ namespace delcarations
#include <ibamr/namespaces.h>
using namespace IBAMR;
using namespace IBTK;
using namespace libMesh;
using namespace std;

// Elasticity model data.
namespace ModelData
{
// Problem parameters.
static const double R = 0.25;
static const double w = 0.0625;
static const double gamma = 0.0;
static const double mu = 1.0;

// Coordinate mapping function.
void
coordinate_mapping_function(
    Point& X,
    const Point& s,
    void* /*ctx*/)
{
    X(0) = (R+      s(1))*cos(s(0)/R)+0.5;
    X(1) = (R+gamma+s(1))*sin(s(0)/R)+0.5;
    return;
}// coordinate_mapping_function

// Stress tensor function.
bool smooth_case = false;
void
PK1_stress_function(
    TensorValue<double>& PP,
    const TensorValue<double>& dX_ds,
    const Point& /*X*/,
    const Point& /*s*/,
    Elem* const /*elem*/,
    NumericVector<double>& /*X_vec*/,
    const std::vector<NumericVector<double>*>& /*system_data*/,
    const double& /*time*/,
    void* /*ctx*/)
{
    PP = (mu/w)*dX_ds;
    if (smooth_case)
    {
        PP(0,1) = 0.0;
        PP(1,1) = 0.0;
    }
    return;
}// PK1_stress_function
}
using namespace ModelData;

namespace Postprocessing
{
void
dump_hier_data(
    const int iteration_num,
    const double loop_time,
    const double dt,
    const string hier_dump_dirname,
    const int hier_dump_interval,
    Pointer<PatchHierarchy<NDIM> > patch_hierarchy,
    const ComponentSelector& hier_data,
    Mesh& mesh,
    EquationSystems& equation_systems)
{
    string file_name;

    // Write Cartesian data.
    file_name = hier_dump_dirname + "/" + "hier_data.";
    char temp_buf[128];
    sprintf(temp_buf, "%05d.samrai.%05d", iteration_num, SAMRAI_MPI::getRank());
    file_name += temp_buf;

    Pointer<HDFDatabase> hier_db = new HDFDatabase("hier_db");
    hier_db->create(file_name);
    patch_hierarchy->putToDatabase(hier_db->putDatabase("PatchHierarchy"), hier_data);
    hier_db->putDouble("loop_time", loop_time);
    hier_db->putDouble("dt", dt);
    hier_db->putInteger("iteration_num", iteration_num);
    hier_db->putInteger("hier_dump_interval", hier_dump_interval);
    hier_db->close();

    // Write Lagrangian data.
    file_name = hier_dump_dirname + "/" + "fe_mesh.";
    sprintf(temp_buf, "%05d", iteration_num);
    file_name += temp_buf;
    file_name += ".xda";
    mesh.write(file_name);

    file_name = hier_dump_dirname + "/" + "fe_equation_systems.";
    sprintf(temp_buf, "%05d", iteration_num);
    file_name += temp_buf;
    equation_systems.write(file_name, (EquationSystems::WRITE_DATA | EquationSystems::WRITE_ADDITIONAL_DATA));
    return;
}// dump_hier_data
}
using namespace Postprocessing;

int
main(
    int argc,
    char* argv[])
{
    if (argc == 1)
    {
        pout << "USAGE:  " << argv[0] << " <input filename> <restart dir> <restore number> [options]\n"
             << "OPTIONS: PETSc command line options; use -help for more information\n";
        return -1;
    }

    // Initialize libMesh, PETSc, MPI, and SAMRAI.
    LibMeshInit init(argc, argv);
    SAMRAI_MPI::setCallAbortInSerialInsteadOfExit();
    SAMRAI_MPI::setCommunicator(PETSC_COMM_WORLD);
    SAMRAIManager::startup();

    // Process command line options and enable logging.
    const string input_filename = argv[1];
    string restart_read_dirname;
    int restore_num = 0;
    bool is_from_restart = false;
    if (argc >= 4)
    {
        // Check whether this appears to be a restarted run.
        FILE* fstream = (SAMRAI_MPI::getRank() == 0 ? fopen(argv[2], "r") : NULL);
        if (SAMRAI_MPI::bcast(fstream != NULL ? 1 : 0, 0) == 1)
        {
            restart_read_dirname = argv[2];
            restore_num = atoi(argv[3]);
            is_from_restart = true;
        }
        if (fstream != NULL)
        {
            fclose(fstream);
        }
    }

    // Create input database and parse all data in input file.
    Pointer<Database> input_db = new InputDatabase("input_db");
    InputManager::getManager()->parseInputFile(input_filename, input_db);

    // Create a simple FE mesh with periodic boundary conditions in the "x"
    // direction.
    //
    // Note that periodic boundary conditions must be registered with each FE
    // system below once they are created after the call to
    // IBFEHierarchyIntegrator::initializeHierarchyIntegrator().
    Mesh mesh(NDIM);
    const double R = 0.25;
    const double w = 0.0625;
    const int M = input_db->getIntegerWithDefault("M", 4);
    string elem_type = input_db->getStringWithDefault("elem_type", "QUAD4");
    MeshTools::Generation::build_square(mesh,
                                        28*M, M,
                                        0.0, 2.0*M_PI*R,
                                        0.0, w,
                                        Utility::string_to_enum<ElemType>(elem_type));
    VectorValue<double> boundary_translation(2.0*M_PI*R, 0.0, 0.0);
    PeriodicBoundary pbc(boundary_translation);
    pbc.myboundary = 3;
    pbc.pairedboundary = 1;

    // Configure stress tensor options.
    smooth_case = input_db->getBoolWithDefault("smooth_case", smooth_case);

    // Create the FE data manager that manages mappings between the FE mesh and
    // the Cartesian grid.
    const string quad_type = input_db->getStringWithDefault("quad_type", "QGAUSS");
    const string quad_order = input_db->getStringWithDefault("quad_order", "SIXTH");
    AutoPtr<QBase> qrule = QBase::build(Utility::string_to_enum<QuadratureType>(quad_type),NDIM,Utility::string_to_enum<Order>(quad_order));
    const string weighting_fcn = input_db->getStringWithDefault("weighting_fcn", "IB_4");
    const bool use_consistent_mass_matrix = input_db->getBoolWithDefault("use_consistent_mass_matrix", true);
    FEDataManager* const fe_data_manager = FEDataManager::getManager("IBFE Manager", weighting_fcn, weighting_fcn, qrule.get(), use_consistent_mass_matrix);
    const int mesh_level_number = input_db->getInteger("MAX_LEVELS")-1;
    EquationSystems equation_systems(mesh);
    fe_data_manager->setEquationSystems(&equation_systems, mesh_level_number);

    // Process "Main" section of the input database.
    Pointer<Database> main_db = input_db->getDatabase("Main");

    // Configure logging options.
    const string log_file_name = main_db->getStringWithDefault("log_file_name","IB.log");
    const bool log_all_nodes = main_db->getBoolWithDefault("log_all_nodes",false);
    if (log_all_nodes)
    {
        PIO::logAllNodes(log_file_name);
    }
    else
    {
        PIO::logOnlyNodeZero(log_file_name);
    }

    // Configure visualization options.
    const int viz_dump_interval = main_db->getIntegerWithDefault("viz_dump_interval",0);
    const bool viz_dump_data = viz_dump_interval > 0;
    string viz_dump_dirname;
    bool uses_visit = false;
    bool uses_exodus = false;
    int visit_number_procs_per_file = 1;
    string exodus_filename;
    if (viz_dump_data)
    {
        Array<string> viz_writer;
        if (main_db->keyExists("viz_writer"))
        {
            viz_writer = main_db->getStringArray("viz_writer");
        }
        for (int i = 0; i < viz_writer.getSize(); i++)
        {
            if (viz_writer[i] == "VisIt"   ) uses_visit  = true;
            if (viz_writer[i] == "ExodusII") uses_exodus = true;
        }

        if (main_db->keyExists("viz_dump_dirname"))
        {
            viz_dump_dirname = main_db->getString("viz_dump_dirname");
            if (viz_dump_dirname.empty())
            {
                TBOX_ERROR("viz_dump_interval > 0, but `viz_dump_dirname' is empty" << endl);
            }
        }
        else
        {
            TBOX_ERROR("viz_dump_interval > 0, but key `viz_dump_dirname' not specifed in input file" << endl);
        }

        if (uses_visit)
        {
            visit_number_procs_per_file = main_db->getIntegerWithDefault("visit_number_procs_per_file",visit_number_procs_per_file);
        }

        if (uses_exodus)
        {
            exodus_filename = main_db->getStringWithDefault("exodus_filename","output.ex2");
            ostringstream os;
            os << viz_dump_dirname << "/" << exodus_filename;
            exodus_filename = os.str();
        }
    }

    // Configure restart options.
    const int restart_interval = main_db->getIntegerWithDefault("restart_interval",0);
    const bool write_restart = restart_interval > 0;
    string restart_write_dirname;
    if (write_restart)
    {
        if (main_db->keyExists("restart_write_dirname"))
        {
            restart_write_dirname = main_db->getString("restart_write_dirname");
            if (restart_write_dirname.empty())
            {
                TBOX_ERROR("restart_interval > 0, but `restart_write_dirname' is empty" << endl);
            }
        }
        else
        {
            TBOX_ERROR("restart_interval > 0, but key `restart_write_dirname' not specifed in input file" << endl);
        }
    }

    // Configure timing options.
    const int timer_dump_interval = main_db->getIntegerWithDefault("timer_dump_interval",0);
    const bool write_timer_data = timer_dump_interval > 0;
    if (write_timer_data)
    {
        TimerManager::createManager(input_db->getDatabase("TimerManager"));
    }

    // Configure load balancing options.
    const bool use_nonuniform_load_balancer = main_db->getBoolWithDefault("use_nonuniform_load_balancer", false);

    // Configure postprocessing options.
    const int hier_dump_interval = main_db->getIntegerWithDefault("hier_dump_interval",0);
    const bool write_hier_data = hier_dump_interval > 0;
    string hier_dump_dirname;
    if (write_hier_data)
    {
        if (main_db->keyExists("hier_dump_dirname"))
        {
            hier_dump_dirname = main_db->getString("hier_dump_dirname");
            if (hier_dump_dirname.empty())
            {
                TBOX_ERROR("hier_dump_interval > 0, but `hier_dump_dirname' is empty" << endl);
            }
        }
        else
        {
            TBOX_ERROR("hier_dump_interval > 0, but key `hier_dump_dirname' not specifed in input file" << endl);
        }
    }

    // Process restart data if this is a restarted run.
    if (is_from_restart)
    {
        RestartManager::getManager()->openRestartFile(restart_read_dirname, restore_num, SAMRAI_MPI::getNodes());
    }

    // Create major algorithm and data objects which comprise the application.
    // These objects are configured from the input database and, if this is a
    // restarted run, from the restart database.
    Pointer<CartesianGridGeometry<NDIM> > grid_geometry = new CartesianGridGeometry<NDIM>(
        "CartesianGeometry", input_db->getDatabase("CartesianGeometry"));

    Pointer<PatchHierarchy<NDIM> > patch_hierarchy = new PatchHierarchy<NDIM>(
        "PatchHierarchy", grid_geometry);

    Pointer<INSStaggeredHierarchyIntegrator> navier_stokes_integrator = new INSStaggeredHierarchyIntegrator(
        "INSStaggeredHierarchyIntegrator", input_db->getDatabase("INSStaggeredHierarchyIntegrator"), patch_hierarchy);

    Pointer<IBFEHierarchyIntegrator> time_integrator = new IBFEHierarchyIntegrator(
        "IBFEHierarchyIntegrator", input_db->getDatabase("IBFEHierarchyIntegrator"),
        patch_hierarchy, navier_stokes_integrator, fe_data_manager);

    Pointer<StandardTagAndInitialize<NDIM> > error_detector = new StandardTagAndInitialize<NDIM>(
        "StandardTagAndInitialize", time_integrator, input_db->getDatabase("StandardTagAndInitialize"));

    Pointer<BergerRigoutsos<NDIM> > box_generator = new BergerRigoutsos<NDIM>();

    Pointer<LoadBalancer<NDIM> > load_balancer = new LoadBalancer<NDIM>(
        "LoadBalancer", input_db->getDatabase("LoadBalancer"));

    Pointer<GriddingAlgorithm<NDIM> > gridding_algorithm = new GriddingAlgorithm<NDIM>(
        "GriddingAlgorithm", input_db->getDatabase("GriddingAlgorithm"), error_detector, box_generator, load_balancer);

    // Configure the IBFE solver.
    time_integrator->registerInitialCoordinateMappingFunction(&coordinate_mapping_function);
    time_integrator->registerPK1StressTensorFunction(&PK1_stress_function);
    if (use_nonuniform_load_balancer)
    {
        time_integrator->registerLoadBalancer(load_balancer);
    }

    // Create initial conditions specification objects.  These objects also are
    // used to specify exact solution values for error analysis.
    Pointer<CartGridFunction> u_init = new muParserCartGridFunction("u_init", input_db->getDatabase("VelocityInitialConditions"), grid_geometry);
    navier_stokes_integrator->registerVelocityInitialConditions(u_init);

    Pointer<CartGridFunction> p_init = new muParserCartGridFunction("p_init", input_db->getDatabase("PressureInitialConditions"), grid_geometry);
    navier_stokes_integrator->registerPressureInitialConditions(p_init);

    // Create boundary condition specification objects (when necessary).
    const IntVector<NDIM>& periodic_shift = grid_geometry->getPeriodicShift();
    const bool has_physical_boundaries = periodic_shift.min() == 0;
    blitz::TinyVector<RobinBcCoefStrategy<NDIM>*,NDIM> u_bc_coefs;
    for (unsigned int d = 0; d < NDIM; ++d)
    {
        if (periodic_shift(d) == 0)
        {
            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();

            if (input_db->isDatabase(bc_coefs_db_name))
            {
                u_bc_coefs[d] = new muParserRobinBcCoefs(bc_coefs_name, input_db->getDatabase(bc_coefs_db_name), grid_geometry);
            }
            else
            {
                TBOX_ERROR("dimension " << d << "is nonperiodic.\n" <<
                           "boundary condition data must be provided in database named " << bc_coefs_db_name << endl);
            }
        }
    }
    if (has_physical_boundaries) time_integrator->registerVelocityPhysicalBcCoefs(u_bc_coefs);

    // Create body force function specification objects (when necessary).
    if (input_db->keyExists("ForcingFunction"))
    {
        tbox::Pointer<CartGridFunction> f_fcn = new muParserCartGridFunction(
            "f_fcn", input_db->getDatabase("ForcingFunction"), grid_geometry);
        time_integrator->registerEulBodyForceFunction(f_fcn);
    }

    // Setup visualization plot file writers.
    Pointer<VisItDataWriter<NDIM> > visit_data_writer;
    if (uses_visit)
    {
        visit_data_writer = new VisItDataWriter<NDIM>(
            "VisIt Writer",
            viz_dump_dirname, visit_number_procs_per_file);
        time_integrator->registerVisItDataWriter(visit_data_writer);
    }
    AutoPtr<ExodusII_IO> exodus_io(uses_exodus ? new ExodusII_IO(mesh) : NULL);

    // Initialize hierarchy configuration and data on all patches.
    time_integrator->initializeHierarchyIntegrator(gridding_algorithm);
    for (unsigned int k = 0; k < equation_systems.n_systems(); ++k)
    {
        System& system = equation_systems.get_system(k);
        system.get_dof_map().add_periodic_boundary(pbc);
    }
    double dt_now = time_integrator->initializeHierarchy();

    // Close the restart manager.
    RestartManager::getManager()->closeRestartFile();

    // Print the input database contents to the log file.
    plog << "Input database:" << endl;
    input_db->printClassData(plog);

    // Setup data used to determine the accuracy of the computed solution.
    VariableDatabase<NDIM>* var_db = VariableDatabase<NDIM>::getDatabase();

    const Pointer<SideVariable<NDIM,double> > u_var = navier_stokes_integrator->getVelocityVar();
    const Pointer<VariableContext> u_ctx = navier_stokes_integrator->getCurrentContext();

    const int u_idx = var_db->mapVariableAndContextToIndex(u_var, u_ctx);
    const int u_cloned_idx = var_db->registerClonedPatchDataIndex(u_var, u_idx);

    const Pointer<CellVariable<NDIM,double> > p_var = navier_stokes_integrator->getPressureVar();
    const Pointer<VariableContext> p_ctx = navier_stokes_integrator->getCurrentContext();

    const int p_idx = var_db->mapVariableAndContextToIndex(p_var, p_ctx);
    const int p_cloned_idx = var_db->registerClonedPatchDataIndex(p_var, p_idx);
    visit_data_writer->registerPlotQuantity("P error", "SCALAR", p_cloned_idx);

    // Indicate the Eulerian data components to be saved for postprocessing.
    const int U_data_idx = var_db->mapVariableAndContextToIndex(
        navier_stokes_integrator->getVelocityVar(), navier_stokes_integrator->getCurrentContext());
    const int P_data_idx = var_db->mapVariableAndContextToIndex(
        navier_stokes_integrator->getPressureVar(), navier_stokes_integrator->getCurrentContext());
    const int P_extrap_data_idx = var_db->mapVariableAndContextToIndex(
        navier_stokes_integrator->getExtrapolatedPressureVar(), navier_stokes_integrator->getCurrentContext());
    hier::ComponentSelector hier_data_comps;
    hier_data_comps.setFlag(U_data_idx);
    hier_data_comps.setFlag(P_data_idx);
    hier_data_comps.setFlag(P_extrap_data_idx);

    // Write out initial data.
    double loop_time = time_integrator->getIntegratorTime();
    int iteration_num = time_integrator->getIntegratorStep();
    if (viz_dump_data && iteration_num%viz_dump_interval == 0)
    {
        pout << "\nWriting visualization files...\n\n";
        if (uses_visit)
        {
            visit_data_writer->writePlotData(patch_hierarchy, iteration_num, loop_time);
        }
        if (uses_exodus)
        {
            exodus_io->write_timestep(exodus_filename, equation_systems, iteration_num/viz_dump_interval+1, loop_time);
        }
    }
    if (write_hier_data && iteration_num%hier_dump_interval == 0)
    {
        pout << "\nWriting data files for postprocessing...\n\n";
        dump_hier_data(iteration_num, loop_time, dt_now, hier_dump_dirname, hier_dump_interval,
                       patch_hierarchy, hier_data_comps, mesh, equation_systems);
    }

    // Main time step loop.
    const double loop_time_end = time_integrator->getEndTime();
    while (!MathUtilities<double>::equalEps(loop_time,loop_time_end) && time_integrator->stepsRemaining())
    {
        iteration_num = time_integrator->getIntegratorStep() + 1;

        pout <<                                                       endl;
        pout << "++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
        pout << "At beginning of timestep # " << iteration_num - 1 << endl;
        pout << "Simulation time is " << loop_time                 << endl;

        double dt_new = time_integrator->advanceHierarchy(dt_now);
        loop_time += dt_now;
        dt_now = dt_new;

        pout <<                                                       endl;
        pout << "At end       of timestep # " << iteration_num - 1 << endl;
        pout << "Simulation time is " << loop_time                 << endl;
        pout << "++++++++++++++++++++++++++++++++++++++++++++++++" << endl;
        pout <<                                                       endl;

        // Compute velocity and pressure error norms.
        const int coarsest_ln = 0;
        const int finest_ln = patch_hierarchy->getFinestLevelNumber();
        for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
        {
            Pointer<PatchLevel<NDIM> > level = patch_hierarchy->getPatchLevel(ln);
            if (!level->checkAllocated(u_cloned_idx)) level->allocatePatchData(u_cloned_idx, 0.0);
            if (!level->checkAllocated(p_cloned_idx)) level->allocatePatchData(p_cloned_idx, 0.0);
        }

        u_init->setDataOnPatchHierarchy(u_cloned_idx, u_var, patch_hierarchy, loop_time);
        p_init->setDataOnPatchHierarchy(p_cloned_idx, p_var, patch_hierarchy, loop_time-0.5*dt_now);

        HierarchyMathOps hier_math_ops("HierarchyMathOps", patch_hierarchy);
        hier_math_ops.setPatchHierarchy(patch_hierarchy);
        hier_math_ops.resetLevels(coarsest_ln, finest_ln);
        const double volume  = hier_math_ops.getVolumeOfPhysicalDomain();
        const int wgt_cc_idx = hier_math_ops.getCellWeightPatchDescriptorIndex();
        const int wgt_sc_idx = hier_math_ops.getSideWeightPatchDescriptorIndex();

        HierarchySideDataOpsReal<NDIM,double> hier_sc_data_ops(patch_hierarchy, coarsest_ln, finest_ln);
        hier_sc_data_ops.subtract(u_cloned_idx, u_idx, u_cloned_idx);
        for (int ln = patch_hierarchy->getFinestLevelNumber(); ln > 0; --ln)
        {
            Pointer<CoarsenOperator<NDIM> > coarsen_operator = grid_geometry->lookupCoarsenOperator(u_var, "CONSERVATIVE_COARSEN");
            CoarsenAlgorithm<NDIM> calg;
            calg.registerCoarsen(u_cloned_idx, u_cloned_idx, coarsen_operator);
            Pointer<PatchLevel<NDIM> > level = patch_hierarchy->getPatchLevel(ln);
            Pointer<PatchLevel<NDIM> > coarser_level = patch_hierarchy->getPatchLevel(ln-1);
            calg.createSchedule(coarser_level, level)->coarsenData();
        }
        pout << "Error in u at time " << loop_time << ":\n"
             << "  L1-norm:  " << hier_sc_data_ops.L1Norm(u_cloned_idx,wgt_sc_idx)  << "\n"
             << "  L2-norm:  " << hier_sc_data_ops.L2Norm(u_cloned_idx,wgt_sc_idx)  << "\n"
             << "  max-norm: " << hier_sc_data_ops.maxNorm(u_cloned_idx,wgt_sc_idx) << "\n";

        HierarchyCellDataOpsReal<NDIM,double> hier_cc_data_ops(patch_hierarchy, coarsest_ln, finest_ln);
        const double p_mean = (1.0/volume)*hier_cc_data_ops.integral(p_idx, wgt_cc_idx);
        hier_cc_data_ops.addScalar(p_idx, p_idx, -p_mean);
        const double p_cloned_mean = (1.0/volume)*hier_cc_data_ops.integral(p_cloned_idx, wgt_cc_idx);
        hier_cc_data_ops.addScalar(p_cloned_idx, p_cloned_idx, -p_cloned_mean);
        hier_cc_data_ops.subtract(p_cloned_idx, p_idx, p_cloned_idx);
        for (int ln = patch_hierarchy->getFinestLevelNumber(); ln > 0; --ln)
        {
            Pointer<CoarsenOperator<NDIM> > coarsen_operator = grid_geometry->lookupCoarsenOperator(p_var, "CONSERVATIVE_COARSEN");
            CoarsenAlgorithm<NDIM> calg;
            calg.registerCoarsen(p_cloned_idx, p_cloned_idx, coarsen_operator);
            Pointer<PatchLevel<NDIM> > level = patch_hierarchy->getPatchLevel(ln);
            Pointer<PatchLevel<NDIM> > coarser_level = patch_hierarchy->getPatchLevel(ln-1);
            calg.createSchedule(coarser_level, level)->coarsenData();
        }
        pout << "Error in p at time " << loop_time-0.5*dt_now << ":\n"
             << "  L1-norm:  " << hier_cc_data_ops.L1Norm(p_cloned_idx,wgt_cc_idx)  << "\n"
             << "  L2-norm:  " << hier_cc_data_ops.L2Norm(p_cloned_idx,wgt_cc_idx)  << "\n"
             << "  max-norm: " << hier_cc_data_ops.maxNorm(p_cloned_idx,wgt_cc_idx) << "\n";

        // At specified intervals, write visualization, restart, and
        // postprocessing files, and print out timer data.
        if (viz_dump_data && iteration_num%viz_dump_interval == 0)
        {
            pout << "\nWriting visualization files...\n\n";
            if (uses_visit)
            {
                visit_data_writer->writePlotData(patch_hierarchy, iteration_num, loop_time);
            }
            if (uses_exodus)
            {
                exodus_io->write_timestep(exodus_filename, equation_systems, iteration_num/viz_dump_interval+1, loop_time);
            }
        }

        if (write_restart && iteration_num%restart_interval == 0)
        {
            pout << "\nWriting restart files...\n\n";
            RestartManager::getManager()->writeRestartFile(restart_write_dirname, iteration_num);
        }

        if (write_hier_data && iteration_num%hier_dump_interval == 0)
        {
            pout << "\nWriting data files for postprocessing...\n\n";
            dump_hier_data(iteration_num, loop_time, dt_now, hier_dump_dirname, hier_dump_interval,
                           patch_hierarchy, hier_data_comps, mesh, equation_systems);
        }

        if (write_timer_data && iteration_num%timer_dump_interval == 0)
        {
            pout << "\nWriting timer data...\n\n";
            TimerManager::getManager()->print(plog);
        }
    }

    for (unsigned int d = 0; d < NDIM; ++d)
    {
        if (periodic_shift(d) == 0) delete u_bc_coefs[d];
    }

    // Shutdown SAMRAI.
    SAMRAIManager::shutdown();
    return 0;
}// main

{'fraction_non_alphanumeric': 0.10908989036848887, 'fraction_numerical': 0.006518125618661636, 'mean_word_length': 3.1818181818181817, 'pattern_counts': {'":': 4, '<': 489, '<?xml version=': 0, '>': 435, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6794324', 'n_tokens_mistral': 17439, 'n_tokens_neox': 16642, 'n_words': 3385}