Mercurial > hg > kwantix
changeset 21:4c5bac1f2612
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 |
| files | Makefile data/sw_eqs_init.m main-sw-rk4.cpp main.cpp results/plotresults.m results/plotresults_vtk.m |
| diffstat | 6 files changed, 214 insertions(+), 225 deletions(-) [+] |
line wrap: on
line diff
--- a/Makefile +++ b/Makefile @@ -1,9 +1,9 @@ CPP = g++ -LINKING = -lgsl -lgslcblas +LINKING = -lgsl -lgslcblas -static CFLAGS = -O3 OPTIONS = -Wall -pedantic -W -Werror -Wconversion -Wshadow \ -Wpointer-arith -Wcast-qual -Wcast-align -Wwrite-strings \ - -fshort-enums -fno-common -Wfatal-errors + -fshort-enums -fno-common -Wfatal-errors -static OBJECTFILES = interpolator.o linalg.o error.o utils.o rbf.o bvp.o \ diff_op.o ddm.o func.o main.o interp_values.o HEADERFILES = include/linalg.hpp include/error.hpp include/utils.hpp \
--- a/data/sw_eqs_init.m +++ b/data/sw_eqs_init.m @@ -1,5 +1,5 @@ -n = 15 -m = 12; +n = 20 +m = 6; r = linspace(0,1,n); rrintr = 0;
--- a/main-sw-rk4.cpp +++ b/main-sw-rk4.cpp @@ -3,7 +3,7 @@ #include <map> #include <set> #include <sstream> - + #include <boost/shared_ptr.hpp> #include "include/linalg.hpp" @@ -35,8 +35,10 @@ const double g = 9.8; // m/s^2 +const double pi = 3.141592653589793238462643383279502; + template<typename RBF> -class Fgen : public realfunc{ +class Fgen { public: void set_interps(const interpolator<RBF>& u0, const interpolator<RBF>& v0, @@ -59,7 +61,7 @@ using Fgen<RBF>::h; using Fgen<RBF>::dt; public: - double at(const point& p) const; + bvp::interp_values at() const; }; template<typename RBF> @@ -69,7 +71,7 @@ using Fgen<RBF>::h; using Fgen<RBF>::dt; public: - double at(const point& p) const; + bvp::interp_values at() const; }; template<typename RBF> @@ -79,32 +81,7 @@ using Fgen<RBF>::h; using Fgen<RBF>::dt; public: - double at(const point& p) const; -}; - -class iter_neumann : public bdry_diff_op{ -public: - //If this is u's, set 1, for v, set 2. - iter_neumann(size_t u_or_v_in); - double at(const realfunc &f, const point &p, const vector &n) const; -private: - size_t u_or_v; - double at(const realfunc &f, const point &p) const {return f(p);}; -}; - -template<typename RBF> -class Gu_or_v : public realfunc{ -public: - //specify the OTHER, whether u or v. - Gu_or_v(size_t u_or_v_in, - const shared_ptr<domain> Omega_in) : u_or_v(u_or_v_in), - Omega(Omega_in) {}; - double at(const point& p) const; - void set_other(const interpolator<RBF>& other_in){other = other_in;}; -private: - size_t u_or_v; - const shared_ptr<domain> Omega; - interpolator<RBF> other; + bvp::interp_values at() const; }; class zero_func : public realfunc{ @@ -120,8 +97,7 @@ size_t n, string which_interp); template<typename RBF> -void bdry_iter(interpolator<RBF> &u, interpolator<RBF> &v, - const boost::shared_ptr<domain> Omega); +void set_bdry_conds(interpolator<RBF> &u, interpolator<RBF> &v); //********************** Main ****************************************** @@ -132,7 +108,7 @@ try{ using namespace std; using boost::shared_ptr; - + map<point, double> h_init = utils::read_pd_map("data/h_init.map"); shared_ptr<domain> Omega(new domain("data/circ_intr.matrix", @@ -140,115 +116,119 @@ "data/circ_nrml.matrix")); shared_ptr<Id_op> I(new Id_op); - shared_ptr<iter_neumann> - u_bdry_op(new iter_neumann(1)), - v_bdry_op(new iter_neumann(2)); + shared_ptr<dirichlet_op> D(new dirichlet_op); zero_func Z; + //FIXME: Turns this into a purely interpolating interpolator, + //without a BVP. shared_ptr<linear_BVP2> - u_bvp_init(new linear_BVP2(Omega, I, u_bdry_op, Z, Z)), - v_bvp_init(new linear_BVP2(Omega, I, v_bdry_op, Z, Z)); + u_bvp_init(new linear_BVP2(Omega, I, D, Z, Z)), + v_bvp_init(new linear_BVP2(Omega, I, D, Z, Z)); - - double dt = 1e-2; - + //init the interps. using namespace rbf; - conical::set_n(5); thin_plate_spline::set_n(6); + conical::set_n(7); thin_plate_spline::set_n(6); c_infty_rbf::set_epsilon(0.01); cout << "Initialising... please wait." << endl; interpolator<RBF_TYPE> - u1, u0(u_bvp_init), - v1, v0(v_bvp_init), - h1, h0(Omega, h_init); - u1 = u0; v1 = v0; h1 = h0; + u0(u_bvp_init), + v0(v_bvp_init), + h0(Omega, h_init); + + //Precompute some RBFs + u0.precompute_ev(); v0.precompute_ev(), h0.precompute_ev(); + u0.precompute_d1(); v0.precompute_d1(), h0.precompute_d1(); //Intermediate interpolators for RK4 std::vector<interpolator<RBF_TYPE> > - k1(3,h0), k2(3,h0), k3(3,h0), k4(3,h0); + // FIXME: make the copy ctor for interps make different copies + // of the bvp. + k1(3,h0), k2(3,h0), k3(3,h0), k4(3,h0); Fu<RBF_TYPE> f_u; Fv<RBF_TYPE> f_v; Fh<RBF_TYPE> f_h; + //Negative because of the way the F's below are written. + double dt = -0.01; + f_u.set_dt(dt); f_v.set_dt(dt); f_h.set_dt(dt); - //main loop - size_t maxiter = 3; + //main loop, timestepping with RK4. + size_t maxiter = 5; for(size_t i = 1; i <= maxiter; i++){ cout << "Now on iteration #" << i << endl; if(i % 1 == 0){ - save_interp(Omega,u0,i,"u"); - save_interp(Omega,v0,i,"v"); - save_interp(Omega,h0,i,"h"); + save_interp(Omega,u0,i,"u"); + save_interp(Omega,v0,i,"v"); + save_interp(Omega,h0,i,"h"); } //k1 f_u.set_interps(u0,v0,h0); f_v.set_interps(u0,v0,h0); f_h.set_interps(u0,v0,h0); + + k1[0].interpolate(f_u.at()); + k1[1].interpolate(f_v.at()); + k1[2].interpolate(f_h.at()); + + set_bdry_conds(k1[0],k1[1]); - //debug - k1[0].set_f(f_u); - cout << "Finished k1[0]" << endl; - k1[1].set_f(f_v); - cout << "Finished k1[1]" << endl; - k1[2].set_f(f_h); - cout << "Finished k1[2]" << endl; + cout << "k1 done!" << endl; - bdry_iter(k1[0],k1[1],Omega); + //debug Euler + u0 = k1[0]; v0 = k1[1]; h0 = k1[2]; continue; //k2 f_u.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2)); f_v.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2)); f_h.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2)); + + k2[0].interpolate(f_u.at()); + k2[1].interpolate(f_v.at()); + k2[2].interpolate(f_h.at()); - //debug - k2[0].set_f(f_u); - cout << "Finished k2[0]" << endl; - k2[1].set_f(f_v); - cout << "Finished k2[1]" << endl; - k2[2].set_f(f_h); - cout << "Finished k2[2]" << endl; + set_bdry_conds(k2[0],k2[1]); - bdry_iter(k2[0],k2[1],Omega); + cout << "k2 done!" << endl; //k3 f_u.set_interps(u0+(k2[0]/2), v0+(k2[1]/2), h0+(k2[2]/2)); f_v.set_interps(u0+(k2[0]/2), v0+(k2[1]/2), h0+(k2[2]/2)); f_h.set_interps(u0+(k2[0]/2), v0+(k2[1]/2), h0+(k2[2]/2)); + + k3[0].interpolate(f_u.at()); + k3[1].interpolate(f_v.at()); + k3[2].interpolate(f_h.at()); - k3[0].set_f(f_u); - k3[1].set_f(f_v); - k3[2].set_f(f_h); - - bdry_iter(k3[0],k3[1],Omega); + set_bdry_conds(k3[0],k3[1]); + + cout << "k3 done!" << endl; //k4 f_u.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]); f_v.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]); f_h.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]); + + k4[0].interpolate(f_u.at()); + k4[1].interpolate(f_v.at()); + k4[2].interpolate(f_h.at()); - k4[0].set_f(f_u); - k4[1].set_f(f_v); - k4[2].set_f(f_h); + set_bdry_conds(k4[0],k4[1]); - bdry_iter(k4[0],k4[1],Omega); + cout << "k4 done!" << endl; //Grand finale - u1.set_f(u0 + (k1[0] + 2*k2[0] + 2*k3[0] + k4[0])/6); - v1.set_f(v0 + (k1[1] + 2*k2[1] + 2*k3[1] + k4[1])/6); - h1.set_f(h0 + (k1[2] + 2*k2[2] + 2*k3[2] + k4[2])/6); - - //Enforce boundary conditions iteratively - bdry_iter(u1,v1,Omega); - - u0 = u1; - v0 = v1; - h0 = h1; + u0.interpolate((u0 + (k1[0] + 2*k2[0] + 2*k3[0] + k4[0])/6).at()); + v0.interpolate((v0 + (k1[1] + 2*k2[1] + 2*k3[1] + k4[1])/6).at()); + h0.interpolate((h0 + (k1[2] + 2*k2[2] + 2*k3[2] + k4[2])/6).at()); + + set_bdry_conds(u0,v0); } return 0; @@ -297,54 +277,32 @@ template<typename RBF> -double Fu<RBF>::at(const point& p) const +interp_values Fu<RBF>::at() const { - return dt*(u(p)*u.d(p,1) + v(p)*u.d(p,2) + g*h.d(p,1)); -} - -template<typename RBF> -double Fv<RBF>::at(const point& p) const -{ - return dt*(u(p)*v.d(p,1) + v(p)*v.d(p,2) + g*h.d(p,2)); + return dt*(u()*u.d(1) + v()*u.d(2) + g*h.d(1)); } template<typename RBF> -double Fh<RBF>::at(const point& p) const -{ - return dt*(u(p)*h.d(p,1) + h(p)*u.d(p,1) + v(p)*h.d(p,2) + h(p)*v.d(p,2)); -} - -iter_neumann::iter_neumann(size_t u_or_v_in) +interp_values Fv<RBF>::at() const { - if(u_or_v_in < 1 or u_or_v_in > 2){ - error_handling::badArgument exc; - exc.reason = "Argument to iter_neumann constructor must be 1 or 2"; - exc.line = __LINE__; - exc.file = __FILE__; - throw exc; - } - u_or_v = u_or_v_in; -} - -double iter_neumann::at(const realfunc &f, const point &p, - const vector &n) const -{ - return n(u_or_v)*f(p); + return dt*(u()*v.d(1) + + v()*v.d(2) + + g*h.d(2)); } template<typename RBF> -double Gu_or_v<RBF>::at(const point& p) const +interp_values Fh<RBF>::at() const { - point n = Omega -> get_normals().at(p); - return -n(u_or_v)*other(p); + return dt*(u()*h.d(1) + h()*u.d(1) + v()*h.d(2) + h()*v.d(2)); } - - + template<typename RBF> void save_interp(const shared_ptr<domain> Omega, const interpolator<RBF>& u, size_t n, string which_interp) { + + using namespace std; string filename = "results/" + which_interp; if(n < 10000) @@ -361,63 +319,63 @@ ss >> n_string; filename += n_string + ".map"; ofstream ofs(filename.c_str()); + + { + slice s(1,2); + matrix M(Omega -> get_interior().size() + + Omega -> get_boundary().size(), 3); + size_t k = 1; + for(set<point>::const_iterator i = Omega -> get_interior().begin(); + i != Omega -> get_interior().end(); i++){ + M(k,s) = *i; + M(k,3) = u(*i); + k++; + } + for(set<point>::const_iterator i = Omega -> get_boundary().begin(); + i != Omega -> get_boundary().end(); i++){ + M(k,s) = *i; + M(k,3) = u(*i); + k++; + } + ofs << M; + } - slice s(1,2); - matrix M(Omega -> get_interior().size() + - Omega -> get_boundary().size(), 3); - size_t k = 1; - for(set<point>::const_iterator i = Omega -> get_interior().begin(); - i != Omega -> get_interior().end(); i++){ - M(k,s) = *i; - M(k,3) = u(*i); - k++; - } - for(set<point>::const_iterator i = Omega -> get_boundary().begin(); - i != Omega -> get_boundary().end(); i++){ - M(k,s) = *i; - M(k,3) = u(*i); - k++; - } - ofs << M; + + if(false){ //debug + double r, th, N, M; + N = 30; M = 70; + linalg::matrix out(static_cast<size_t>(N),static_cast<size_t>(M)); + size_t i,j; + for(r = 0, j=1; j <= N; r += 1/(N-1),j++){ + for(th = 0,i=1; i <= M; th += 2*pi/(M-1),i++){ + linalg::vector p(2); + p(1) = r*cos(th); + p(2) = r*sin(th); + out(j,i) = u(p); + } + } + + ofs << out; + }//debug } + template<typename RBF> -linalg::vector at_bdry(const interpolator<RBF>& u, - const boost::shared_ptr<domain> Omega) +void set_bdry_conds(interpolator<RBF> & u, interpolator<RBF> & v) { - using namespace linalg; - - vector out(Omega -> get_boundary().size()); - std::set<point>::const_iterator I = Omega -> get_boundary().begin(); - - for(size_t i = 1; i <= out.size(); i++){ - out(i) = u(*I); - I++; + bvp::bdry_values u_bdry(u.at()), v_bdry(v.at()); + bvp::normals N=u.get_normals(), M=v.get_normals(); + if(N != M){ + failure exc; + exc.reason="u and v have incompatible domains. This should never" + " happen."; + exc.line=__LINE__; exc.file=__FILE__; throw exc; } - return out; + //Project the vector (u,v) onto the perp of (N(1),N(2)) = (-N(2),N(1)). + bvp::bdry_values u_bdry_new = u_bdry*N(2)*N(2) - v_bdry*N(1)*N(2); + bvp::bdry_values v_bdry_new = v_bdry*N(1)*N(1) - u_bdry*N(1)*N(2); + + u.set_bdry_values(u_bdry_new); + v.set_bdry_values(v_bdry_new); } - -template<typename RBF> -void bdry_iter(interpolator<RBF> &u, interpolator<RBF> &v, - const boost::shared_ptr<domain> Omega) -{ - Gu_or_v<RBF_TYPE> gu(2, Omega), gv(1, Omega); - double err = 1; - do{ - linalg::vector u_old, u_new, v_old, v_new; - u_old = at_bdry(u,Omega); - gu.set_other(v); - u.set_g(gu); - u_new = at_bdry(u,Omega); - - v_old = at_bdry(v,Omega); - gv.set_other(u); - v.set_g(gv); - v_new = at_bdry(v,Omega); - - double err_u = norm(u_new - u_old); - double err_v = norm(v_new - v_old); - err = (err_u > err_v? err_u : err_v); - }while(err > 1e-2); -}
--- a/main.cpp +++ b/main.cpp @@ -1,4 +1,5 @@ #include <iostream> +#include <iomanip> #include <fstream> #include <map> #include <set> @@ -18,7 +19,7 @@ using namespace linalg; using namespace bvp; -#define RBF_TYPE rbf::thin_plate_spline +#define RBF_TYPE rbf::conical //Will solve the system // @@ -35,6 +36,8 @@ const double g = 9.8; // m/s^2 +const double pi = 3.141592653589793238462643383279502; + template<typename RBF> class Fgen { public: @@ -125,8 +128,8 @@ //init the interps. using namespace rbf; - conical::set_n(5); thin_plate_spline::set_n(6); - c_infty_rbf::set_epsilon(10); + conical::set_n(7); thin_plate_spline::set_n(6); + c_infty_rbf::set_epsilon(0.01); cout << "Initialising... please wait." << endl; @@ -143,27 +146,30 @@ std::vector<interpolator<RBF_TYPE> > // FIXME: make the copy ctor for interps make different copies // of the bvp. - k1(3,h0), k2(3,h0), k3(3,h0), k4(3,h0); + k1(3,u0), k2(3,u0), k3(3,u0), k4(3,u0); + + //debug + cout << "Size of interpolator is " << sizeof(k1[0]) << endl; Fu<RBF_TYPE> f_u; Fv<RBF_TYPE> f_v; Fh<RBF_TYPE> f_h; //Negative because of the way the F's below are written. - double dt = -0.01; + double dt = -0.0001; f_u.set_dt(dt); f_v.set_dt(dt); f_h.set_dt(dt); //main loop, timestepping with RK4. - size_t maxiter = 1000; + size_t maxiter = 5000; for(size_t i = 1; i <= maxiter; i++){ cout << "Now on iteration #" << i << endl; if(i % 1 == 0){ - save_interp(Omega,u0,i,"u"); - save_interp(Omega,v0,i,"v"); - save_interp(Omega,h0,i,"h"); + save_interp(Omega,u0,i,"u"); + save_interp(Omega,v0,i,"v"); + save_interp(Omega,h0,i,"h"); } //k1 @@ -171,13 +177,21 @@ f_v.set_interps(u0,v0,h0); f_h.set_interps(u0,v0,h0); - k1[0].interpolate(f_u.at()); - k1[1].interpolate(f_v.at()); - k1[2].interpolate(f_h.at()); +// k1[0].interpolate(f_u.at()); +// k1[1].interpolate(f_v.at()); +// k1[2].interpolate(f_h.at()); + + //debug Euler + k1[0].interpolate(u0.at() + f_u.at()); + k1[1].interpolate(v0.at() + f_v.at()); + k1[2].interpolate(h0.at() + f_h.at()); set_bdry_conds(k1[0],k1[1]); cout << "k1 done!" << endl; + + //debug Euler + u0 = k1[0]; v0 = k1[1]; h0 = k1[2]; continue; //k2 f_u.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2)); @@ -298,38 +312,49 @@ { using namespace std; string filename = "results/" + which_interp; - if(n < 10000) - filename += "0"; - if(n < 1000) - filename += "0"; - if(n < 100) - filename += "0"; - if(n < 10) - filename += "0"; stringstream ss; + ss << std::setw(5) << setfill('0') << n; string n_string; - ss << n; ss >> n_string; filename += n_string + ".map"; - ofstream ofs(filename.c_str()); + ofstream ofs(filename.c_str()); + { + slice s(1,2); + matrix M(Omega -> get_interior().size() + + Omega -> get_boundary().size(), 3); + size_t k = 1; + for(set<point>::const_iterator i = Omega -> get_interior().begin(); + i != Omega -> get_interior().end(); i++){ + M(k,s) = *i; + M(k,3) = u(*i); + k++; + } + for(set<point>::const_iterator i = Omega -> get_boundary().begin(); + i != Omega -> get_boundary().end(); i++){ + M(k,s) = *i; + M(k,3) = u(*i); + k++; + } + ofs << M; + } - slice s(1,2); - matrix M(Omega -> get_interior().size() + - Omega -> get_boundary().size(), 3); - size_t k = 1; - for(set<point>::const_iterator i = Omega -> get_interior().begin(); - i != Omega -> get_interior().end(); i++){ - M(k,s) = *i; - M(k,3) = u(*i); - k++; - } - for(set<point>::const_iterator i = Omega -> get_boundary().begin(); - i != Omega -> get_boundary().end(); i++){ - M(k,s) = *i; - M(k,3) = u(*i); - k++; - } - ofs << M; + + if(false){ //debug + double r, th, N, M; + N = 30; M = 70; + linalg::matrix out(static_cast<size_t>(N),static_cast<size_t>(M)); + size_t i,j; + for(r = 0, j=1; j <= N; r += 1/(N-1),j++){ + for(th = 0,i=1; i <= M; th += 2*pi/(M-1),i++){ + linalg::vector p(2); + p(1) = r*cos(th); + p(2) = r*sin(th); + out(j,i) = u(p); + } + } + + ofs << out; + }//debug }
--- a/results/plotresults.m +++ b/results/plotresults.m @@ -1,9 +1,11 @@ -r = linspace(0,1,40); -th = linspace(-pi,pi,60); +r = linspace(0,1,80); +th = linspace(0,2*pi,60); [rr,thth] = meshgrid(r,th); [xx,yy] = pol2cart(thth,rr); -for i = 1:1:1000 +skip = 1; + +for i = skip:skip:1000 name = "h"; if(i < 10000); name = strcat(name, "0"); @@ -21,11 +23,14 @@ filename = strcat(name,".map"); load(filename); h = eval(name); - zz = griddata(h(:,1),h(:,2),h(:,3),xx,yy); - surf(xx,yy,zz); - xlim([-1,1]); - ylim([-1,1]); -## zlim([0,0.2]); - print(strcat(name,".png"),"-dpng"); + x = h(:,1); y = h(:,2); z = h(:,3); + [xi,yi,zi] = griddata(x,y,z,xx,yy); + zi = zi(:,1:end-1); xi = xi(:,1:end-1); yi = yi(:,1:end-1); + vtk_clear; + vtk_surf(xi,yi,zi); + vtk_xlim([-1,1]); + vtk_ylim([-1,1]); + vtk_zlim([0,0.2]); + vtk_print(strcat(name,".png"),"-dpng"); clear eval(name); endfor