Commit 8a849b2b authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

worked on transport

parent 05fa842c
......@@ -3,8 +3,35 @@
using namespace std;
void MonteCarloTransport::method() {
Date Start;
vout(1) << " Start Monte Carlo method on level " << l;
if (!baseLevel) vout(1) << " and level " << l - 1;
vout(1) << " with " << dM << " samples" << endl;
stochFields->setCurrentLevel(l);
for (int m = M; m < M + dM; m++) {
stochFields->generate();
}
}
void MonteCarloTransport::solveTransport(double &Q, double &cost, Vector &u) {
Solver solver;
Preconditioner *BJ = GetPC("PointBlockJacobi");
Solver IM(BJ);
DGTimeIntegrator timeIntegrator(solver, solver, solver, IM);
TimeSeries timeSeries;
// transportAssemble.PrintInfo(u);
timeIntegrator(*assemble, timeSeries, u);
}
void MonteCarloTransport::solveElliptic(double &Q, double &cost, Vector &u) {
Preconditioner *preconditioner = GetPC(*graphs, *assemble, "LIB_PS");
Solver solver(preconditioner, "GMRES");
Newton newton(solver);
newton(*assemble, u);
}
\ No newline at end of file
......@@ -3,6 +3,7 @@
#include "MonteCarlo.h"
#include "assemble/DGTransportAssemble.h"
#include "TimeIntegrator.h"
class MonteCarloTransport : public MonteCarlo {
public:
......@@ -13,8 +14,8 @@ public:
StochasticFields *stochFields, MatrixGraphs *graphs,
DGTransportAssemble *assemble) :
assemble(assemble),
MonteCarlo(l, dM, baseLevel, meshes, stochFields, graphs) {
MonteCarlo(l, dM, baseLevel, meshes, stochFields, graphs),
assemble(assemble) {
transfer = new MatrixTransfer(*assemble);
transfer->Construct((*graphs)[l - plevel], (*graphs)[l - plevel - 1]);
......@@ -26,18 +27,27 @@ public:
}
void upscale(const Vector &uc, Vector &uc_up) override {
uc_up = (*transfer) * uc;
}
void updateUSum(const Vector &uc_up, const Vector &uf) {
u_sum += uf;
u_sum_diff += uf;
u_sum_diff -= uc_up;
}
void updateUAvg() {
Scalar M_inverted = 1.0 / M;
u_avg = M_inverted * u_sum;
u_avg_diff = M_inverted * u_sum_diff;
}
void method() override;
void solveTransport(double &Q, double &cost, Vector &u);
void solveElliptic(double &Q, double &cost, Vector &u);
};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment