mesh of a right paving stone and a cube

frederichecht · May 3, 2023, 9:18am

One idea is to introduce a auxiliary variable c = ||v||
when you can see the problem as minimisation under constraint and use IPOPT to solve the non
linear problem for fully implicite schema in time

but otherwise you can explicite c a do a semi implicite schema in time.

start you 1d or 2d cas to see if the scheme converge (no idea)
try do solve a problem with no dependance in time,

frederichecht · May 3, 2023, 10:02am

Juste a small example

aweni.edp (582 Bytes)

with a fixed point iteration,
remark, the you can do Newton iteration, in 10 mn.

aweni · May 3, 2023, 11:52am

Thank you very much Professor
I would like to know what this part of the code is for so that I can do it for the 2d
case before moving on to the time dependent case.

up-= u;
real err= up.linfty;
cout << iter << " c "<< c<< " err = "<< err<< endl;
if(err < 1e-6) break;

frederichecht · May 3, 2023, 12:38pm

// here up = previous u 
up-= u; // in up = up -u 
real err= up.linfty; //  nom of the diff up-u ..
cout << iter << " c "<< c<< " err = "<< err<< endl;
if(err < 1e-6) break;

aweni · May 3, 2023, 12:44pm

Sorry for the 3d i want to write.
I would like to understand better so that I don’t have to deal with this kind of problem anymore.

up-= u;
real err= up.linfty;
cout << iter << " c "<< c<< " err = "<< err<< endl;
if(err < 1e-6) break;

frederichecht · May 3, 2023, 12:58pm

the 3d version
aweni.edp (605 Bytes)

marchywka · May 3, 2023, 12:58pm

AFAIK, linifinty is just a fancy way of saying maximum absolute error
but you can presumably pick a criteria you like.

The static case is intersting ( or would be interesting )f the signs
change ) as the period is related to the amplitude. Wavelengths
tend to be fixed by geometry and boundary conditions too.

Since you speciffied initial conditions over omega you already have an
expression for the time derivative which is probably not entirely zero.

You could write the thing as v=v_0 + dt*( \laplacian v + cv)
with c derived from v_0 and hopefully in the limit of small dt results make sense
of as I have done do a 2d+1 case where the “z” dimension in the mesh is actually time.

aweni · May 3, 2023, 2:16pm

this is what I wrote and what I obtain justifies well what I obtained in the theory :

mesh Th=square(10,10);
// solve partial_t u -Delta u + c u = f , l = || u||_2
// u = g on gamma
// arg min int2d(Th) 1/2 grad(u) grad(u) + 1/2 int2d(Th) u^2
macro grad(u) [dx(u),dy(u)] ) //
func g = 0;
func f= 0.;
// solution u = 1
int kk=0;
real c = 0, dt= 0.1,Tfinal=10;
real[int] instT(floor(Tfinal/dt)+1);// vecteur qui contiendra les instants de temps tn
real[int] NORML2(floor(Tfinal/dt)+1);// vecteur qui contiendra les différents valeurs de E(tn)
fespace Vh(Th,P1);

