-
Notifications
You must be signed in to change notification settings - Fork 1
/
MultiIntegrator.h
185 lines (150 loc) · 4.31 KB
/
MultiIntegrator.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
#ifndef __MultiIntegrator_h__
#define __MultiIntegrator_h__ 1
#include "kernel.h"
#include "Array.h"
using namespace Array;
extern unsigned Ngrids;
extern const char *subintegrator;
class MultiProblem : public ProblemBase {
protected:
DynVector<unsigned> gridindex;
public:
enum Field {Nfields};
unsigned Ngrids, grid, nfields;
vector3 mY;
vector y0; // saved data
virtual unsigned getNfields(unsigned g) {return (unsigned) nfields;};
void InitialConditions() {}
virtual void InitialConditions(unsigned Ngrids);
virtual unsigned getnfields(unsigned g)=0;
virtual void Project(unsigned toG)=0;
virtual void Prolong(unsigned toG)=0;
virtual int Rescale()=0;
virtual void doifnewy0() {};
};
void MultiProblem::InitialConditions(unsigned Ngrids0)
{
Ngrids=Ngrids0;
Allocate(mY,Ngrids);
}
class MultiIntegrator : public IntegratorBase {
public:
vector y0;
private:
protected:
bool new_y0;
Array1<RK *> Integrator;
Array1<DynVector<unsigned> > nY; //NB: ProblemBase has a variable NY.
MultiProblem *MProblem;
unsigned Nfields;
vector3 mY;
unsigned grid;
public:
const char *Name() {return "MultiIntegrator";}
virtual void Allocator(ProblemBase& problem,size_t Align=0);
void setGrid(unsigned g) {
grid=g;
MProblem->grid=g;
}
void TimestepDependence() {
for (unsigned g=0; g < Ngrids ; g++) {
setGrid(g);
Integrator[g]->SetTime(t,dt);
Integrator[g]->TimestepDependence();
}
}
Solve_RC Solve();
};
extern MultiProblem *MProblem;
void MultiIntegrator::Allocator(ProblemBase& problem, size_t Align)
{
const vector2& YP=problem.YVector();
DynVector<size_t>* NYP=problem.Sizes();
align=Align;
ny=0;
NY.SetDynVector(*NYP);
size_t nfields=YP.Size();
for(size_t i=0; i < nfields; i++)
ny += NY[i];
Dimension(Y,YP);
Dimension(y,ny,(Var *) (YP[0]));
Allocate(y0,ny,align);
Ngrids=::Ngrids;
if(Ngrids < 2) msg(ERROR,"Need more grids");
MProblem=::MProblem;
SetProblem(problem);
Allocate(Integrator,Ngrids);
Dimension(mY,MProblem->mY);
Set(MProblem->y0,y0);
Dimension(MProblem->y0,ny);
Allocate(nY,Ngrids);
for (unsigned g=0; g < Ngrids; g++) {
setGrid(g);
RK *integrator=dynamic_cast<RK *>(Vocabulary->NewIntegrator(subintegrator));
if(!integrator) msg(ERROR,"subintegrator must be an RK integrator");
Integrator[g]=integrator;
Integrator[g]->SetProblem(problem);
Integrator[g]->SetParam(*this);
Nfields=MProblem->getNfields(g);
mY[g].Allocate(Nfields);
for (unsigned F=0; F < Nfields; ++F) {
nY[g][F]=Problem->Size(Nfields*g+F);
Dimension(mY[g][F],nY[g][F],Problem->YVector()[Nfields*g+F]);
}
Integrator[g]->Allocator(mY[g],&nY[g],Problem->ErrorMask(),align);
}
new_y0=true;
// Assumes that sub-integrators are all the same order
// TODO: this could just take the lowest order and be fine.
order=Integrator[0]->Order();
pgrow=(order > 0) ? 0.5/order : 0;
pshrink=(order > 1) ? 0.5/(order-1) : pgrow;
// this should also give an option for rescaling.
// Add rescaling options to MultiIntegrator vocab or something?
}
Solve_RC MultiIntegrator::Solve() {
Solve_RC flag;
errmax=0.0;
if(new_y0)
for(unsigned i=0; i < ny; ++i) y0[i]=y[i]; // Save the initial data.
else
for(unsigned i=0; i < ny; ++i) y[i]=y0[i]; // Reload the initial data.
TimestepDependence();
unsigned lastgrid=Ngrids-1;
for (unsigned g=0; g <= lastgrid; g++) {
setGrid(g);
Integrator[g]->initialize0();
Integrator[g]->iSource();
Integrator[g]->Predictor(0,Integrator[g]->Ny());
if(Integrator[g]->Corrector(0,Integrator[g]->Ny())) {
double err=Integrator[g]->Errmax();
if(err > errmax) errmax=err;
flag=dynamic ? CheckError() : SUCCESSFUL;
new_y0=(flag != UNSUCCESSFUL);
} else {
flag=NONINVERTIBLE;
new_y0=false;
break;
}
if(new_y0) {
if(g < lastgrid)
MProblem->Project(g+1);
}
}
if(new_y0 && MProblem->Rescale() > 0) {
flag=NONINVERTIBLE;
new_y0=false;
}
if (new_y0) {
//cout << "newy0" << endl;
MProblem->doifnewy0();
for (unsigned g=lastgrid; g > 0; g--)
MProblem->Prolong(g-1);
for (unsigned g=0; g < Ngrids; g++) {
setGrid(g);
MProblem->Stochastic(mY[g],t,dt);
}
}
return flag;
}
#endif