*
Source:    cd_imp_tps.cpp
                  Conveccion - Difusion en 2-d 
                  Metodo implicito 
                  Nucleo de placa delgada + Pol de grado 6 
Compile:
                   c++ -o cditps cd_imp_tps.cpp -lm 

Example:
                   ./cditps 30 1 0.1 10 1   0.001 0.1 
Output:
-----------------------------------------------------
Convection-Difussion - Implicit - TPS
n           = 30 , total(nxn) =900
a           = 1.000000
b           = 0.100000
convectivo  = 10.000000
difusivo    = 1.000000
stept       = 0.001000
tmax        = 0.100000

t=0.100000, e_rms       = 0.000676
t=0.100000, e_inf       = 0.002434
-----------------------------------------------------

Los siguientes archivos de datos que
se pueden ver con GNUPLOT
          
            splot  'real'
            splot  'numerica'
            splot  'realnumerica'
            splot  'realt'
            splot  'numericat'
            splot  'realnumericat'

real   -> corresponde a la solucion analitica de la EDP en t=0
realt  -> corresponde a la solucion analitica de la EDP en t=tmax
numerica   -> corresponde a la solucion numerica de la EDP en t=0
numericat   -> corresponde a la solucion numerica de la EDP en t=tmax
realnumerica ->     analitica - numerica  en t=0
realnumericat ->     analitica - numerica  en t=tmax

Autor: Jose Antonio Muñoz Gómez
Supervisión: Pedro González Casanova


*/
#include<math.h>
#include<stdlib.h>
#include<stdio.h>

#define INTERIOR       0
#define FRONTERA_REAL  1

#include "../libs/jag_mem.cpp"
#include "../libs/jag_phi.cpp"
#include "../libs/jag_vector.cpp"
#include "../libs/jag_cond.cpp"
#include "../libs/qSort.cpp"
#include "../libs/jag_file.cpp"

#include "cd_edp_ut.cpp"


#include "../libs/varios1.cpp"



//--------------------------------------------------------------------
void reconstruye(double *x,double *y, double *lambda,double c,int N, double *sol_numerica)
{
int     i,j;
double xc,yc;
double s;

 for(i=0;i<N;i++)
 {
  s=0.0;
   for(j=0;j<N;j++)
   {
      s += lambda[j]*phitps42(x[i]-x[j],y[i]-y[j],0,0);
   }

   s = s + lambda[N+0];
   s = s + lambda[N+1]*x[i];
   s = s + lambda[N+2]*y[i];
   s = s + lambda[N+3]*x[i]*y[i];
   s = s + lambda[N+4]*x[i]*x[i];
   s = s + lambda[N+5]*y[i]*y[i];

   sol_numerica[i] = s;
 }
}
//--------------------------------------------------------------------
//Proyeccion
void mult_matrix_vec(double **P, double *lambda,int m,double *z)
{
 double s;    
  for(int i=0;i<m;i++)    
  {
      s = 0.0; 
   for(int j=0;j<m;j++)    
   {
     s = s + P[i][j]*lambda[j];
   }         
  s = s + P[i][m+0]*lambda[m+0];
  s = s + P[i][m+1]*lambda[m+1];
  s = s + P[i][m+2]*lambda[m+2];
  s = s + P[i][m+3]*lambda[m+3];
  s = s + P[i][m+4]*lambda[m+4];
  s = s + P[i][m+5]*lambda[m+5];

  z[i] = s; 
 } 
      
}
//--------------------------------------------------------------------
//Correcion
void corrige_bc(double *x,double *y, int ni, int m, double t,double *f)
{
  for(int i=ni;i<m;i++)        
     f[i] = u_t(x[i],y[i],t);
}
//--------------------------------------------------------------------

