How to solve Possion equation have a parameter h for different values of h using loop

M_Usman · September 20, 2023, 3:55am

Dear, I want to solve the Poisson equation involving h for different values of h given in Pram (in the following code). Then, I used a loop to pick each value from Pram as a value for h and then solve the problem. I want to compute AvgNu at different values of h. Kindly see my code and guide me to achieve the required results.

real Mesh=60;
real ra=0.1;
border C(t=0,1){x=0.5+racos(2pit);y=0.2sqrt(3)+rasin(2pit);}
mesh Th=buildmesh(C(Mesh));
plot(Th,wait=1);
real Pram=[1,2,3];
for(int i;i<=3;i++)
{
real h=Pram[i];
fespace Vh(Th,P2);
func f= xy;
Vh U,V;
problem Poisson(U,V,solver=LU)=
int2d(Th)(dx(U)*dx(V) + dy(U)*dy(V))

hint2d(Th)( fV )

on(C,U=1) ;
Poisson;
real AvgNu;
macro grad(u) [dx(u),dy(u)]//
AvgNu[i]=-int1d(Th,C)((grad(U))'*[N.x,N.y]);
cout<<"Avg = "<<AvgNu[i]<<endl;
}

frederichecht · September 20, 2023, 9:27am

They is no problem

use the following writing:

problem Poisson(U,V,solver=LU)=
int2d(Th)(dx(U)*dx(V) + dy(U)*dy(V))
- int2d(Th)( h* f * V )
on(C,U=1) ;

in Poisson problem

M_Usman · September 20, 2023, 12:51pm

Dear Professor @frederichecht thanks for your response. I think you don’t understand my point. I want to solve the Poisson equation for different values of h using the loop. For example,
for i=1, the code should solve the Poisson problem for h = 1 and compute the value for AvgNu,
for i=2, the code should solve the Poisson problem for h = 2 and compute the value for AvgNu, and
for i=3, the code should solve the Poisson problem for h = 3 and compute the value for AvgNu.

frederichecht · September 20, 2023, 1:26pm

No , for( i = 0; i<3; ++i) read the doc please.

M_Usman · September 20, 2023, 3:41pm

That means we can’t solve using the loop.

frederichecht · September 21, 2023, 8:23am

the corrected code with only factorisation (init=0) because the matrix is constante here.

// M_Usman Muhammad Usman
real Mesh=60;
real ra=0.1;
border C(t=0,1){x=0.5+ra*cos(2*pi*t);y=0.2*sqrt(3)+ra*sin(2*pi*t);}
mesh Th=buildmesh(C(Mesh));
plot(Th,wait=1);
real[int]  Pram=[1,2,3],AvgNu(3);
real h;
fespace Vh(Th,P2);
func f= x*y;
Vh U,V;
int i=0;
problem Poisson(U,V,solver=LU,init=i)=
    int2d(Th)(dx(U)*dx(V) + dy(U)*dy(V))
    - int2d(Th)( h*f*V )
    + on(C,U=1) ;
    verbosity=0; 
for( i=0;i<3;i++)
{
  h=Pram[i];

  Poisson;
  plot(U,wait=1);
  macro grad(u) [dx(u),dy(u)]//
AvgNu[i]=-int1d(Th,C)((grad(U))'*[N.x,N.y]);
cout<<i << "Avg = "<<AvgNu[i]<<endl;
}

the output :

ff++-13/bin/FreeFem++ 'Usman.edp' -glut '/usr/local/ff++-13/bin/ffglut' ;exit;
-- FreeFem++ v4.13 (Mer 20 sep 2023 11:48:02 CEST - git v4.13-141-g86b66a0e)
   file : Usman.edp
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : // M_Usman Muhammad Usman
    2 : real Mesh=60;
    3 : real ra=0.1;
    4 : border C(t=0,1){x=0.5+ra*cos(2*pi*t);y=0.2*sqrt(3)+ra*sin(2*pi*t);}
    5 : mesh Th=buildmesh(C(Mesh));
    6 : plot(Th,wait=1);
    7 : real[int]  Pram=[1,2,3],AvgNu(3);
    8 : real h;
    9 : fespace Vh(Th,P2);
   10 : func f= x*y;
   11 : Vh U,V;
   12 : int i=0;
   13 : problem Poisson(U,V,solver=LU,init=i)=
   14 :     int2d(Th)(dx(U)*dx(V) + dy(U)*dy(V))
   15 :     - int2d(Th)( h*f*V )
   16 :     + on(C,U=1) ;
   17 :     verbosity=0; 
   18 : for( i=0;i<3;i++)
   19 : {
   20 :   h=Pram[i];
   21 : 
   22 :   Poisson;
   23 :   plot(U,wait=1);
   24 :   macro grad(u) [dx(u),dy(u)] )  //
   25 : AvgNu[i]=-int1d(Th,C)((grad(U)  [dx(U),dy(U)])'*[N.x,N.y]);
   26 : cout<<i << "Avg = "<<AvgNu[i]<<endl;
   27 : } sizestack + 1024 =2168  ( 1144 )

  --  mesh:  Nb of Triangles =    638, Nb of Vertices 350
0Avg = 0.00545748
1Avg = 0.010915
2Avg = 0.0163724

M_Usman · September 21, 2023, 12:31pm

Thank you, dear Professor @frederichecht . For nonlinear PDEs, I got the same value for AvgNu. Maybe I did something wrong in the following code. Kindly have a look and guide me to overcome this error.

real Mesh=60;
real ra=0.1;
border C(t=0,1){x=0.5+racos(2pit);y=0.2sqrt(3)+rasin(2pit);}
mesh Th=buildmesh(C(Mesh));
plot(Th,wait=1);
real a=1;
real err=1;
real[int] Pram=[1,2,3],AvgNu(3);
real b;
fespace Vh(Th,P2);
func f= xy;
Vh U,U0,V;
U0=0;
int i=0;
problem Poisson(U,V,solver=LU,init=i)=
int2d(Th)(dx(U)dx(V) + dy(U)dy(V))
-int2d(Th)(abfU0UV)
+on(C,U=1) ;
verbosity=0;
for( i=0;i<3;i++)
{
b=Pram[i];
for(int ij=0; (ij)<=10; ij++)
while (err>1e-4)
{
Poisson;
err=int2d(Th)(U-U0)^2;
U0 = U;
endl;
}
macro grad(u) [dx(u),dy(u)]//
AvgNu[i]=-int1d(Th,C)((grad(U))'[N.x,N.y]);
cout<<i << "Avg = "<<AvgNu[i]<<endl;
}

frederichecht · September 21, 2023, 2:40pm

For non linear equation first you have to choose which algorithm you use to solve your problem
fixed point , Newton, …

so you put a loop in a loop…

M_Usman · September 21, 2023, 3:06pm

I solved the nonlinear problem as:
for(int ij=0; (ij)<=10; ij++)
while (err>1e-4)
{
Poisson;
err=int2d(Th)(U-U0)^2;
U0 = U;
endl;
}
and then put it in the loop to compute the solution for various values of b. The got the solution successfully but the solution is same for each value of b.

frederichecht · September 21, 2023, 3:20pm

Normal, your problem do not change in the iteration process, and the err must set before big value.

M_Usman · September 21, 2023, 3:49pm

I already set err=1 here, then I used it in the iterative process.

M_Usman · September 21, 2023, 6:01pm

Dear Professor @frederichecht , I also defined err=1 just after the value b and before starting the second loop. But, I got the same values for AvgNu. Where I have to define the value err so that I get the correct results for AvgNu. I am waiting for your positive response.

M_Usman · September 22, 2023, 12:16pm

Dear Professor @frederichecht , what do you meant by big value.

frederichecht · September 22, 2023, 12:40pm

bigger than 1e-4 to enter , in the loop. so 1 for exemple…

M_Usman · September 22, 2023, 1:10pm

Dear Professor @frederichecht , I followed your instructions (code given below) but found the same results for AvgNu.

real Mesh=100;
real ra=0.1;
border C(t=0,1){x=0.5+racos(2pit);y=0.2sqrt(3)+rasin(2pit);}
mesh Th=buildmesh(C(Mesh));
plot(Th,wait=1);
real a=1;
real[int] Pram=[1,20,30],AvgNu(3),E=[1,1,1];
real b;
fespace Vh(Th,P2);
func f= xy;
Vh U,U0,V;
U0=0;
int i=0;
problem Poisson(U,V,solver=LU,init=i)=
int2d(Th)(dx(U)dx(V) + dy(U)dy(V))
-int2d(Th)(abfU0UV)
+on(C,U=1) ;
verbosity=0;
for( i=0;i<3;i++)
{
b=Pram[i];
real err=1;
for(int ij=0; (ij)<=10; ij++)
while (err>1e-7)
{
cout<<err<<endl;
Poisson;
err=int2d(Th)(U-U0)^2;
U0 = U;
endl;
}
macro grad(u) [dx(u),dy(u)]//
AvgNu[i]=-int1d(Th,C)((grad(U))'[N.x,N.y]);
cout<<i << "Avg = "<<AvgNu[i]<<endl;
}

M_Usman · September 22, 2023, 2:27pm

Dear Professor @frederichecht , kindly see this and guide me to fix the query.

frederichecht · September 28, 2023, 9:09am

I have no idea of your No linear problem, because you only put linear problem.
So it is difficult to help you.

Topic		Replies	Views
Problem vs varf General Discussion	3	46	September 30, 2024
Let the exact solution is unknown for the Poisson equation. How to compute the L2 error estimate u_h - u_h/2 General Discussion	4	316	August 14, 2023
Type mismatch when using varf General Discussion	0	274	August 20, 2021
Solution to the Poisson Boltzmann equation General Discussion	3	122	June 28, 2024
Lagrange Multipliers and Adaptive Meshing General Discussion	1	617	September 16, 2020

How to solve Possion equation have a parameter h for different values of h using loop

Related topics