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;
}