int main(int argc, char **argv)
{
 double  t,tini,tmax,stept;    
 int     n,N,ni;

 
 double *sol_numerica;
 double *sol_analitica;
 double *sol_analitica_numerica; 
 
 double **A,**M,**P,*f;
 double *lambda; 
 int     *pivot_index;
 double *x,*y; 
 int    *tipo;
 int     Pol; 
 double  u,ux,uy,uxx,uyy;
 
 double cond;
 double  e_rms,e_inf;
 double  c;
 double xmin,xmax,ymin,ymax;
 mesh   malla;   


 if(argc!=8){
   printf("ERROR:  requiere  <n> <a> <b> <convect> <diffus> <stept> <tmax>\n\n");
   exit(1);     
 }
 
    n = atoi(argv[1]);
a_cd = atof(argv[2]);  
b_cd = atof(argv[3]);  
Convectivo = atof(argv[4]);  
Difusivo       = atof(argv[5]);  
stept = atof(argv[6]);  
tmax  = atof(argv[7]);  

 tini = 0.0;
 Pol  = 6;
 printf("-----------------------------------------------------\n");
 printf("n           = %d , total(nxn) =%d\n",n,n*n);
 printf("a           = %f\n",a_cd);   
 printf("b           = %f\n",b_cd);   
 printf("convectivo  = %f\n",Convectivo);   
 printf("difusivo    = %f\n",Difusivo);   
 printf("stept       = %f\n",stept);   
 printf("tmax        = %f\n",tmax);    


//Creo el grid cartesiano
 xmin =  0;
 ymin =  0;
 xmax =  1;
 ymax =  1;

  malla  = make_mesh(xmin,xmax,ymin,ymax,n);
  
  N = malla.N;
  x = malla.x;
  y = malla.y;
tipo= malla.tipo;  
 ni = malla.ni; 
 
 
//Reservo la memoria dinamica
  A           = make_dmatrix(N+Pol,N+Pol);
  M           = make_dmatrix(N+Pol,N+Pol);
  P           = make_dmatrix(N+Pol,N+Pol);
  f           = make_dvector(N+Pol);
  lambda      = make_dvector(N+Pol);
  pivot_index = make_ivector(N+Pol);
  
  sol_analitica          = make_dvector(N);
  sol_numerica           = make_dvector(N);
  sol_analitica_numerica = make_dvector(N);   
   
  
//Creo la matriz de interpolacion:  condicion inicial

//Relleno de zeros
  fill(f,N+Pol,0.0);


 for(int i=0;i<(N+Pol);i++)
   fill(A[i],N+Pol,0.0);


//Creo la matriz de interpolacion  
for(int i=0;i<N;i++)
{
    for(int j=0;j<N;j++)
    {
        A[i][j] = phitps42(x[i]-x[j],y[i]-y[j],0,0);
    }
  f[i]=u_t0(x[i],y[i]); 

 A[i][N+0] = 1.0;
 A[i][N+1] = x[i];
 A[i][N+2] = y[i];
 A[i][N+3] = x[i]*y[i];
 A[i][N+4] = x[i]*x[i];
 A[i][N+5] = y[i]*y[i];

 A[N+0][i] = 1.0;
 A[N+1][i] = x[i];
 A[N+2][i] = y[i];
 A[N+3][i] = x[i]*y[i];
 A[N+4][i] = x[i]*x[i];
 A[N+5][i] = y[i]*y[i];

}

//Factorizo la matriz
   Factor(A,N+Pol,&cond,pivot_index);

   printf("num.condi   = %e\n",cond);  
 
//Resuelvo mediante LU
  copy(f,N+Pol,lambda);
  Solve(A,N+Pol,pivot_index,lambda);
 

//Calculo el error cuadratico medio y el error maximo
//sobre el mismo grid
  for(int i=0;i<N;i++)
    sol_analitica[i] = u_t0(x[i],y[i]);

  reconstruye(x,y,lambda,c,N,sol_numerica);

  e_rms = calc_erms(sol_analitica,sol_numerica,N);
  e_inf = calc_einf(sol_analitica,sol_numerica,N);
  
  for(int i=0;i<malla.N;i++) 
      sol_analitica_numerica[i] = sol_analitica[i]-sol_numerica[i];     
    
  save_gnu_solucion("numerica",x,y,sol_numerica,N);  
  save_gnu_solucion("real",x,y,sol_analitica,N);     
  save_gnu_solucion("realnumerica",x,y,sol_analitica_numerica,N);      


 printf("t=%f, e_rms       = %e\n",0.0,e_rms);
 printf("t=%f, e_inf       = %e\n",0.0,e_inf); 




//Creo la matriz-M  
 for(int i=0;i<(N+Pol);i++)
   fill(M[i],N+Pol,0.0);

  for(int i=0;i<ni;i++)
  {
   for(int j=0;j<N;j++)
   {
      u = phitps42(x[i]-x[j],y[i]-y[j],c,0);   
     ux    = phitps42(x[i]-x[j],y[i]-y[j],c,1);   
     uxx  = phitps42(x[i]-x[j],y[i]-y[j],c,2);   
     uy    = phitps42(x[i]-x[j],y[i]-y[j],c,3);        
     uyy    = phitps42(x[i]-x[j],y[i]-y[j],c,4);        

     M[i][j] = u  -   0.5*stept*Difusivo*(uxx+uyy) - 0.5*stept*Convectivo*(ux+uy);
   } 
    M[i][N+0] = 1.0;
    M[i][N+1] = x[i];
    M[i][N+2] = y[i];
    M[i][N+3] = x[i]*y[i];
    M[i][N+4] = x[i]*x[i];
    M[i][N+5] = y[i]*y[i];

    M[i][N+0] -= 0.5*stept*Difusivo*(0.0);
    M[i][N+1] -= 0.5*stept*Difusivo*(0.0);
    M[i][N+2] -= 0.5*stept*Difusivo*(0.0);
    M[i][N+3] -= 0.5*stept*Difusivo*(0.0);
    M[i][N+4] -= 0.5*stept*Difusivo*(2.0);
    M[i][N+5] -= 0.5*stept*Difusivo*(2.0);

    M[i][N+0] -= 0.5*stept*Convectivo*(0.0);
    M[i][N+1] -= 0.5*stept*Convectivo*(1.0);
    M[i][N+2] -= 0.5*stept*Convectivo*(1.0);
    M[i][N+3] -= 0.5*stept*Convectivo*(x[i]+y[i]);
    M[i][N+4] -= 0.5*stept*Convectivo*(2.0*x[i]);
    M[i][N+5] -= 0.5*stept*Convectivo*(2.0*y[i]);

    M[N+0][i] = 1.0;
    M[N+1][i] = x[i];
    M[N+2][i] = y[i];
    M[N+3][i] = x[i]*y[i];
    M[N+4][i] = x[i]*x[i];
    M[N+5][i] = y[i]*y[i];
  }

  for(int i=ni;i<N;i++)
  {
   for(int j=0;j<N;j++)
   {
      u = phitps42(x[i]-x[j],y[i]-y[j],0,0);   
     M[i][j] = u;
   } 
    M[i][N+0] = 1.0;
    M[i][N+1] = x[i];
    M[i][N+2] = y[i];
    M[i][N+3] = x[i]*y[i];
    M[i][N+4] = x[i]*x[i];
    M[i][N+5] = y[i]*y[i];

    M[N+0][i] = 1.0;
    M[N+1][i] = x[i];
    M[N+2][i] = y[i];
    M[N+3][i] = x[i]*y[i];
    M[N+4][i] = x[i]*x[i];
    M[N+5][i] = y[i]*y[i];
  }


    
//Creo la matriz-P  
  for(int i=0;i<N;i++)
  {
   for(int j=0;j<N;j++)
   {
          u     = phitps42(x[i]-x[j],y[i]-y[j],c,0);   
     ux    = phitps42(x[i]-x[j],y[i]-y[j],c,1);   
     uxx  = phitps42(x[i]-x[j],y[i]-y[j],c,2);   
     uy    = phitps42(x[i]-x[j],y[i]-y[j],c,3);        
     uyy  = phitps42(x[i]-x[j],y[i]-y[j],c,4);        

     P[i][j] = u  +   0.5*stept*Difusivo*(uxx+uyy) + 0.5*stept*Convectivo*(ux+uy);
   } 
    P[i][N+0] = 1.0;
    P[i][N+1] = x[i];
    P[i][N+2] = y[i];
    P[i][N+3] = x[i]*y[i];
    P[i][N+4] = x[i]*x[i];
    P[i][N+5] = y[i]*y[i];

    P[i][N+0] += 0.5*stept*Difusivo*(0.0);
    P[i][N+1] += 0.5*stept*Difusivo*(0.0);
    P[i][N+2] += 0.5*stept*Difusivo*(0.0);
    P[i][N+3] += 0.5*stept*Difusivo*(0.0);
    P[i][N+4] += 0.5*stept*Difusivo*(2.0);
    P[i][N+5] += 0.5*stept*Difusivo*(2.0);

    P[i][N+0] += 0.5*stept*Convectivo*(0.0);
    P[i][N+1] += 0.5*stept*Convectivo*(1.0);
    P[i][N+2] += 0.5*stept*Convectivo*(1.0);
    P[i][N+3] += 0.5*stept*Convectivo*(x[i]+y[i]);
    P[i][N+4] += 0.5*stept*Convectivo*(2.0*x[i]);
    P[i][N+5] += 0.5*stept*Convectivo*(2.0*y[i]);
}

//Factorizo la matriz-M
   Factor(M,N+Pol,&cond,pivot_index);

   printf("num.condi   = %e\n",cond);  
   

//Itero sobre el tiempo   
  t = tini;
  
  int ask = int ((tmax / stept) / 10.0);
  int cont=0;

  while(t<tmax)
  {
    if(!((cont++)%ask)) 
       printf("%f\n",t);

    //Proyeccion
       fill(f,N+Pol,0.0);
       mult_matrix_vec(P,lambda,N,f);
       
    //Correcion
       corrige_bc(x,y,ni,N,t+stept,f);      
       
    //Solucion     
       copy(f,N+Pol,lambda);
       Solve(M,N+Pol,pivot_index,lambda);
            
    //Avanzo en el tiempo   
       t = t + stept; 
  }

  printf("\n");
  
//Calculo el error cuadratico medio y el error maximo
//sobre el mismo grid
  for(int i=0;i<N;i++)
    sol_analitica[i] = u_t(x[i],y[i],t);

  reconstruye(x,y,lambda,c,N,sol_numerica);

  e_rms = calc_erms(sol_analitica,sol_numerica,N);
  e_inf = calc_einf(sol_analitica,sol_numerica,N);
  
  for(int i=0;i<N;i++)     
     sol_analitica_numerica[i] = sol_analitica[i]-sol_numerica[i];     
    
  save_gnu_solucion("numericat",x,y,sol_numerica,N);  
  save_gnu_solucion("realt",x,y,sol_analitica,N);     
  save_gnu_solucion("realnumericat",x,y,sol_analitica_numerica,N);      

 printf("-----------------------------------------------------\n");     
 printf("Convection-Difussion - Implicit - TPS\n");
 printf("n           = %d , total(nxn) =%d\n",n,n*n);
 printf("a           = %f\n",a_cd);   
 printf("b           = %f\n",b_cd);   
 printf("convectivo  = %f\n",Convectivo);   
 printf("difusivo    = %f\n",Difusivo);   
 printf("stept       = %f\n",stept);   
 printf("tmax        = %f\n\n",tmax);    
 printf("t=%f, e_rms       = %f\n",t,e_rms);
 printf("t=%f, e_inf       = %f\n",t,e_inf);
 printf("-----------------------------------------------------\n");     
 
  free_dmatrix(A,N+Pol);  
  free_ivector(pivot_index);
  free_dvector(lambda);  

  free_dmatrix(M,N+Pol);    
  free_dmatrix(P,N+Pol);  
  free_dvector(f);  
   
 
//Libero la memoria dinamica
  free_dvector(x);
  free_dvector(y);
  free_ivector(tipo);
  free_dvector(sol_analitica);
  free_dvector(sol_numerica);
  free_dvector(sol_analitica_numerica);  
 
return 0;
}

Superior

Tenemos 18 visitantes y ningun miembro en Línea