Mercurial > hg > kwantix

comparison main.cpp @ 21:4c5bac1f2612

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Testing Euler's method
author Jordi Gutiérrez Hermoso <jordigh@gmail.com>
date Fri, 12 Sep 2008 11:13:18 -0500
parents 9cd42f1be76e
children 532145a38fa2
comparison
equal deleted inserted replaced
20:9cd42f1be76e 21:4c5bac1f2612
1 #include <iostream> 1 #include <iostream>
2 #include <iomanip>
2 #include <fstream> 3 #include <fstream>
3 #include <map> 4 #include <map>
4 #include <set> 5 #include <set>
5 #include <sstream> 6 #include <sstream>
6 7
16 #include "include/func.hpp" 17 #include "include/func.hpp"
17 18
18 using namespace linalg; 19 using namespace linalg;
19 using namespace bvp; 20 using namespace bvp;
20 21
21 #define RBF_TYPE rbf::thin_plate_spline 22 #define RBF_TYPE rbf::conical
22 23
23 //Will solve the system 24 //Will solve the system
24 // 25 //
25 // u_t = -u*u_x - v*u_y - g*h_x =: f1(u,v,h) 26 // u_t = -u*u_x - v*u_y - g*h_x =: f1(u,v,h)
26 // 27 //
32 //and h starts out as an undisturbed surface high with a small 33 //and h starts out as an undisturbed surface high with a small
33 //disturbance. We'll do this by RK4. 34 //disturbance. We'll do this by RK4.
34 35
35 36
36 const double g = 9.8; // m/s^2 37 const double g = 9.8; // m/s^2
38
39 const double pi = 3.141592653589793238462643383279502;
37 40
38 template<typename RBF> 41 template<typename RBF>
39 class Fgen { 42 class Fgen {
40 public: 43 public:
41 void set_interps(const interpolator<RBF>& u0, 44 void set_interps(const interpolator<RBF>& u0,
123 v_bvp_init(new linear_BVP2(Omega, I, D, Z, Z)); 126 v_bvp_init(new linear_BVP2(Omega, I, D, Z, Z));
124 127
125 128
126 //init the interps. 129 //init the interps.
127 using namespace rbf; 130 using namespace rbf;
128 conical::set_n(5); thin_plate_spline::set_n(6); 131 conical::set_n(7); thin_plate_spline::set_n(6);
129 c_infty_rbf::set_epsilon(10); 132 c_infty_rbf::set_epsilon(0.01);
130 133
131 cout << "Initialising... please wait." << endl; 134 cout << "Initialising... please wait." << endl;
132 135
133 interpolator<RBF_TYPE> 136 interpolator<RBF_TYPE>
134 u0(u_bvp_init), 137 u0(u_bvp_init),
141 144
142 //Intermediate interpolators for RK4 145 //Intermediate interpolators for RK4
143 std::vector<interpolator<RBF_TYPE> > 146 std::vector<interpolator<RBF_TYPE> >
144 // FIXME: make the copy ctor for interps make different copies 147 // FIXME: make the copy ctor for interps make different copies
145 // of the bvp. 148 // of the bvp.
146 k1(3,h0), k2(3,h0), k3(3,h0), k4(3,h0); 149 k1(3,u0), k2(3,u0), k3(3,u0), k4(3,u0);
150
151 //debug
152 cout << "Size of interpolator is " << sizeof(k1[0]) << endl;
147 153
148 Fu<RBF_TYPE> f_u; 154 Fu<RBF_TYPE> f_u;
149 Fv<RBF_TYPE> f_v; 155 Fv<RBF_TYPE> f_v;
150 Fh<RBF_TYPE> f_h; 156 Fh<RBF_TYPE> f_h;
151 157
152 //Negative because of the way the F's below are written. 158 //Negative because of the way the F's below are written.
153 double dt = -0.01; 159 double dt = -0.0001;
154 160
155 f_u.set_dt(dt); 161 f_u.set_dt(dt);
156 f_v.set_dt(dt); 162 f_v.set_dt(dt);
157 f_h.set_dt(dt); 163 f_h.set_dt(dt);
158 164
159 //main loop, timestepping with RK4. 165 //main loop, timestepping with RK4.
160 size_t maxiter = 1000; 166 size_t maxiter = 5000;
161 for(size_t i = 1; i <= maxiter; i++){ 167 for(size_t i = 1; i <= maxiter; i++){
162 cout << "Now on iteration #" << i << endl; 168 cout << "Now on iteration #" << i << endl;
163 if(i % 1 == 0){ 169 if(i % 1 == 0){
164 save_interp(Omega,u0,i,"u"); 170 save_interp(Omega,u0,i,"u");
165 save_interp(Omega,v0,i,"v"); 171 save_interp(Omega,v0,i,"v");
166 save_interp(Omega,h0,i,"h"); 172 save_interp(Omega,h0,i,"h");
167 } 173 }
168 174
169 //k1 175 //k1
170 f_u.set_interps(u0,v0,h0); 176 f_u.set_interps(u0,v0,h0);
171 f_v.set_interps(u0,v0,h0); 177 f_v.set_interps(u0,v0,h0);
172 f_h.set_interps(u0,v0,h0); 178 f_h.set_interps(u0,v0,h0);
173 179
174 k1[0].interpolate(f_u.at()); 180 // k1[0].interpolate(f_u.at());
175 k1[1].interpolate(f_v.at()); 181 // k1[1].interpolate(f_v.at());
176 k1[2].interpolate(f_h.at()); 182 // k1[2].interpolate(f_h.at());
183
184 //debug Euler
185 k1[0].interpolate(u0.at() + f_u.at());
186 k1[1].interpolate(v0.at() + f_v.at());
187 k1[2].interpolate(h0.at() + f_h.at());
177 188
178 set_bdry_conds(k1[0],k1[1]); 189 set_bdry_conds(k1[0],k1[1]);
179 190
180 cout << "k1 done!" << endl; 191 cout << "k1 done!" << endl;
192
193 //debug Euler
194 u0 = k1[0]; v0 = k1[1]; h0 = k1[2]; continue;
181 195
182 //k2 196 //k2
183 f_u.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2)); 197 f_u.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
184 f_v.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2)); 198 f_v.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
185 f_h.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2)); 199 f_h.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
296 const interpolator<RBF>& u, 310 const interpolator<RBF>& u,
297 size_t n, string which_interp) 311 size_t n, string which_interp)
298 { 312 {
299 using namespace std; 313 using namespace std;
300 string filename = "results/" + which_interp; 314 string filename = "results/" + which_interp;
301 if(n < 10000)
302 filename += "0";
303 if(n < 1000)
304 filename += "0";
305 if(n < 100)
306 filename += "0";
307 if(n < 10)
308 filename += "0";
309 stringstream ss; 315 stringstream ss;
316 ss << std::setw(5) << setfill('0') << n;
310 string n_string; 317 string n_string;
311 ss << n;
312 ss >> n_string; 318 ss >> n_string;
313 filename += n_string + ".map"; 319 filename += n_string + ".map";
314 ofstream ofs(filename.c_str()); 320 ofstream ofs(filename.c_str());
315 321 {
316 slice s(1,2); 322 slice s(1,2);
317 matrix M(Omega -> get_interior().size() + 323 matrix M(Omega -> get_interior().size() +
318 Omega -> get_boundary().size(), 3); 324 Omega -> get_boundary().size(), 3);
319 size_t k = 1; 325 size_t k = 1;
320 for(set<point>::const_iterator i = Omega -> get_interior().begin(); 326 for(set<point>::const_iterator i = Omega -> get_interior().begin();
321 i != Omega -> get_interior().end(); i++){ 327 i != Omega -> get_interior().end(); i++){
322 M(k,s) = *i; 328 M(k,s) = *i;
323 M(k,3) = u(*i); 329 M(k,3) = u(*i);
324 k++; 330 k++;
331 }
332 for(set<point>::const_iterator i = Omega -> get_boundary().begin();
333 i != Omega -> get_boundary().end(); i++){
334 M(k,s) = *i;
335 M(k,3) = u(*i);
336 k++;
337 }
338 ofs << M;
325 } 339 }
326 for(set<point>::const_iterator i = Omega -> get_boundary().begin(); 340
327 i != Omega -> get_boundary().end(); i++){ 341
328 M(k,s) = *i; 342 if(false){ //debug
329 M(k,3) = u(*i); 343 double r, th, N, M;
330 k++; 344 N = 30; M = 70;
331 } 345 linalg::matrix out(static_cast<size_t>(N),static_cast<size_t>(M));
332 ofs << M; 346 size_t i,j;
347 for(r = 0, j=1; j <= N; r += 1/(N-1),j++){
348 for(th = 0,i=1; i <= M; th += 2*pi/(M-1),i++){
349 linalg::vector p(2);
350 p(1) = r*cos(th);
351 p(2) = r*sin(th);
352 out(j,i) = u(p);
353 }
354 }
355
356 ofs << out;
357 }//debug
333 } 358 }
334 359
335 360
336 template<typename RBF> 361 template<typename RBF>
337 void set_bdry_conds(interpolator<RBF> & u, interpolator<RBF> & v) 362 void set_bdry_conds(interpolator<RBF> & u, interpolator<RBF> & v)