Domain and Pressure

I try to solve the base flow of a hated cylinder (2D)(natural convection) and try to understand the outer domain and pressure relationship. However, the pressure values sometimes increase and decrease when the domain increases.
I use the Newton method for iteration,
the code is;

fespace Mh(Th,[P2,P2,P1,P2]);
fespace Nh(Th,[P2,P2,P1]);
fespace XXh(Th,P2);
fespace Xh(Th,P1);
Mh [ux,uy,p,T]=[0,0,0,0];
Mh [dux,duy, dp, dT];
Mh [psix,psiy,psip, psiT];

macro grad1(v) [dx(v),dy(v)] //
macro div(ux,uy) (dx(ux)+dy(uy)) //
macro grad2(ux,uy) [dx(ux),dy(uy)] //
//mathematical solution with weak form (psi functions)

varf dNS([dux,duy,dp,dT],[psix,psiy,psip,psiT])=
dx(dT)) // heat equation (cylindrical and weak form)

 +int2d(Th)( x^2*psip*div(dux,duy) + x*dux*psip)    // continuity (cylindrical and weak form)

 +int2d(Th)( x^2*psiy*(uy*dy(duy)+duy*dy(uy)
	        + x^2*psix*(uy*dy(dux)+duy*dy(ux)
	       	          +ux*dx(dux)+dux*dx(ux)))       //u.grad(u)(convective term of N-S)

 +int2d(Th)( x^2*Pr*(2*grad2(dux,duy)' *grad2(psix,psiy)+div(duy,dux)*div(psiy,psix))
	   + x*Pr*(2*psix*dx(dux)+psiy*div(duy,dux))
	   + 2*Pr*dux*psix
	   - dp*(x^2*div(psix,psiy) + 2*x*psix)
	   - x^2*1e-8*dp*psip )                          // diffusion and pressure term of N-S

 -int2d(Th)(x^2*Ra*Pr*psiy*dT)                      // bousinessq approximation for buoyancy

//boundary conditions

    + on(axis,dux=0-ux, duy=0-uy, dT=0-T) 
+ on(outer,dux=0-ux, duy=0-uy, dT=0-T)													
+ on(sphere1,  dux=0-ux, duy=0-uy, dT=0-T)
+ on(sphere2,  dux=0-ux, duy=0-uy, dT=1-T);

//Solution for main variables
varf NS([dux,duy,dp,dT],[psix,psiy,psip,psiT])=
dx(T)) // heat equation (cylindrical and weak form)

 + int2d(Th)( x^2*psip*div(ux,uy) + x*ux*psip)     // continuity (cylindrical and weak form)

 +int2d(Th)( x^2*psix*(uy*dy(ux)
	          + x^2*psiy*(uy*dy(uy)
	       	             +ux*dx(uy)))	           //u.grad(u)(convective term of N-S)
 +int2d(Th)( Pr*(x^2*(2*grad2(ux,uy)' *grad2(psix,psiy)+div(uy,ux)*div(psiy,psix))
	   + x*(2*psix*dx(ux)+psiy*div(uy,ux))
	   + 2*ux*psix)
       - p*(x^2*div(psix,psiy) + 2*x*psix )
       - x^2*1e-8*p*psip )                        // diffusion and pressure term of N-S

 -int2d(Th)(x^2*Ra*Pr*psiy*T);                  //bousinessq approximation for buoyancy

//iteration part
real iterNUM=1; //number of iterations
real errNS=1; //Navier-Stokes equ error
real tolST=1e-8; //target value for iteration
real normux0 = 1; //initial velocity
real normuy0 = 1; //initial velocity
real normp0=1; //initial pressure
real normT0=1; //initial temperature
real normNS=1; //initial N-S
real deltaerr=1; //delta of error
real normp=1; //target for pressure
real normT=1; //target for temperature
real normux; //target for velocity
real normuy; //target for velocity

while (normNS > tolST)
cout << " Iteration number =" << iterNUM <<endl;
cout << " Solving Navier-Stokes Equations" <<endl;
cout << “======================================================” <<endl;

matrix A = dNS(Mh, Mh,solver=sparsesolver);//matrix solution
real[int] b1 = NS(0,Mh); //constant part coming from initial guess
real[int] b2 = dNS(0,Mh); //constant part coming from Dirichlet BC's
real[int] b = b1 - b2;
real[int] q = (Mh.ndof);
q = A^-1*b ;
dux[] =q;

ux[] = ux[]- dux[];
errNS= normNS;
//plot([ux,uy], wait= false, value=true, fill=true);
// updating the L2 norms
normux = int2d(Th)(x*((dux)*(dux) ));  // L2 norm
normuy = int2d(Th)(x*( (duy)*(duy) ));  // L2 norm
normp = int2d(Th)(x*(dp)*(dp));                      // L2 norm
normT = int2d(Th)(x*(dT)*(dT));				// L2 norm
normNS = sqrt(normux+normp+normT+normuy); 
deltaerr= normNS-errNS;

cout << " NS residual = " << normNS <<endl;	
cout << " velocity error X= " << sqrt(normux) <<endl;
cout << " velocity error Y= " << sqrt(normuy) <<endl;
cout << " pressure error= " << sqrt(normp) <<endl;
cout << " temperature error= " << sqrt(normT) <<endl;

cout << " NS total error= " << errNS <<endl;
cout << "======================================================" <<endl;

cout << "======================================================" <<endl;
cout << " Iteration number= " << iterNUM <<endl;
cout << " Stokes residual= " << normNS <<endl;
cout << "======================================================" <<endl;

Can somebody help so that I can understand the reason?

I haven’t looked at everything, but part of the issue must be that your operators are wrong. Cylindrical coordinates changes things (eg weak form of laplacian)

I do not understand, Can you explain a bit?