123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566 |
- #include "horizontal_problem.hpp"
- #include "debug.hpp"
- namespace Kernel{
- void
- HorizontalProblem::init(const Geometry* geometry_){
- geometry=geometry_;
- sub_M=new double[geometry->nX-1];
- diag_M=new double[geometry->nX];
- sup_M=new double[geometry->nX-1];
- F=new double[geometry->nX];
- hydr=new double[geometry->nX];
- }
- void
- HorizontalProblem::compute_error(){
- total_error=0;
- for(size_t ix=0;ix<geometry->nX;++ix){
- double e=abs(previous_hydr[ix]-hydr[ix])/n1_init_hydr;
- error_x[ix]=e;
- total_error+=e;
- }
- }
- void
- HorizontalProblem::run(){
- #ifdef MYDEBUG
- cout<<" [HorizontalProblem::run] start"<<endl;
- #endif
- /*for(size_t ix=0;ix<geometry->nX;++ix){
- hydr[ix]=previous_hydr[ix];
- }*/
- solve_system();
- compute_error();
- #ifdef MYDEBUG
- cout<<" [HorizontalProblem::run] end"<<endl;
- #endif
- }
- void
- HorizontalProblem::solve_system(){
- for(size_t ix=1;ix<geometry->nX-1;++ix){
- double dx=geometry->dX;
- double hbot=geometry->hbot[ix];
- double lx=l[ix];
- double a=-Physics::tilde_a(lx,hbot)*Physics::k0/(dx*dx);
- double c=Physics::tilde_c(lx,hbot)*Physics::k0/dx;
- double dhbot=geometry->dhbot[ix];
- sub_M[ix-1]=a-(dhbot*c)/2;
- diag_M[ix]=1-2*a;
- sup_M[ix]=a+(dhbot*c)/2;
- F[ix]=Pl[ix]/(Physics::rho*Physics::g)+lx;
- }
- diag_M[0]=1;
- sup_M[0]=-1;
- F[0]=0;
- diag_M[geometry->nX-1]=1;
- sub_M[geometry->nX-2]=-1;
- F[geometry->nX-1]=0;
- Thomas(geometry->nX,sub_M,diag_M,sup_M,F,hydr);
- }
- }
|