Function fixMassLoss2PhaseFlows

src/level_set/LevelSetUtilities.cpp:284–350

} // setLSData

void
fixMassLoss2PhaseFlows(double /*current_time*/,
                       double new_time,
                       bool /*skip_synchronize_new_state_data*/,
                       int /*num_cycles*/,
                       void* ctx)
{
    LevelSetMassLossFixer* mass_fixer = static_cast<LevelSetMassLossFixer*>(ctx);
#if !defined(NDEBUG)
    TBOX_ASSERT(mass_fixer);
#endif
    const LevelSetContainer& ls_container = mass_fixer->getLevelSetContainer();
    Pointer<AdvDiffHierarchyIntegrator> adv_diff_integrator = ls_container.getAdvDiffHierarchyIntegrator();
    const int integrator_step = adv_diff_integrator->getIntegratorStep();
    const int mass_correction_interval = mass_fixer->getCorrectionInterval();

    if (integrator_step % mass_correction_interval != 0) return;

    Pointer<PatchHierarchy<NDIM>> patch_hier = adv_diff_integrator->getPatchHierarchy();
    Pointer<HierarchyMathOps> hier_math_ops = adv_diff_integrator->getHierarchyMathOps();

    const int hier_finest_ln = patch_hier->getFinestLevelNumber();
    const double vol_target = mass_fixer->getTargetVolume();
    const double ncells = ls_container.getInterfaceHalfWidth();

    // NOTE: In practice the level set mass loss would be fixed during the postprocess integrate hierarchy stage.
    // Hence the application time would be the new time and the variable context would be the new context.
    VariableDatabase<NDIM>* var_db = VariableDatabase<NDIM>::getDatabase();
    const int ls_idx =
        var_db->mapVariableAndContextToIndex(ls_container.getLevelSetVariable(), adv_diff_integrator->getNewContext());

    // Carry out the Newton iterations
    double rel_error = 1.0e12;
    int current_iter = 0;

    const double min_rel_error = mass_fixer->getErrorRelTolerance();
    const int max_its = mass_fixer->getMaxIterations();

    double q = 0.0;
    HierarchyCellDataOpsReal<NDIM, double> hier_cc_ops(patch_hier, 0, hier_finest_ln);
    while (rel_error > min_rel_error && current_iter < max_its)
    {
        std::vector<double> integrals = compute_heaviside_integrals(hier_math_ops, ls_idx, ncells);
        const double& vol_phase1 = integrals[0]; // Target the gas volume
        const double& integral_delta = integrals[2];

        rel_error = std::abs(vol_phase1 / vol_target - 1.0);

        if (mass_fixer->enableLogging())
        {
            plog << "fixMassLoss2PhaseFlows():: current iter = " << current_iter << " , rel error  = " << rel_error
                 << std::endl;
        }

        const double delta_q = (vol_phase1 - vol_target) / integral_delta;
        hier_cc_ops.addScalar(ls_idx, ls_idx, delta_q);

        q += delta_q;