worked on PDESolver

mlmc/src/montecarlo/PDESolver.cpp

@@ -4,13 +4,38 @@
 
 void EllipticPDESolver::run(IAssemble *assemble, SampleSolution &solution) {
     newton(*assemble, solution.U);
-    // Todo functional
-    // Todo plotting
+}
+
+void EllipticPDESolver::computeQ(IAssemble *assemble, SampleSolution &solution) {
+    IStochasticEllipticAssemble *ellAssemble =
+        dynamic_cast<IStochasticEllipticAssemble *>(assemble);
+    if (quantity == "L2") solution.Q = ellAssemble->L2(solution.U);
+    else if (quantity == "H1") solution.Q = ellAssemble->H1(solution.U);
+    else if (quantity == "Energy") solution.Q = ellAssemble->Energy(solution.U);
+    else if (quantity == "Inflow") solution.Q = ellAssemble->InflowOutflow(solution.U).first;
+    else if (quantity == "Outflow") solution.Q = ellAssemble->InflowOutflow(solution.U).second;
+    else Exit("Quantity of interest not implemented")
+}
+
+void EllipticPDESolver::computeCost(SampleSolution &solution) {
+    if (costMeasure == "size") solution.Cost = solution.U.size();
+//    else if (costMeasure == "time") solution.Cost = solution.U.size(); // Todo
+    else Exit("Cost measure not implemented")
+}
+
+void EllipticPDESolver::plotSolution(SampleSolution &solution) {
 }
 
 void TransportPDESolver::run(IAssemble *assemble, SampleSolution &solution) {
     TimeSeries timeSeries;
     timeInt(dynamic_cast<TAssemble *>(assemble), timeSeries, solution.U);
-    // Todo functional
-    // Todo plotting
+}
+
+void TransportPDESolver::computeQ(IAssemble *assemble, SampleSolution &solution) {
+}
+
+void TransportPDESolver::computeCost(SampleSolution &solution) {
+}
+
+void TransportPDESolver::plotSolution(SampleSolution &solution) {
 }

mlmc/src/montecarlo/PDESolver.hpp

@@ -11,17 +11,32 @@ class PDESolver {
 protected:
     int verbose = 1;
 
+    string quantity = "L2";
+
+    string costMeasure = "size";
+
     virtual void run(IAssemble *assemble, SampleSolution &solution) = 0;
 
+    virtual void computeQ(IAssemble *assemble, SampleSolution &solution) = 0;
+
+    virtual void computeCost(SampleSolution &solution) = 0;
+
+    virtual void plotSolution(SampleSolution &solution) = 0;
+
 public:
     PDESolver() {
         config.get("PDESolverVerbose", verbose);
+        config.get("Quantity", quantity);
+        config.get("CostMeasure", costMeasure);
     }
 
     void Run(IAssemble *assemble, SampleSolution &solution) {
         mout.StartBlock("PDESolver");
         vout(1) << solution.id.Str() << endl;
         run(assemble, solution);
+        computeQ(assemble, solution);
+        computeCost(solution);
+        plotSolution(solution);
         vout(1) << "Q=" << solution.Q << " cost=" << solution.Cost << endl;
         mout.EndBlock(verbose == 0);
     }
@@ -36,6 +51,12 @@ private:
 protected:
     void run(IAssemble *assemble, SampleSolution &solution) override;
 
+    void computeQ(IAssemble *assemble, SampleSolution &solution) override;
+
+    void computeCost(SampleSolution &solution) override;
+
+    void plotSolution(SampleSolution &solution) override;
+
 public:
     EllipticPDESolver() :
         PDESolver(),
@@ -52,11 +73,17 @@ private:
 protected:
     void run(IAssemble *assemble, SampleSolution &solution) override;
 
+    void computeQ(IAssemble *assemble, SampleSolution &solution) override;
+
+    void computeCost(SampleSolution &solution) override;
+
+    void plotSolution(SampleSolution &solution) override;
+
+public:
     TransportPDESolver() :
         PDESolver(),
         solver(Solver(GetPC("PointBlockGaussSeidel"), "GMRES")),
         timeInt(TimeIntegrator(solver)) {}
 };
 
-
 #endif //PDESOLVER_HPP