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.
Juste a small example
aweni.edp (582 Bytes)
with a fixed point iteration,
remark, the you can do Newton iteration, in 10 mn.
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;
// 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;
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;
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.
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)
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 ;
}
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
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.
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?
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)
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.
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)
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[];
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.
How does the output compare to your expectations?
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
Do the results look right?