Vh u=x^2+y^2+1 ,up,v ;
for (real t=0;t<=Tfinal ;t+=dt)// boucle du temps t
{

for (int iter= 0; iter < 100; ++iter)
{
up= u;
c = sqrt(int2d(Th)(u*u));
solve Pb(u,v) = int2d(Th) ((1/dt)uv)

int2d(Th) ( [dx(u),dy(u)]'* [dx(v),dy(v)] - cuv)

int2d(Th) ((1/dt)upv) - int2d(Th)(fv) + on(1,2,3,4, u=g);
up[]-= u[];
real err= up[].linfty;
cout << iter << " c "<< c<< " err = "<< err<< endl;
if(err < 1e-6) break;
plot(u, wait=1,dim=3);
}
//Calcul of E(t)= |u|_L^2(Omega)+|nabla u|_L^2(Omega)
instT[kk]= t;
NORML2[kk] = int2d(Th)( uu)

int2d(Th)([dx(u),dy(u)]'* [dx(u),dy(u)]);
cout <<“t” << instT(kk)<< endl;

cout <<"E(t) " << NORML2(kk) << endl;
//save of E(t) and t

if ( !(kk % 10)) // kk%10= le reste de la division kk/10
{
{
ofstream fout(“energie.txt”);
for(int i=0;i <= kk ; i++)
fout << instT(i)<<" " << NORML2(i)<<endl ;
}
}
kk+=1 ;

}

aweni · May 3, 2023, 2:24pm

You could write the thing as v=v_0 + dt*( \laplacian v + cv)
with c derived from v_0 and hopefully in the limit of small dt results make sense
of as I have done do a 2d+1 case where the “z” dimension in the mesh is actually time.

I don’t understand it is the L^2 norm of v .
I don’t see how it is the derivative of v0 and the derivative is with respect to which variable

marchywka · May 3, 2023, 4:59pm

The L2 norm is in “c” which I thought was used above. where you had a time
derivaive. Let me look at the code you just posted when I get a chance.

marchywka · May 3, 2023, 6:30pm

Can you post your code as an attached file? Copy/paste makes a lot ofthings to clean up.
Thanks.
What does it do in the limit of dt → zero?

aweni · May 3, 2023, 11:06pm

I will send you the file.
My objective was to see that the energy |v\|_L^2 decreases towards 0 when the time t tends towards infinity.

(Attachment ss is missing)

marchywka · May 3, 2023, 11:57pm

It says attachment is missing or something similar.

Typically you can test your discretization by letting the steps get smaller until
the results don’t change. Some of these things are unstable with steps too
big also.

aweni · May 4, 2023, 7:36am

ok sorry. I will resend it

ss.edp (1.57 KB)

marchywka · May 4, 2023, 10:26am

I had to edit a lot of your upload so I did not do a diff but there are a number of
issues. I used one of the idp files to move your parameters to the command line,
see the uploaded fie.
First, you appear to be using “up” for convergence testing and time stepping.
Then, dt is in the wrong place. and you can see the huge difference
by changing the time step although you expect a lot of this when
the time step is too big anywya I thought your initial conditions were given and you
appear to iterate those away. I’ll see if I can change it a bit if I get time later.

I always end up with simple algebra mistakes.
The integration by parts for the weak form is almost magic- it seems
to remove a lot of terms when more things vary spatially- but it is also
another source of sign mistakes in algebra

 FreeFem++ forum7.edp  | grep "t=" | grep " E(t) " | more
   43 :  cout <<"t=" << instT(kk)<<" E(t) " << NORML2(kk) << endl;
t=0 E(t) 0.034371
t=0.001 E(t) 0.00063618
t=0.002 E(t) 1.16154e-05
t=0.003 E(t) 2.11678e-07
t=0.004 E(t) 3.85664e-09
t=0.005 E(t) 5.8702e-10
t=0.006 E(t) 5.63973e-10
t=0.007 E(t) 5.4183e-10
t=0.008 E(t) 5.20557e-10
t=0.009 E(t) 5.00119e-10
marchywka@happy:/home/documents/cpp/proj/freefem/play$ FreeFem++ forum7.edp -dt .0005  | grep "t=" | grep " E(t) " | more
   43 :  cout <<"t=" << instT(kk)<<" E(t) " << NORML2(kk) << endl;
t=0 E(t) 0.244795
t=0.0005 E(t) 0.0336824
t=0.001 E(t) 0.00455287
t=0.0015 E(t) 0.000610981
t=0.002 E(t) 8.1772e-05
t=0.0025 E(t) 1.09334e-05
t=0.003 E(t) 1.46133e-06
t=0.0035 E(t) 1.95292e-07
t=0.004 E(t) 2.60976e-08
t=0.0045 E(t) 3.48745e-09
t=0.005 E(t) 2.33179e-09
t=0.0055 E(t) 2.28533e-09
t=0.006 E(t) 2.23979e-09

forum7.edp (1.7 KB)

marchywka · May 4, 2023, 11:16am

What do you think of this ( note that this is probably one or more sign mistakes
but if you have a theory you can compare it to that lol). Note
that the result is not too sensitive to dt although ti may still be wrong.
( note I’m just doing this for practice and any feedback would help ).

FreeFem++ forum7x.edp -dt .0005  | grep "t=" | grep " E(t) " | tail
t=0.095 E(t) 0.0322261
t=0.0955 E(t) 0.0315783
t=0.096 E(t) 0.0309436
t=0.0965 E(t) 0.0303217
t=0.097 E(t) 0.0297124
t=0.0975 E(t) 0.0291153
t=0.098 E(t) 0.0285303
t=0.0985 E(t) 0.0279571
t=0.099 E(t) 0.0273955
t=0.0995 E(t) 0.0268452
marchywka@happy:/home/documents/cpp/proj/freefem/play$ FreeFem++ forum7x.edp -dt .0002  | grep "t=" | grep " E(t) " | tail
t=0.098 E(t) 0.0285142
t=0.0982 E(t) 0.028283
t=0.0984 E(t) 0.0280536
t=0.0986 E(t) 0.0278261
t=0.0988 E(t) 0.0276005
t=0.099 E(t) 0.0273767
t=0.0992 E(t) 0.0271547
t=0.0994 E(t) 0.0269345
t=0.0996 E(t) 0.0267161
t=0.0998 E(t) 0.0264995
marchywka@happy:/home/documents/cpp/proj/freefem/play$ FreeFem++ forum7x.edp -dt .001  | grep "t=" | grep " E(t) " | tail
t=0.09 E(t) 0.0394554
t=0.091 E(t) 0.0378915
t=0.092 E(t) 0.0363899
t=0.093 E(t) 0.034948
t=0.094 E(t) 0.0335636
t=0.095 E(t) 0.0322342
t=0.096 E(t) 0.0309577
t=0.097 E(t) 0.029732
t=0.098 E(t) 0.0285551
t=0.099 E(t) 0.0274249

forum7x.edp (2.0 KB)

 diff forum7*
11c11
<      macro grad(u) [dx(u),dy(u)] )  //
---
>      macro grad(u) [dx(u),dy(u)]   //
17c17
< real Tfinal= getARGV("-tfinal",0.01);
---
> real Tfinal= getARGV("-tfinal",0.1);
23a24
> 	Vh ut=u;
30,32c31,42
<         solve Pb(u,v) = int2d(Th) ((1/dt)*u*v) 
<           + int2d(Th) ( [dx(u),dy(u)]'* [dx(v),dy(v)] - c*u*v) 
<           -  int2d(Th) ((1/dt)*up*v) - int2d(Th)(f*v) + on(1,2,3,4, u=g);
---
> 	real cold=c;
> 
> //        solve Pb(u,v) = int2d(Th) ((1/dt)*u*v) 
>  //         + int2d(Th) ( [dx(u),dy(u)]'* [dx(v),dy(v)] - c*u*v) 
>  //         -  int2d(Th) ((1/dt)*up*v) - int2d(Th)(f*v) + on(1,2,3,4, u=g);
> 	solve Pb(u,v) = 	 
>            int2d(Th) (dt*( grad(u)'* grad(v) )) 
>           + int2d(Th) (dt*(  c*u*v)) 
>           -int2d(Th)( ut*v) // TODO check sign lol 
>           +int2d(Th)( u*v)
>   		+ on(1,2,3,4, u=g);
>     c =  sqrt(int2d(Th)(u*u));
38a49
>     ut[]=u[];

aweni · May 4, 2023, 12:23pm

First of all, I would like to thank you for the help you have given me. It encourages me to learn more about FreeFem++.
.
I am not yet used to using the new version of Freefem
so I noticed that you have integrated a lot of things in your code.
My first concern is to know the usefulness of
include “macro_ddm.idp”
and why define the parameters dt and Tfinal in this form.

marchywka · May 4, 2023, 12:32pm

How does the output compare to your expectations?

aweni · May 4, 2023, 12:35pm

I would also like to know how to generate a file containing E(t) and t in order to plot E(t) with the new Freefem ++4.12

marchywka · May 4, 2023, 1:00pm

Do the results look right?

Topic		Replies	Views
Different result than example, could it be due to an installation error? FreeFEM installation	0	296	November 26, 2021
Mesh generation - boundary conditions - command cube Feature Request	1	804	January 14, 2020
Program running General Discussion	12	58	October 11, 2024
About tetrahedron definition General Discussion	3	38	November 21, 2024
Transient heat transfer problem with FreeFem++ General Discussion	1	924	November 18, 2020

mesh of a right paving stone and a cube

Related topics