I’m try to parallelize the stabilization term of pressure in Boussinesq problem, because when we start the problem in 3D we need the parallelization in order to reduce the time. The problem occurs when we use 2 or more processor, because the mean of the pressure is not zero, so we think to create the matrix associated to stabilization term in global and then restrict this matrix in local and introduce to our system, but we do not know how to do it.
Please share your code and what you want to do.
I want to parallelize this code and this is what I try and do not work correctly when I prove with more than 1 processor. When I put for example VhP1L it refers to fspace VhP1L(ThL,P1) in the local mesh.
matrix PIg = interpolate(VhP1L,VhP1dcL);
matrix IPg = interpolate(VhP1dcL,VhP1L);
matrix IPPIg = IPg*PIg;
matrix IPhP1dcP1L = IdP1dcL + (-1.)*IPPIg;
int[int] cs2=[3];
matrix D4X4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=1);
matrix D4Y4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=2);
matrix LPSpres;
VhP0 mk,tau,h,tKcod,tKTcod;
varf med(unused,v)=intN(Th)(1.*v);
real[int] medk=med(0,VhP0);
mk=sqrt(medk);
//COEF. ESTAB. CODINA
varf tauK(unused,v)=intN(Th)((u1tmp^2+u2tmp^2)*v);
varf tauKT(unused,v)=intN(Th)(((u1dcX)^2+(u1dcY)^2+(u2dcX)^2+(u2dcY)^2)*v);
real[int] tK=tauK(0,VhP0);
real[int] tKT=tauKT(0,VhP0);
tKcod = sqrt(tK);
tKTcod = sqrt(tKT);
tau = ((cc1*(Pr+((CShTriangle)^2)(tKTcod/mk))/hTriangle^2) + (cc2*(tKcod/mk)/hTriangle))^(-1.);
/////////////////////////////
//MATRIZ ESTAB. PRES.
/////////////////////////////
varf termPres(pp,q) = intN(Th)(tauppq);
matrix TermP = termPres(VhP1dcL,VhP1dcL);
matrix DDxx;
{
matrix Maux;
Maux = IPhP1dcP1L * D4X4P2L;
DDxx = Maux’ *TermP;
DDxx = DDxx *Maux;
}
matrix DDyy;
{
matrix Maux;
Maux = IPhP1dcP1L * D4Y4P2L;
DDyy = Maux’ * TermP;
DDyy = DDyy * Maux;
}
LPSpres = (DDxx + DDyy);
This I want to parallelize, but if I can’t I think that maybe I can calculate the LPSpres in global and then restrict to local introduce into my final matrix of the system.
What issue are you facing trying to parallelise the LPSpres matrix?
The issue is that when I parallelize the LPSpres matrix, in the form that I have included before, when I use more than one processor the pressure solution of my problem does not have mean zero, but when I only use one processor I have mean zero and I suppose that the issue probably can be in the parallelization.
Please share a runnable code.
I have divided the code in some .idp files so I have included the folder and you have to run the file EF.edp with ff-mpirun -np _ EF.edp in this folder.
2D_paralelo.zip (114.6 KB)
I have included all in one file in order to do easily.
verbosity =0.;
// MACROS 2D
macro Grad(u) [dx(u), dy(u)] //EOM
macro div(u1,u2) (dx(u1)+dy(u2)) //EOM
macro gradugradv(u1,u2,v1,v2) (dx(u1)*dx(v1)+dy(u1)*dy(v1)+dx(u2)*dx(v2)+dy(u2)*dy(v2))//EOM
macro ugradv(u1,u2,v1,v2) [[u1,u2]'*[dx(v1),dy(v1)],[u1,u2]'*[dx(v2),dy(v2)]] //EOM
macro ugradvw(u1,u2,v1,v2,w1,w2) (ugradv(u1,u2,v1,v2)'*[w1,w2]) //EOM
macro ugradT(u1,u2,temp) ([u1,u2]'*Grad(temp)) //EOM
macro antisimetric(u1,u2,v1,v2,w1,w2) (0.5*(ugraduw(u1,u2,v1,v2,w1,w2)-ugraduw(u1,u2,w1,w2,v1,v2))) //EOM
macro mod(u1x,u1y,u2x,u2y) (sqrt(u1x^2+u1y^2+u2x^2+u2y^2)) //EOM
// Si usamos vectores velocidad U,V,W en 2D con U=[U, UY], V=[V, VY], W=[W, WY]
// Las derivadas de U son UdX, UdY, UYdX y UYdY, si se quiere usar en versión vectorial (notacion mas especifica)
macro Div(U) (dx(U) + dy(U#Y)) //EOM
macro GradUGradV(U,V) (dx(U)*dx(V) +dy(U)*dy(V) + dx(U#Y)*dx(V#Y) + dy(U#Y)*dy(V#Y))//EOM
macro UGradV(U,V) [[U,U#Y]'*[dx(V),dy(V)], [U,U#Y]'*[dx(V#Y),dy(V#Y)]] //EOM
macro UGradVW(U,V,W) (UGradV(U,V)'*[W, W#Y]) //EOM
macro UGradT(U,temp) ([U,U#Y]'*Grad(temp)) //EOM
macro Antisimetric(U,V,W) (0.5*(UGradVW(U,V,W) - UGradVW(U,W,V))) //EOM
macro MOD(U) (sqrt(U#dX^2+U#dY^2+U#YdX^2+U#YdY^2)) //EOM
//// MACROS DEL PARALELO
macro pause() mpiBarrier(mpiCommWorld)// EOM
macro reduceSolution(uL,u,D,map)
{
real[int] aux(u[].n);aux=0;
uL[].*=D;
aux(map)=uL[];
u[]=0;
mpiAllReduce(aux,u[],mpiCommWorld,mpiSUM);
}
//EOM
NewMacro solvesystem(A,b,u)
{
Vh4P2L [u#L, u#LY, u#LT, u#LP];
MatAVh4P2 = A;
set(MatAVh4P2,sparams=sparamsv);
u#L[]=MatAVh4P2^-1*b;
reduceSolution(u#L,u,MatAVh4P2.D,mapVh4P2);
}
EndMacro
//MESH
//Malla cavidad -> Estructurada
int suelo = 1;
int paredd = 2;
int techo = 3;
int paredi = 4;
int[int] labs = [suelo,paredd,techo,paredi];
mesh Th=square(32,32,flags=3,label = labs);
mesh ThL;
func fx = 0.5*(1 + (tanh(2*(2*x-1)))/(tanh(2)));
func fy = 0.5*(1 + (tanh(2*(2*y-1)))/(tanh(2)));
Th = movemesh(Th,[fx,fy]);
ThL = Th;
//FESPACE
fespace Vh4P2(Th, [P2, P2, P2, P2]); //Cambiar P1 en la presión por P2 cuando se añada el termino de stab de presion
fespace Vh2P2(Th, [P2, P2]);
fespace Vh2P1(Th, [P1,P1]);
fespace Vh4P1(Th, [P1, P1, P1, P1]); //Cambiar P0 por P1 cuando se añada termino de estab presion
fespace Vh4P1dc(Th, [P1dc, P1dc, P1dc, P1dc]);
fespace VhP0(Th, P0);
fespace VhP1dc(Th, P1dc);
fespace VhP1(Th, P1);
fespace VhP2(Th, P2);
// fespace Vh2P2(Th, [P2,P2]);
fespace Vh4P2L(ThL, [P2, P2, P2, P2]); //Cambiar P1 en la presión por P2 cuando se añada el termino de stab de presion
fespace Vh2P2L(ThL, [P2, P2]);
fespace Vh2P1L(ThL, [P1,P1]);
fespace Vh4P1L(ThL, [P1, P1, P1, P1]); //Cambiar P0 por P1 cuando se añada termino de estab presion
fespace Vh4P1dcL(ThL, [P1dc, P1dc, P1dc, P1dc]);
fespace VhP0L(ThL, P0);
fespace VhP1dcL(ThL, P1dc);
fespace VhP1L(ThL, P1);
fespace VhP2L(ThL, P2);
//PARAMETERS
// Parametros del PETSC
int petsc = 1; //Para usar el paralelo = 1
macro dimension()2 //EOM
include "macro_ddm.idp";
load "PETSc";
Mat MatAVh4P2;
Mat MatAVhP1dc;
int[int] mapVh4P2,mapVhP1dc;
string sparamsv = "-pc_type lu -pc_factor_mat_solver_type mumps";
if(petsc){
int[int] myN2o;
macro ThLN2O() myN2o // EOM
buildDmesh(ThL);
mapVh4P2 = restrict(Vh4P2L,Vh4P2,myN2o);
{
macro def(i) [i, iY, iT, iP] //
macro init(i) [i, i, i, i] // EOM
createMat(ThL, MatAVh4P2, [P2, P2, P2, P2]);
}
mapVhP1dc = restrict(VhP1dcL,VhP1dc,myN2o);
{
macro def(i) [i] //
macro init(i) [i] //
createMat(ThL,MatAVhP1dc,P1dc);
}
}
real[int] auxD = MatAVhP1dc.D;
VhP0 mk,tau,h,tKcod,tKTcod;
varf med(unused,v)=intN(Th)(1.*v);
real[int] medk=med(0,VhP0);
mk[]=sqrt(medk);
h = hTriangle;
//Constantes varias
real CS=0.1;
real cc1=16.;
real cc2=sqrt(cc1);
real Pr=0.71; //Prandlt aire
real dt=1.e-2;//2.e-3
real dtt = 1./dt;
real epspen=1.e-6;//1.e-8
real Ramin=1.e3;//1.e5
real Ramax=1.e5;//1.e6
int ni=100;
int nIterations=2000;
real epsError = 1.e-6;//1.e-8
//INTERPOLATION MATRIX
//For global spaces
// Matrices de filtrado.
matrix IPhP2P1, IPh4P24P1, IPhP1dcP1;
{
matrix IdP2,Id4P2,IdP1dc,IdP1; // IdEFX: Matriz Identidad con dim(EFX) grados de libertad
{
VhP2 faux2=1.;
VhP1dc faux1=1.;
VhP1 fauxp1 = 1.;
IdP2 = faux2[];
IdP1dc = faux1[];
IdP1 = fauxp1[];
Id4P2=[[IdP2, 0, 0,0],[0, IdP2, 0, 0],[0, 0, IdP2,0],[0, 0, 0, IdP2]];
}
matrix PIg = interpolate(VhP1,VhP1dc); //(Id-πh) P1dc->P1->P1dc
matrix IPg = interpolate(VhP1dc,VhP1);
matrix IPPIg = IPg*PIg;
IPhP1dcP1 = IdP1dc + (-1.)*IPPIg;
matrix PI=interpolate(VhP1,VhP2);
matrix IP=interpolate(VhP2,VhP1);
matrix IPPI=IP*PI;
IPhP2P1 = IdP2 + (-1.)*IPPI; // (Id-πh) P2->P1->P2
matrix PI4=interpolate(Vh4P1,Vh4P2);
matrix IP4=interpolate(Vh4P2,Vh4P1);
matrix IPPI4=IP4*PI4;
IPh4P24P1 = Id4P2 + (-1)*IPPI4; // (Id-πh) 4P2->4P1->4P2
}
// Matrices de derivadas
matrix DX4P2, DY4P2; // ∂x, ∂y 4P2 -> 4P1dc en las 4 componentes
matrix DXP2, DYP2; // ∂x, ∂y P2 -> P1dc (El elemento finito solo tiene una componente)
matrix D1X4P2, D2X4P2, D3X4P2, D4X4P2; // ∂x 3P2 -> P1dc en la componente 1, 2 y 3 de 3P2
matrix D1Y4P2, D2Y4P2, D3Y4P2, D4Y4P2; // ∂y 3P2 -> P1dc en la componente 1, 2 y 3 de 3P2
{
matrix D1,D2,D3,D4; // Cada una de las componentes de EFX que derivo
int[int] c0=[0,-1,-1,-1];
int[int] c1=[-1,1,-1,-1];
int[int] c2=[-1,-1,2,-1];
int[int] c3=[-1,-1,-1,3];
D1 = interpolate(Vh4P1dc,Vh4P2,U2Vc=c0,op=1);
D2 = interpolate(Vh4P1dc,Vh4P2,U2Vc=c1,op=1);
D3 = interpolate(Vh4P1dc,Vh4P2,U2Vc=c2,op=1);
D4 = interpolate(Vh4P1dc,Vh4P2,U2Vc=c3,op=1);
DX4P2 = D1 + D2 + D3 + D4;
D1 = interpolate(Vh4P1dc,Vh4P2,U2Vc=c0,op=2);
D2 = interpolate(Vh4P1dc,Vh4P2,U2Vc=c1,op=2);
D3 = interpolate(Vh4P1dc,Vh4P2,U2Vc=c2,op=2);
D4 = interpolate(Vh4P1dc,Vh4P2,U2Vc=c3,op=2);
DY4P2 = D1 + D2 + D3 + D4;
/////////////////////////////////////////////////////////
//PARA LA ESTABILIDAD DE LA PRESION CAMBIAR CUANDO SE LA AÑADA
int[int] cs2=[0];
D1X4P2 = interpolate(VhP1dc,Vh4P2,U2Vc=cs2,op=1);
D1Y4P2 = interpolate(VhP1dc,Vh4P2,U2Vc=cs2,op=2);
cs2=[1];
D2X4P2 = interpolate(VhP1dc,Vh4P2,U2Vc=cs2,op=1);
D2Y4P2 = interpolate(VhP1dc,Vh4P2,U2Vc=cs2,op=2);
cs2=[2];
D3X4P2 = interpolate(VhP1dc,Vh4P2,U2Vc=cs2,op=1);
D3Y4P2 = interpolate(VhP1dc,Vh4P2,U2Vc=cs2,op=2);
cs2=[3];
D4X4P2 = interpolate(VhP1dc,Vh4P2,U2Vc=cs2,op=1);
D4Y4P2 = interpolate(VhP1dc,Vh4P2,U2Vc=cs2,op=2);
////////////////////////////////////////////////////////
DXP2 = interpolate(VhP1dc,VhP2,op=1);
DYP2 = interpolate(VhP1dc,VhP2,op=2);
}
// Matrices de derivada, con filtrado
matrix IPhD4X4P2 = IPhP1dcP1 * D4X4P2; // (I-πh) de ∂x componente 3 de 3P2 (3P2 (derivo componente 3) -> P1dc (filtro)-> P1 -> P1dc)
matrix IPhD4Y4P2 = IPhP1dcP1 * D4Y4P2; // (I-πh) de ∂y componente 3 de 3P2 (3P2 (derivo componente 3) -> P1dc (filtro)-> P1 -> P1dc)
//Matrices de derivada, con filtrada para una componente solo
matrix IPhDXP2 = IPhP1dcP1*DXP2;
matrix IPhDYP2 = IPhP1dcP1*DYP2;
//Matrices de filtrado con derivada
matrix DXP2IPh = DXP2*IPhP2P1;
matrix DYP2IPh = DYP2*IPhP2P1;
//For local spaces
// Matrices de filtrado.
matrix IPhP2P1L, IPh4P24P1L, IPhP1dcP1L;
{
matrix IdP2L,Id4P2L,IdP1dcL,IdP1L; // IdEFX: Matriz Identidad con dim(EFX) grados de libertad
{
VhP2L faux2=1.;
VhP1dcL faux1=1.;
IdP2L = faux2[];
IdP1dcL = faux1[];
VhP1L fauxp1 = 1.;
IdP1L = fauxp1[];
Id4P2L=[[IdP2L, 0, 0, 0],[0, IdP2L, 0, 0],[0, 0, IdP2L, 0],[0, 0, 0, IdP2L]];
}
matrix PIg = interpolate(VhP1L,VhP1dcL); //(Id-πh) P1dc->P1->P1dc
matrix IPg = interpolate(VhP1dcL,VhP1L);
matrix IPPIg = IPg*PIg;
IPhP1dcP1L = IdP1dcL + (-1.)*IPPIg;
matrix PI=interpolate(VhP1L,VhP2L);
matrix IP=interpolate(VhP2L,VhP1L);
matrix IPPI=IP*PI;
IPhP2P1L = IdP2L + (-1.)*IPPI; // (Id-πh) P2->P1->P2
matrix PI4=interpolate(Vh4P1L,Vh4P2L);
matrix IP4=interpolate(Vh4P2L,Vh4P1L);
matrix IPPI4=IP4*PI4;
IPh4P24P1L = Id4P2L + (-1)*IPPI4; // (Id-πh) 3P2->3P1->3P2
}
// Matrices de derivadas
matrix DX4P2L, DY4P2L, DZ4P2L; // ∂x, ∂y, ∂z 4P2 -> 4P1dc en las 4 componentes
matrix DXP2L, DYP2L, DZP2L; // ∂x, ∂y, ∂z P2 -> P1dc (El elemento finito solo tiene una componente)
matrix D1X4P2L, D2X4P2L, D3X4P2L, D4X4P2L; // ∂x 4P2 -> P1dc en la componente 1, 2, 3 y 4 de 4P2
matrix D1Y4P2L, D2Y4P2L, D3Y4P2L, D4Y4P2L; // ∂y 4P2 -> P1dc en la componente 1, 2, 3 y 4 de 4P2
matrix D1Z4P2L, D2Z4P2L, D3Z4P2L, D4Z4P2L; // ∂z 4P2 -> P1dc en la componente 1, 2, 3 y 4 de 4P2
{
matrix D1,D2,D3,D4; // Cada una de las componentes de EFX que derivo
int[int] c0=[0,-1,-1,-1];
int[int] c1=[-1,1,-1,-1];
int[int] c2=[-1,-1,2,-1];
int[int] c3=[-1,-1,-1,3];
D1 = interpolate(Vh4P1dcL,Vh4P2L,U2Vc=c0,op=1);
D2 = interpolate(Vh4P1dcL,Vh4P2L,U2Vc=c1,op=1);
D3 = interpolate(Vh4P1dcL,Vh4P2L,U2Vc=c2,op=1);
D4 = interpolate(Vh4P1dcL,Vh4P2L,U2Vc=c3,op=1);
DX4P2L = D1 + D2 + D3 + D4;
D1 = interpolate(Vh4P1dcL,Vh4P2L,U2Vc=c0,op=2);
D2 = interpolate(Vh4P1dcL,Vh4P2L,U2Vc=c1,op=2);
D3 = interpolate(Vh4P1dcL,Vh4P2L,U2Vc=c2,op=2);
D4 = interpolate(Vh4P1dcL,Vh4P2L,U2Vc=c3,op=2);
DY4P2L = D1 + D2 + D3 + D4;
D1 = interpolate(Vh4P1dcL,Vh4P2L,U2Vc=c0,op=3);
D2 = interpolate(Vh4P1dcL,Vh4P2L,U2Vc=c1,op=3);
D3 = interpolate(Vh4P1dcL,Vh4P2L,U2Vc=c2,op=3);
D4 = interpolate(Vh4P1dcL,Vh4P2L,U2Vc=c3,op=3);
DZ4P2L = D1 + D2 + D3 + D4;
/////////////////////////////////////////////////////////
int[int] cs2=[0];
D1X4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=1);
D1Y4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=2);
D1Z4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=3);
cs2=[1];
D2X4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=1);
D2Y4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=2);
D2Z4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=3);
cs2=[2];
D3X4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=1);
D3Y4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=2);
D3Z4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=3);
cs2=[3];
D4X4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=1);
D4Y4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=2);
D4Z4P2L = interpolate(VhP1dcL,Vh4P2L,U2Vc=cs2,op=3);
////////////////////////////////////////////////////////
DXP2L = interpolate(VhP1dcL,VhP2L,op=1);
DYP2L = interpolate(VhP1dcL,VhP2L,op=2);
DZP2L = interpolate(VhP1dcL,VhP2L,op=3);
}
// Matrices de derivada, con filtrado
matrix IPhD4X4P2L = IPhP1dcP1L * D4X4P2L; // (I-πh) de ∂x componente 4 de 4P2 (4P2 (derivo componente 4) -> P1dc (filtro)-> P1 -> P1dc)
matrix IPhD4Y4P2L = IPhP1dcP1L * D4Y4P2L; // (I-πh) de ∂y componente 4 de 4P2 (4P2 (derivo componente 4) -> P1dc (filtro)-> P1 -> P1dc)
matrix IPhD4Z4P2L = IPhP1dcP1L * D4Z4P2L; // (I-πh) de ∂z componente 4 de 4P2 (4P2 (derivo componente 4) -> P1dc (filtro)-> P1 -> P1dc)
// Global a local
//Data
//Definimos el intervalo de Rayleigh que vamos a calcular
real[int] vRa(ni+1);
for(int i=0;i<=ni;i++){
vRa(i)=Ramin+i*(Ramax-Ramin)/(ni);
//vRa(i)=10^(3+i*(1.)/50);
}
//BOussinesq problem
//BUCLE RESOLVER EL PROBLEMA
for (int i=0;i<=ni;i++)
{
real Ra=vRa(i);
if(mpirank==0){cout<<"i: "<<i<<"-----"<<Ra<<endl;}
//lifting for temperature
VhP2 G = 1-x; //Levantamiento para la temperatura
Vh4P1dc [Gfv1,Gfv2,GfT,Gfp];
VhP2 Gfaux;
VhP1 Gpaux;
Vh4P2 [GG1,GG2,GGT,GGP]=[Gfaux,Gfaux,G,Gfaux]; //Levantamiento temperatura
//Fixed matrix
//// Construimos la matrices fijas para EF
varf vNS([uu1,uu2,tt,pp],[v1,v2,zt,q])=
intN(ThL)((1./dt)*(uu1*v1+uu2*v2+tt*zt)
-div(v1,v2)*pp+div(uu1,uu2)*q
+Pr*gradugradv(uu1,uu2,v1,v2)
+Grad(tt)'*Grad(zt)
+1.e-8*pp*q
-Pr*Ra*tt*v2)
//Segundo miembro
+intN(ThL)(Pr*Ra*G*v2
-Grad(G)'*Grad(zt))
+on(suelo,paredd,paredi,techo, uu1=0., uu2=0.)
+on(paredi, tt=0.)+on(paredd, tt=0.);
//
matrix AvNS;
AvNS=vNS(Vh4P2L,Vh4P2L);
real[int] bNSf(Vh4P2L.ndof);
bNSf=vNS(0,Vh4P2L);
//Definimos la solucion y la solucion en la etapa anterior
Vh4P2 [u,uY,uT,uP]; //Solucion
//[u1,u2,t,p] = [UX,UY,t,p];
u[]=0.;
// Preliminares a la iteracion
real itercontrol1 = 2.*epsError + 1.;
real iterrelative1 = itercontrol1;
VhP2 u1tmp,u2tmp; //Solucion en la etapa anterior
{
u1tmp=u;
u2tmp=uY;
}
//Derivadas en la etapa anterior para el VMS-SMago
VhP2 udcAux;
VhP1dc u1dcX,u1dcY,u2dcX,u2dcY,Gr;
udcAux[] = IPhP2P1*u1tmp[];
u1dcX[] = DXP2*udcAux[];
u1dcY[] = DYP2*udcAux[];
udcAux[] = IPhP2P1*u2tmp[];
u2dcX[] = DXP2*udcAux[];
u2dcY[] = DYP2*udcAux[];
Gr=mod(u1dcX,u1dcY,u2dcX,u2dcY);
//Bucle para llegar a la estabilizacion
for(int ii=0;ii<=nIterations && (iterrelative1 > epsError); ii++)
{
///////////////////////////////////
//Matriz término de Conveccion N-S
///////////////////////////////////
varf vNSb([uu1,uu2,tt,pp],[v1,v2,zt,q])=
intN(ThL)(ugradv(u,uY,uu1,uu2)'*[v1,v2]
+([u,uY]'*Grad(tt))*zt)
//Seg miembro
+intN(ThL)((1./dt)*(u*v1+uY*v2+uT*zt)
-([u,uY]'*Grad(G))*zt)
;
matrix AvNSb=vNSb(Vh4P2L,Vh4P2L);
real[int] bNSb(Vh4P2L.ndof);
bNSb = vNSb(0,Vh4P2L);
////////////////////////////////////
////MATRIZ SMAGORINSKY VEL FILTRADA
////////////////////////////////////
varf VMSSma([u1l,u2l,unusedTT,unusedpp],[v1l,v2l,zl,ql]) =
intN(ThL)(((CS*hTriangle)^2)*Gr*(u1l*v1l+u2l*v2l));
matrix M = VMSSma(Vh4P1dcL,Vh4P1dcL);
matrix Sma;
matrix DXX;
matrix DYY;
{
matrix Maux;
Maux = DX4P2L*IPh4P24P1L;
DXX = (Maux')*M;
DXX = DXX*Maux;
Maux = DY4P2L*IPh4P24P1L;
DYY = (Maux')*M;
DYY = DYY*Maux;
}
Sma = DXX + DYY;
////////////////////////////////////
////MATRIZ SMAGORINSKY TEMP FILTRADA
////////////////////////////////////
varf VMSSmaT([un1l,un2l,TT,unusedpp],[v1l,v2l,zl,ql]) =
intN(ThL)((1./Pr)*((CS*hTriangle)^2)*Gr*(TT*zl))
-intN(ThL)((1./Pr)*((CS*hTriangle)^2)*Gr*(GfT*zl));
//FALTA SEGUNDO MIEMBRO DEL LEVANTAMIENTO DE LA TEMPERATURA
matrix MT = VMSSmaT(Vh4P1dcL,Vh4P1dcL);
matrix SmaT;
matrix DXXT;
matrix DYYT;
real[int] bTx(Vh4P2L.ndof),bTy(Vh4P2L.ndof);
{
matrix Maux,MauxG;
Maux = DX4P2L*IPh4P24P1L;
MauxG = DX4P2*IPh4P24P1;
DXXT = (Maux')*MT;
DXXT = DXXT*Maux;
Gfv1[]=MauxG*GG1[];
real[int] bTs=VMSSmaT(0,Vh4P1dcL);
bTx=(Maux')*bTs;
Maux = DY4P2L*IPh4P24P1L;
MauxG = DY4P2*IPh4P24P1;
DYYT = (Maux')*MT;
DYYT = DYYT*Maux;
Gfv1[]=MauxG*GG1[];
bTs=VMSSmaT(0,Vh4P1dcL);
bTy=(Maux')*bTs;
}
SmaT = DXXT + DYYT;
////////////////////////////////////
/////MATRIZ LPS
////////////////////////////////////
matrix LPSpres;
//COEF. ESTAB. CODINA
varf tauK(unused,v)=intN(Th)((u1tmp^2+u2tmp^2)*v);
varf tauKT(unused,v)=intN(Th)(((u1dcX)^2+(u1dcY)^2+(u2dcX)^2+(u2dcY)^2)*v);
real[int] tK=tauK(0,VhP0);
real[int] tKT=tauKT(0,VhP0);
tKcod[] = sqrt(tK);
tKTcod[] = sqrt(tKT);
tau = ((cc1*(Pr+((CS*hTriangle)^2)*(tKTcod/mk))/hTriangle^2) + (cc2*(tKcod/mk)/hTriangle))^(-1.);
/////////////////////////////
//MATRIZ ESTAB. PRES.
/////////////////////////////
varf termPres(pp,q) = intN(Th)(tau*pp*q);
matrix TermP = termPres(VhP1dcL,VhP1dcL);
matrix DDxx;
{
matrix Maux;
Maux = IPhD4X4P2L;
DDxx = Maux' *TermP;
DDxx = DDxx *Maux;
}
matrix DDyy;
{
matrix Maux;
Maux = IPhD4Y4P2L;
DDyy = Maux' * TermP;
DDyy = DDyy * Maux;
}
LPSpres = (DDxx + DDyy);
///////////////////////////
//Matriz final
///////////////////////////
matrix Afinal;
Afinal=AvNSb + AvNS;
Afinal+= Sma;
Afinal+= SmaT;
Afinal+= LPSpres;
////////////////////////
//Segundo miembro
///////////////////////
real[int] bNS(Vh4P2L.ndof);
bNS = vNSb(0,Vh4P2L);
bNS += bNSf;
bNS += bTx;
bNS += bTy;
solvesystem(Afinal,bNS,u)
//Pintamos la solucion de cada iteracion
// VhP2 TEMP=t+G;
// plot(TEMP,wait=0,cmm="Solucion en la iteracion "+ii);
VhP2 eu1,eu2;
eu1=u-u1tmp;
eu2=uY-u2tmp;
//Inicio calculo normas
itercontrol1 = sqrt(int2d(Th)(gradugradv(eu1,eu2,eu1,eu2)));
iterrelative1 = itercontrol1;///itertmp1;
if(mpirank==0){cout << itercontrol1 << endl;
cout << "media pres:" << intN(Th)(uP) << endl;}
// cout << "------------------------------------" << endl;
// cout << "ERROR RELATIVO PARA u1: "<<iterrelative1<<endl;
//ACTUALIZ.
u1tmp=u;
u2tmp=uY;
udcAux[] = IPhP2P1*u1tmp[];
u1dcX[] = DXP2*udcAux[];
u1dcY[] = DYP2*udcAux[];
udcAux[] = IPhP2P1*u2tmp[];
u2dcX[] = DXP2*udcAux[];
u2dcY[] = DYP2*udcAux[];
Gr=mod(u1dcX,u1dcY,u2dcX,u2dcY);
}
VhP2 T = uT+G;
plot(T,cmm="Temperatura para Ra="+Ra,fill=1,value=1);
plot(u,cmm="Velocidad x para Ra="+Ra,fill=1,value=1);
plot(uY,cmm="Velocidad y para Ra="+Ra,fill=1,value=1);
//plot([u1,u2],cmm="Modulo velocidad para Ra = "+Ra,wait=1);
plot(uP,cmm="Presion para Ra = "+Ra,fill=1,value=1);
}
It’s extremely difficult to debug your code (too long). Small comment, when you do this:
matrix Afinal;
Afinal=AvNSb + AvNS;
Afinal+= Sma;
Afinal+= SmaT;
Afinal+= LPSpres;
Your boundary conditions may not be enforced properly. Have you double checked that there is no error there?
I know that the code is so long sorry. I have checked this and the boundary condition work correctly because I only included in one “varf”.
Thank so much
Yes, but you use the default tgv value of 10^{30}, so that could be the source of one problem (I’m not saying this is the only problem in the script).
I change all the varf with matrix A = varf(Vh,Vh,tgv=-2) and the problem keep going.
That’s making things worse, you need one with -1, and all the others with -10.
I prove it and it does not work.