/*************************************************************************** * Copyright (C) 2009 by Douglas Castro * * douglas@ime.unicamp.br * * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program; if not, write to the * * Free Software Foundation, Inc., * * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * ***************************************************************************/ #include "timestep.h" double rk1ndv2(int d, int prof, VFl pc[], double Dt, double ex, double xmin[], double xmax[], Fluxo f) { double D = 0.0; int dir; double cvol = 1.0; for(dir = 0; dir < d; dir++) { double dx = fabs(xmax[dir] - xmin[dir]); cvol *= dx; } for(dir = 0; dir < d; dir++) { double dx = fabs(xmax[dir] - xmin[dir]); double Dx = /*ex == dir ? 0.5*(dx + 0.5*dx) :*/ dx; D += (pc[1].p->flxsup[dir] - pc[1].p->flxinf[dir])/Dx; //D += (dx*pc[1].p->flxsup[dir] - dx*pc[1].p->flxinf[dir])/cvol; } double s3 = pc[1].p->fval; /* solucao (nesta celula) antes do primeiro passo de rk */ return s3 - Dt*D; /* primeiro passo de rk */ } double rk2ndv2(int d, int prof, VFl pc[], double Dt, double ex, double xmin[], double xmax[], Fluxo f) { double D = 0.0; int dir; double cvol = 1.0; for(dir = 0; dir < d; dir++) { double dx = fabs(xmax[dir] - xmin[dir]); cvol *= dx; } for(dir = 0; dir < d; dir++) { double Dx = fabs(xmax[dir] - xmin[dir]); D += (pc[1].p->flxsup[dir] - pc[1].p->flxinf[dir])/Dx; } double s3 = pc[1].p->fval; /* solucao (nesta celula) depois de aplicar rk1_nd */ double s2 = pc[1].p->fv[0]; /* solucao (nesta celula) antes de aplicar rk1_nd */ return (s2 + s3 - Dt*D)/2.0; /* segundo passo de rk */ } void avancaTempoMeio(int d, int prof, VFl pc[], double xmin[], double xmax[], double Dt, Fluxo f) { if(pc[1].p == NULL){ return; } int ex = prof%d; int qual; for (qual = 0; qual <= 1; qual++) { //calcula a posicao da celula filho double xminf[d], xmaxf[d]; limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf); VFl pcf[d+1]; /* Constroi o pacote do filho {qual}: */ divide_pacote_menor(d, pc, ex, qual, pcf); avancaTempoMeio(d, prof + 1, pcf, xminf, xmaxf, Dt, f); } assert((pc[1].p->fil[0] == NULL) == (pc[1].p->fil[1] == NULL)); if(pc[1].p->fil[0] == NULL) { /* primeiro passo de Runge-Kutta 2 */ pc[1].p->fv[0] = rk1ndv2( d, prof, pc, Dt, ex, xmin, xmax, f ); /* */ } } void avancaTempoUm(int d, int prof, VFl pc[], double xmin[], double xmax[], double Dt, Fluxo f) { if(pc[1].p == NULL){return;} int ex = prof%d; int qual; for (qual = 0; qual <= 1; qual++) { //calcula a posicao da celula filho double xminf[d], xmaxf[d]; limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf); VFl pcf[d+1]; /* Constrói o pacote do filho {qual}: */ divide_pacote_menor(d, pc, ex, qual, pcf); avancaTempoUm(d, prof + 1, pcf, xminf, xmaxf, Dt, f); } assert((pc[1].p->fil[0] == NULL) == (pc[1].p->fil[1] == NULL)); if(pc[1].p->fil[0] == NULL) { /* segundo passo de Runge-Kutta 2 */ double tn2 = rk2ndv2( d, prof, pc, Dt, ex, xmin, xmax, f);/* */ pc[1].p->fv[1] = tn2; pc[1].p->fval = tn2; pc[1].p->fv[0] = tn2; } } /** ------------------ Runge-Kutta 3 ordem ------------------ */ double rk1nd3o(int d, No *p, double Dt, double xmin[], double xmax[]) { double D = 0.0; int dir; for(dir = 0; dir < d; dir++) { double Dx = fabs(xmax[dir] - xmin[dir]); D += (p->flxsup[dir] - p->flxinf[dir])/Dx; } double v0 = p->fval; // solucao (nesta celula) antes do 1 passo de rk return (v0 - Dt*D); // primeiro passo de rk } double rk2nd3o(int d, No *p, double Dt, double xmin[], double xmax[]) { double D = 0.0; int dir; for(dir = 0; dir < d; dir++) { double Dx = fabs(xmax[dir] - xmin[dir]); D += (p->flxsup[dir] - p->flxinf[dir])/Dx; } double v1 = p->fval; // solucao (nesta celula) depois de aplicar rk1nd3o double v0 = p->fv[0]; // solucao (nesta celula) antes de aplicar rk1nd3o return ((3.0/4.0)*v0 + (1.0/4.0)*v1 - Dt*(1.0/4.0)*D); // segundo passo de rk } double rk3nd3o(int d, No *p, double Dt, double xmin[], double xmax[]) { double D = 0.0; int dir; for(dir = 0; dir < d; dir++) { double Dx = fabs(xmax[dir] - xmin[dir]); D += (p->flxsup[dir] - p->flxinf[dir])/Dx; } double v2 = p->fval; // solucao (nesta celula) depois de aplicar rk2nd3o double v0 = p->fv[0]; // solucao (nesta celula) antes de aplicar rk1nd3o return ( (1.0/3.0)*v0 + (2.0/3.0)*v2 - Dt*(2.0/3.0)*D ); // 3 passo de rk } void avancaT13(int d, int prof, No *p, double xmin[], double xmax[], double Dt, Fluxo f) { if(p == NULL){ return; } int ex = prof%d; int qual; for (qual = 0; qual <= 1; qual++) { //calcula a posicao da celula filho double xminf[d], xmaxf[d]; limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf); No *pcf = p->fil[qual]; avancaT13(d, prof + 1, pcf, xminf, xmaxf, Dt, f); } assert((p->fil[0] == NULL) == (p->fil[1] == NULL)); if(p->fil[0] == NULL) { // primeiro passo de Runge-Kutta 3 p->fv[0] = rk1nd3o(d, p, Dt, xmin, xmax); } } void avancaT23(int d, int prof, No *p, double xmin[], double xmax[], double Dt, Fluxo f) { if(p == NULL){return;} int ex = prof%d; int qual; for (qual = 0; qual <= 1; qual++) { //calcula a posicao da celula filho double xminf[d], xmaxf[d]; limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf); No *pcf = p->fil[qual]; avancaT23(d, prof + 1, pcf, xminf, xmaxf, Dt, f); } assert((p->fil[0] == NULL) == (p->fil[1] == NULL)); if(p->fil[0] == NULL) { // segundo passo de Runge-Kutta 3 p->fv[1] = rk2nd3o( d, p, Dt, xmin, xmax); } } void avancaT33(int d, int prof, No *p, double xmin[], double xmax[], double Dt, Fluxo f) { if(p == NULL){return;} int ex = prof%d; int qual; for (qual = 0; qual <= 1; qual++) { //calcula a posicao da celula filho double xminf[d], xmaxf[d]; limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf); No *pcf = p->fil[qual]; avancaT33(d, prof + 1, pcf, xminf, xmaxf, Dt, f); } assert((p->fil[0] == NULL) == (p->fil[1] == NULL)); if(p->fil[0] == NULL) { // segundo passo de Runge-Kutta 2 double tn2 = rk2nd3o( d, p, Dt, xmin, xmax); //p->fv[1] = tn2; //p->fval = tn2; p->fv[2] = tn2; //p->fv[0] = tn2; } } /** ------------------ Runge-Kutta 3 ordem ------------------ */ /* double rk1nd3o(int d, No *p, double Dt, double xmin[], double xmax[]) { double D = 0.0; int dir; for(dir = 0; dir < d; dir++) { double Dx = fabs(xmax[dir] - xmin[dir]); D += (p->flxsup[dir] - p->flxinf[dir])/Dx; } double v0 = p->fval; // solucao (nesta celula) antes do primeiro passo de rk p->fv[0] = v0; return v0 - Dt*D; // primeiro passo de rk } double rk2nd3o(int d, No *p, double Dt, double xmin[], double xmax[]) { double D = 0.0; int dir; for(dir = 0; dir < d; dir++) { double Dx = fabs(xmax[dir] - xmin[dir]); D += (p->flxsup[dir] - p->flxinf[dir])/Dx; } double v1 = p->fval; // solucao (nesta celula) depois de aplicar rk1nd3o double v0 = p->fv[0]; // solucao (nesta celula) antes de aplicar rk1nd3o p->fv[1] = v1; return (0.75*v0 + 0.25*v1 - Dt*0.25*D); // segundo passo de rk } double rk3nd3o(int d, No *p, double Dt, double xmin[], double xmax[]) { double D = 0.0; int dir; for(dir = 0; dir < d; dir++) { double Dx = fabs(xmax[dir] - xmin[dir]); D += (p->flxsup[dir] - p->flxinf[dir])/Dx; } double v2 = p->fval; // solucao (nesta celula) depois de aplicar rk2nd3o double v0 = p->fv[0]; // solucao (nesta celula) antes de aplicar rk1nd3o return ( (1.0/3.0)*v0 + (2.0/3.0)*v2 - Dt*(2.0/3.0)*D ); // 3 passo de rk } void avancaT13(int d, int prof, No *p, double xmin[], double xmax[], double Dt, Fluxo f) { if(p == NULL){ return; } int ex = prof%d; int qual; for (qual = 0; qual <= 1; qual++) { //calcula a posicao da celula filho double xminf[d], xmaxf[d]; limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf); No *pcf = p->fil[qual]; // Constroi o pacote do filho {qual}: //divide_pacote_menor(d, pc, ex, qual, pcf); avancaT13(d, prof + 1, pcf, xminf, xmaxf, Dt, f); } assert((p->fil[0] == NULL) == (p->fil[1] == NULL)); if(p->fil[0] == NULL) { // primeiro passo de Runge-Kutta 3 p->fval = rk1nd3o(d, p, Dt, xmin, xmax); } } void avancaT23(int d, int prof, No *p, double xmin[], double xmax[], double Dt, Fluxo f) { if(p == NULL){return;} int ex = prof%d; int qual; for (qual = 0; qual <= 1; qual++) { //calcula a posicao da celula filho double xminf[d], xmaxf[d]; limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf); No *pcf = p->fil[qual]; // Constrói o pacote do filho {qual}: //divide_pacote_menor(d, pc, ex, qual, pcf); avancaT23(d, prof + 1, pcf, xminf, xmaxf, Dt, f); } assert((p->fil[0] == NULL) == (p->fil[1] == NULL)); if(p->fil[0] == NULL) { // segundo passo de Runge-Kutta 3 p->fval = rk2nd3o( d, p, Dt, xmin, xmax); } } void avancaT33(int d, int prof, No *p, double xmin[], double xmax[], double Dt, Fluxo f) { if(p == NULL){return;} int ex = prof%d; int qual; for (qual = 0; qual <= 1; qual++) { //calcula a posicao da celula filho double xminf[d], xmaxf[d]; limites_celula( d, ex, qual, xmin, xmax, xminf, xmaxf); No *pcf = p->fil[qual]; avancaT33(d, prof + 1, pcf, xminf, xmaxf, Dt, f); } assert((p->fil[0] == NULL) == (p->fil[1] == NULL)); if(p->fil[0] == NULL) { // segundo passo de Runge-Kutta 2 double tn2 = rk2nd3o( d, p, Dt, xmin, xmax); p->fv[1] = tn2; p->fval = tn2; p->fv[2] = tn2; p->fv[0] = tn2; } } */