Rutina principal y auxiliares (comprimido)

 


 

/*
Source:      ad_imp_mq.cpp
                  Pura Conveccion en 2-d  (una gaussiana desplazandose)
                  Metodo implicito 
                  Nucleo Multicuadrico
Compile:
                 c++ -o adimq  ad_imp_mq.cpp -lm 

Example:
                     ./adimq 30 0.01 0.1 0.01 0.4 

Output;
-----------------------------------------------------
Advection - Implicit - MQ
n           = 30 , total(nxn) =900
sigma       = 0.010000
c           = 0.100000
stept       = 0.010000
tmax        = 0.400000
t=0.400000, e_rms       = 0.004325
t=0.400000, e_inf       = 0.036591
-----------------------------------------------------
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 "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;
 
 double cond;
 double  e_rms,e_inf;
 double  c;
 double xmin,xmax,ymin,ymax;
 mesh   malla;   


 if(argc!=6){
   printf("ERROR:  requiere  <n>  <sigma> <c> <stept> <tmax>\n\n");
   exit(1);     
 }
 
    n = atoi(argv[1]);
SIGMA=atof(argv[2]);
    c = atof(argv[3]); 
stept = atof(argv[4]);  
tmax  = atof(argv[5]);  

 tini = 0.0;

 printf("-----------------------------------------------------\n");
 printf("n           = %d , total(nxn) =%d\n",n,n*n);
 printf("sigma       = %f\n",SIGMA);   
 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);   
     uy = phimq2(x[i]-x[j],y[i]-y[j],c,3);        
     M[i][j] = u + 0.5*stept*(ux+uy);

     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);   
     uy = phimq2(x[i]-x[j],y[i]-y[j],c,3);        
     P[i][j] = u - 0.5*stept*(ux+uy);
   } 

//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("Advection - Implicit - MQ\n");
 printf("n           = %d , total(nxn) =%d\n",n,n*n);
 printf("sigma       = %f\n",SIGMA);   
 printf("c           = %f\n",c);  
 printf("stept       = %f\n",stept);   
 printf("tmax        = %f\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 165 visitantes y ningun miembro en Línea