Hi,
i have to calculate in 3D solutions u and w of two problems:
du/dt - \Delta u = f and dw/dt - Delta w +w =f , (x,y,z) \in (0,1)^3, t \in [0,1]
where f(x,y,z)= 0.25*sin(2*pi*x)
and calculate the norme H^1 of the difference between u and w.
There is what i try. Is it correct?
Kind regards
’ ’ ’ freefem
load “UMFPACK64”
defaulttoUMFPACK64();
load “iovtk”
load “msh3”
verbosity=1;
func real f (real x)
{
return 0.25sin(2pi*x);
};
real t=0.00;
real dt=0.1;
mesh3 Th=cube(20,20,20);
fespace Vh(Th,P1,periodic=[[1,x,z],[3,x,z],[2,y,z],[4,y,z],[5,x,y],[6,x,y]]);
fespace Wh(Th,P2,periodic=[[1,x,z],[3,x,z],[2,y,z],[4,y,z],[5,x,y],[6,x,y]]);
Vh uh=0.0;
Vh wh=0.0;
Vh vh;
Vh oldU=0.0;
Vh oldW=0.0;
real T=1.0;// arbitraire
int[int] Order= [1];
Vh fh=f;
// Nonlinear part is calculated in previous iteration (nith) and becomes a part of the right hand side
problem het(uh,vh,init=t) =
int3d(Th)((dx(uh)dx(vh)+dy(uh)dy(vh)+dz(uh)dz(vh))dt)
+int3d(Th)(uhvh)
-int3d(Th)(oldUvh)
-int3d(Th)(fhvhdt)
;
problem hom(wh,vh,init=t) =
int3d(Th)((dx(wh)dx(vh)+dy(wh)dy(vh)+dz(wh)dz(vh))dt)
+int3d(Th)(whvh)
-int3d(Th)(oldWvh)
+int3d(Th)(whvhdt)
-int3d(Th)(fhvhdt)
;
for (t=0;t<=T;t+=dt) {
fh=f(t);
oldU=uh;
het;
}
for (t=0;t<=T;t+=dt) {
fh=f(t);
oldW=wh;
hom;
}
real errH1 = sqrt(int3d(Th)(square(uh - wh)) + int3d(Th)(square(dx(uh)-dx(wh)) + square(dy(uh)-dy(wh))+ square(dz(uh)-dz(wh))));
cout << " errH1 = "<< errH1 <<endl;