Rutina principal y auxiliares (comprimido)

 


 

 

*
Source:    cd_imp_mq.cpp
                  Conveccion - Difusion en 2-d 
                  Metodo implicito 
                  Nucleo Multicuadrico
Compile:
                   c++ -o cdimq cd_imp_mq.cpp -lm 

Example:
                    ./cdimq 30 1 0.1 10 1  0.2 0.001 0.1

Output:
-----------------------------------------------------
Convection-Difussion - Implicit - MQ
n           = 30 , total(nxn) =900
a           = 1.000000
b           = 0.100000
convectivo  = 10.000000
difusivo    = 1.000000
c           = 0.200000
stept       = 0.001000
tmax        = 0.100000

t=0.100000, e_rms       = 0.000314
t=0.100000, e_inf       = 0.001146
-----------------------------------------------------

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]*phimq2(x[i]-x[j],y[i]-y[j],c,0);
   }
   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];
   }         
  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;
 double  u,ux,uy,uxx,uyy;
 
 double cond;
 double  e_rms,e_inf;
 double  c;
 double xmin,xmax,ymin,ymax;
 mesh   malla;   



 if(argc!=9){
   printf("ERROR:  requiere  <n> <a> <b> <convect> <diffus> <c> <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]);  
c  = atof(argv[6]);  
stept = atof(argv[7]);  
tmax  = atof(argv[8]);  

 tini = 0.0;
 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("c           = %f\n",c);      
 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,N);
  M           = make_dmatrix(N,N);
  P           = make_dmatrix(N,N);
  f           = make_dvector(N);
  lambda      = make_dvector(N);
  pivot_index = make_ivector(N);
  
  sol_analitica = make_dvector(malla.N);
  sol_numerica  = make_dvector(malla.N);
  sol_analitica_numerica = make_dvector(malla.N);   
   
  
//Creo la matriz de interpolacion:  condicion inicial
for(int i=0;i<N;i++)
{
    for(int j=0;j<N;j++)
    {
        A[i][j] = phimq2(x[i]-x[j],y[i]-y[j],c,0);
    }
  f[i]=u_t0(x[i],y[i]); 
}

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

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

//Calculo el error cuadratico medio y el error maximo
//sobre el mismo grid
  for(int i=0;i<malla.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;i++)
   for(int j=0;j<N;j++)
   {
      u = phimq2(x[i]-x[j],y[i]-y[j],c,0);   
     ux = phimq2(x[i]-x[j],y[i]-y[j],c,1);   
     uxx = phimq2(x[i]-x[j],y[i]-y[j],c,2);   
     uy = phimq2(x[i]-x[j],y[i]-y[j],c,3);        
     uyy = phimq2(x[i]-x[j],y[i]-y[j],c,4);        

     M[i][j] = u - 0.5*stept*Convectivo*(ux+uy) - 0.5*stept*Difusivo*(uxx+uyy);

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

       P[i][j] = u + 0.5*stept*Convectivo*(ux+uy) + 0.5*stept*Difusivo*(uxx+uyy);
   } 

//Factorizo la matriz-M
   Factor(M,N,&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
       mult_matrix_vec(P,lambda,N,f);
    //Correcion
       corrige_bc(x,y,ni,N,t+stept,f);      
    //Solucion     
       copy(f,N,lambda);
       Solve(M,N,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<malla.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<malla.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 - MQ\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("c           = %f\n",c);   
 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);  
  free_ivector(pivot_index);
  free_dvector(lambda);  

  free_dmatrix(M,N);    
  free_dmatrix(P,N);  
  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 228 visitantes y ningun miembro en Línea