1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950 |
- #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];
- }
- void
- HorizontalProblem::run(){
- if(Debug::level>0) cout<<" [HorizontalProblem::run] start"<<endl;
- /*for(size_t ix=0;ix<geometry->nX;++ix){
- hydr[ix]=previous_hydr[ix];
- }*/
- solve_system();
- error=error2(previous_hydr,hydr,geometry->nX);
- if(Debug::level>0) cout<<" [HorizontalProblem::run] end"<<endl;
- }
- 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);
- }
- }
|