MonteCarloElliptic.cpp

#include "MonteCarloElliptic.h"


using namespace std;

void MonteCarloElliptic::initialize() {
    finePermeability = 0.0, coarsePermeability = 0.0;
    fineSolution = 0.0, coarseSolution = 0.0;
    fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
}

void MonteCarloElliptic::method() {
    for (int m = M; m < M + dM; m++) {
        string sampleDirName = mkSampleDirName(m, true);
        mkSampleDir(m, sampleDirName);

        stochasticField->SetSampleDir(l, sampleDirName);
        stochasticField->SetPlot(l, finePlot);
        stochasticField->GetFineSample(l, finePermeability);

        assemble->sampleDir = sampleDirName;
        assemble->plot = finePlot;
        assemble->problem->LoadNewSample(&finePermeability);

        solvePDE(l, m, true, fineQ, fineCost, fineSolution);

        if (plotting) assemble->PlotResults(fineSolution);

        if (!onlyFineLevel) {
            sampleDirName = mkSampleDirName(m, false);
            mkSampleDir(m, sampleDirName);

            stochasticField->SetSampleDir(l, sampleDirName);
            stochasticField->SetPlot(l, coarsePlot);
            stochasticField->GetCoarseSample(l, finePermeability, coarsePermeability);

            assemble->sampleDir = sampleDirName;
            assemble->plot = coarsePlot;
            assemble->problem->LoadNewSample(&coarsePermeability);

            solvePDE(l, 0, false, coarseQ, coarseCost, coarseSolution);

            if (plotting) assemble->PlotResults(coarseSolution);

        } else {
            coarseQ = 0.0, coarseCost = 0.0;
        }
        updateSums(fineCost, fineQ - coarseQ, fineQ);
    }
}

void MonteCarloElliptic::solvePDE(int l,
                                  int m,
                                  bool fineLevel,
                                  double &valueQ,
                                  double &cost,
                                  Vector &solution) {
    Preconditioner *pc = GetPC("SuperLU");
    Solver solver(pc, "GMRES");
    Newton newton(solver);

    newton(*assemble, solution);

    cost = solution.pSize();

    if (functional == "L2") valueQ = assemble->L2(solution);
    if (functional == "Energy") valueQ = assemble->Energy(solution);
    if (functional == "Outflow")
        valueQ = assemble->getInflowOutflow(solution).second;
    if (functional == "Goal") valueQ = assemble->getGoal(solution);
}