@@ -55,216 +55,6 @@ namespace pfasst {
5555 virtual Encapsulation<scalar,time>* create (const EncapType) = 0;
5656 };
5757
58- template <typename scalar, typename time>
59- class EncapSweeper : public ISweeper {
60- vector<scalar> nodes;
61- shared_ptr<EncapFactory<scalar,time>> factory;
62-
63- public:
64-
65- void set_nodes (vector<time> nodes)
66- {
67- this ->nodes = nodes;
68- }
69-
70- const vector<time> get_nodes () const
71- {
72- return nodes;
73- }
74-
75- void set_factory (EncapFactory<scalar,time>* factory)
76- {
77- this ->factory = shared_ptr<EncapFactory<scalar,time>>(factory);
78- }
79-
80- EncapFactory<scalar,time>* get_factory () const
81- {
82- return factory.get ();
83- }
84-
85- virtual void set_state (const Encapsulation<scalar,time>* q0, unsigned int m)
86- {
87- throw NotImplementedYet (" sweeper" );
88- }
89-
90- virtual Encapsulation<scalar,time>* get_state (unsigned int m) const
91- {
92- throw NotImplementedYet (" sweeper" );
93- return NULL ;
94- }
95-
96- virtual Encapsulation<scalar,time>* get_tau (unsigned int m) const
97- {
98- throw NotImplementedYet (" sweeper" );
99- return NULL ;
100- }
101-
102- virtual Encapsulation<scalar,time>* get_saved_state (unsigned int m) const
103- {
104- throw NotImplementedYet (" sweeper" );
105- return NULL ;
106- }
107-
108- virtual Encapsulation<scalar,time>* get_end_state ()
109- {
110- return this ->get_state (this ->get_nodes ().size ()-1 );
111- }
112-
113- virtual void evaluate (int m)
114- {
115- throw NotImplementedYet (" sweeper" );
116- }
117-
118- virtual void advance ()
119- {
120- throw NotImplementedYet (" sweeper" );
121- }
122-
123- virtual void integrate (time dt, vector<Encapsulation<scalar,time>*> dst) const
124- {
125- throw NotImplementedYet (" sweeper" );
126- }
127- };
128-
129- template <typename scalar, typename time>
130- class PolyInterpMixin : public pfasst ::ITransfer {
131- matrix<time> tmat, fmat;
132-
133- public:
134-
135- virtual ~PolyInterpMixin () { }
136-
137- virtual void interpolate (ISweeper *dst, const ISweeper *src,
138- bool interp_delta_from_initial,
139- bool interp_initial)
140- {
141- auto * fine = dynamic_cast <EncapSweeper<scalar,time>*>(dst);
142- auto * crse = dynamic_cast <const EncapSweeper<scalar,time>*>(src);
143-
144- if (tmat.size1 () == 0 )
145- tmat = pfasst::compute_interp<time>(fine->get_nodes (), crse->get_nodes ());
146-
147- int nfine = fine->get_nodes ().size ();
148- int ncrse = crse->get_nodes ().size ();
149-
150- auto * crse_factory = crse->get_factory ();
151- auto * fine_factory = fine->get_factory ();
152-
153- vector<Encapsulation<scalar,time>*> fine_state (nfine), fine_delta (ncrse);
154-
155- for (int m=0 ; m<nfine; m++) fine_state[m] = fine->get_state (m);
156- for (int m=0 ; m<ncrse; m++) fine_delta[m] = fine_factory->create (solution);
157-
158- if (interp_delta_from_initial)
159- for (int m=1 ; m<nfine; m++)
160- fine_state[m]->copy (fine_state[0 ]);
161-
162- auto * crse_delta = crse_factory->create (solution);
163- int m0 = interp_initial ? 0 : 1 ;
164- for (int m=m0; m<ncrse; m++) {
165- crse_delta->copy (crse->get_state (m));
166- if (interp_delta_from_initial)
167- crse_delta->saxpy (-1.0 , crse->get_state (0 ));
168- else
169- crse_delta->saxpy (-1.0 , crse->get_saved_state (m));
170- interpolate (fine_delta[m], crse_delta);
171- }
172- delete crse_delta;
173-
174- if (! interp_initial)
175- fine_delta[0 ]->setval (0.0 );
176-
177- fine->get_state (0 )->mat_apply (fine_state, 1.0 , tmat, fine_delta, false );
178-
179- for (int m=0 ; m<ncrse; m++) delete fine_delta[m];
180- for (int m=m0; m<nfine; m++) fine->evaluate (m);
181- }
182-
183- virtual void restrict (ISweeper *dst, const ISweeper *src, bool restrict_initial)
184- {
185- auto * crse = dynamic_cast <EncapSweeper<scalar,time>*>(dst);
186- auto * fine = dynamic_cast <const EncapSweeper<scalar,time>*>(src);
187-
188- auto dnodes = crse->get_nodes ();
189- auto snodes = fine->get_nodes ();
190-
191- int ncrse = crse->get_nodes ().size ();
192- int nfine = fine->get_nodes ().size ();
193-
194- int trat = (nfine - 1 ) / (ncrse - 1 );
195-
196- int m0 = restrict_initial ? 0 : 1 ;
197- for (int m=m0; m<ncrse; m++) {
198- if (dnodes[m] != snodes[m*trat])
199- throw NotImplementedYet (" coarse nodes must be nested" );
200- this ->restrict (crse->get_state (m), fine->get_state (m*trat));
201- }
202-
203- for (int m=m0; m<ncrse; m++) crse->evaluate (m);
204- }
205-
206- virtual void fas (time dt, ISweeper *dst, const ISweeper *src)
207- {
208- auto * crse = dynamic_cast <EncapSweeper<scalar,time>*>(dst);
209- auto * fine = dynamic_cast <const EncapSweeper<scalar,time>*>(src);
210-
211- int ncrse = crse->get_nodes ().size ();
212- int nfine = fine->get_nodes ().size ();
213-
214- auto * crse_factory = crse->get_factory ();
215- auto * fine_factory = fine->get_factory ();
216-
217- vector<Encapsulation<scalar,time>*> crse_z2n (ncrse-1 ), fine_z2n (nfine-1 ), rstr_z2n (ncrse-1 );
218- for (int m=0 ; m<ncrse-1 ; m++) crse_z2n[m] = crse_factory->create (solution);
219- for (int m=0 ; m<ncrse-1 ; m++) rstr_z2n[m] = crse_factory->create (solution);
220- for (int m=0 ; m<nfine-1 ; m++) fine_z2n[m] = fine_factory->create (solution);
221-
222- // compute '0 to node' integral on the coarse level
223- crse->integrate (dt, crse_z2n);
224- for (int m=1 ; m<ncrse-1 ; m++)
225- crse_z2n[m]->saxpy (1.0 , crse_z2n[m-1 ]);
226-
227- // compute '0 to node' integral on the fine level
228- fine->integrate (dt, fine_z2n);
229- for (int m=1 ; m<nfine-1 ; m++)
230- fine_z2n[m]->saxpy (1.0 , fine_z2n[m-1 ]);
231-
232- // restrict '0 to node' fine integral
233- int trat = (nfine - 1 ) / (ncrse - 1 );
234- for (int m=1 ; m<ncrse; m++)
235- this ->restrict (rstr_z2n[m-1 ], fine_z2n[m*trat-1 ]);
236-
237- // compute 'node to node' tau correction
238- vector<Encapsulation<scalar,time>*> tau (ncrse-1 );
239- for (int m=0 ; m<ncrse-1 ; m++) tau[m] = crse->get_tau (m);
240-
241- tau[0 ]->copy (rstr_z2n[0 ]);
242- tau[0 ]->saxpy (-1.0 , crse_z2n[0 ]);
243-
244- for (int m=1 ; m<ncrse-1 ; m++) {
245- tau[m]->copy (rstr_z2n[m]);
246- tau[m]->saxpy (-1.0 , rstr_z2n[m-1 ]);
247-
248- tau[m]->saxpy (-1.0 , crse_z2n[m]);
249- tau[m]->saxpy (1.0 , crse_z2n[m-1 ]);
250- }
251-
252- for (int m=0 ; m<ncrse-1 ; m++) delete crse_z2n[m];
253- for (int m=0 ; m<ncrse-1 ; m++) delete rstr_z2n[m];
254- for (int m=0 ; m<nfine-1 ; m++) delete fine_z2n[m];
255- }
256-
257- // required for interp/restrict helpers
258- virtual void interpolate (Encapsulation<scalar,time> *dst, const Encapsulation<scalar,time> *src) {
259- throw NotImplementedYet (" mlsdc/pfasst" );
260- }
261-
262- virtual void restrict (Encapsulation<scalar,time> *dst, const Encapsulation<scalar,time> *src) {
263- throw NotImplementedYet (" mlsdc/pfasst" );
264- }
265-
266- };
267-
26858 }
26959}
27060
0 commit comments