1//ff-mpirun  -np 4 Richards-2d.edp -wg -raspart -ns -ffddm_disable_plots
2
3// If you have openmpi you may need to add the option --oversubscribe to allow more processes than the number of cores available on your computer
4
5// for the make check:
6// NBPROC 4
7// PARAM -raspart -ffddm_disable_plots -nt 20
8
9macro dimension 2// EOM            // 2D or 3D
10
11include "ffddm.idp"
12
13int nt = getARGV("-nt", 100);
14
15searchMethod = 1; // safe search algo, use brute force in case of missing point
16
17macro def(i)i// EOM                         // scalar field definition
18macro init(i)i// EOM                        // scalar field initialization
19func Pk = P1;
20
21real Ks=0.01,
22     hg=30,
23     thetas=0.3,
24     eta = 6.55,
25     m = 0.173,
26     n = 2/(1-m);
27
28
29real z0=215;
30real q0=15/3600.;
31real dt=60;
32
33// $A(h) - \partial h / \partial t - div(K(h)(\nabla(h-y)) = f $ dans $ \Omega$
34//  -K(h)(\nabla(h-y)). n = q0 $ sur $ \Gamma_1$
35//  h = h_0$ sur $\Gamma_0$
36//  + condition initial $d_d$
37//   A(h) = h <0 ? C(h) : 0;
38//
39
40
41// remarque z == y
42real xmax = 300, ymax=300, x0=60, y0= 215;
43
44border ba(t=0,ymax)   { x=0; y=ymax-t ;label=2;};  // left
45
46border bb1(t=x0,0)    { x=t; y=ymax ;label=1;};   // top     1
47border bb2(t=xmax,x0) { x=t; y=ymax ;label=2;};  // top   2
48
49border bc1(t=y0,0) { x=xmax; y=ymax-t ;label=2;};  // right
50border bc2(t=ymax,y0) { x=xmax; y=ymax-t ;label=3;};  // right
51
52border bd(t=0,xmax)   { x=t; y=0; label=4;};   // bottom
53
54int Gamma0=3;
55int Gamma1=1;
56
57int nn=20;
58int nn1=nn*x0/xmax,nn2=nn-nn1;
59int ny1=nn*y0/ymax,ny2=nn-ny1;
60plot(ba(nn)+bb1(nn1)+bb2(nn2)+bc1(ny1)+bc2(ny2)+bd(nn),wait=1);
61mesh Th=buildmesh(ba(nn)+bb1(nn1)+bb2(nn2)+bc1(ny1)+bc2(ny2)+bd(nn));
62plot(Th,wait=1);
63
64mesh Thc = Th;
65
66fespace Vh(Th,P1);
67Vh h,v,hhh;
68
69macro theta(h) (thetas*(1+((abs(h)/hg))^n)^(-m))//
70macro dtheta(h) (m*n*thetas*(1+((abs(h)/hg))^n)^(-m-1)*(((abs(h)/hg))^(n-1))/hg)
71//
72macro A(h)  ( (h<=0)* dtheta(h) )
73//
74macro K(h) (Ks*((h<=0)*((theta(h)/thetas)^eta)+ (h>0)))
75//
76
77Vh hd= -y0+(ymax-y); // bof bof ????
78Vh hn=hd,hh;
79Vh Ahdt,Kh;
80
81int nbiso=20;
82real[int] viso(3+(75+110/2)/5);
83
84{int k=0;
85for(int i=-75;i<0;i+=5)
86 viso(k++)=i;
87 viso(k++)=-0.5;
88 viso(k++)=0.;
89 viso(k++)=0.5;
90for(int i=5;i<=110;i+=5*2)
91 viso(k++)=i;
92}
93/*
94problem Richard(h,v) =
95  int2d(Th)( Ahdt * h * v+ Kh* (dx(h)*dx(v)+dy(h)*dy(v)) )
96- int2d(Th)( Ahdt* hn*v - Kh* dy(v) )
97- int1d(Th,Gamma1)(q0*v)
98+ on(3,h=(ymax-y)-y0)
99;
100*/
101
102real pena=1e10;
103
104macro VarfPb(varfName, meshName, VhName)
105
106  VhName mh, mhn;
107  RfromVhi(hi[],VhName,mh[])
108  RfromVhi(hin[],VhName,mhn[])
109  VhName Ahdt = A(mh)/dt;
110  VhName Kh = K(mh);
111
112  varf varfName(h,v) =
113  int2d(meshName)( Ahdt * h * v+ Kh* (dx(h)*dx(v)+dy(h)*dy(v)) )
114- int2d(meshName)( Ahdt* mhn*v - Kh* dy(v) )
115- int1d(meshName,Gamma1)(q0*v)
116 +int1d(meshName,3)(pena*h*v)-int1d(Th,3)(pena*((ymax-y)-y0)* v)
117
118; // EOM
119
120plot(hn,wait=1,cmm=" hd ");
121
122macro Rdefmplot(u)u//EOM
123
124//macro Rwithhpddmkrylov()1//EOM
125
126ffddmbuildDmeshAug(R,Th,mpiCommWorld)
127
128ffddmbuildDfespace(R,R,real,def,init,Pk)
129
130ffddmsetupinit(R,R)//
131
132ffddmset(R,verbosity,1)//
133
134RVhi hi=hd, hin=hd;
135
136real[int] rhs;
137
138// Richard;
139// plot(hd,wait=1,cmm="hd ----");
140real temps=0;
141for(int ii=0;ii<nt;ii++)
142{
143  string scmm="h + temps "+int(temps)/3600+"h "+ ((temps)%3600/60.) + "mn ";
144  for(int k=0;k<3;k++)
145  {
146
147	ffddmsetup(R,R,VarfPb,null)
148
149	ffddmbuildrhs(R,VarfPb,rhs)
150
151	rhs = -rhs;
152
153  //cout << " "<< Kh[].min << " " << Kh[].max << endl;
154 // plot(Ahdt,fill=1,value=1,wait=1,cmm="Ahdt");
155 // plot(Kh,fill=1,value=1,wait=1,cmm="Kh");
156 // plot(Kh,wait=1,cmm="Kh");
157
158	hi[] = RfGMRES(hi[], rhs, -1.e-6, 200, "right");
159
160  // forces a refresh of the FE function h, so that it is defined on the currrent mesh Th
161  // in particular, the underlying array h[] has the correct size
162  h = 0.;
163
164	RfromVhi(hi[],Vh,h[]);
165
166  // plot(h,wait=0,cmm=scmm,viso=viso);
167//  hhh = h <0;
168  }
169  if(ii%10==1) {
170  Th=adaptmesh(Th,h,ratio=1.1);
171  broadcast(processor(0),Th);
172
173  // the mesh has changed ; we create a new decomposition and obtain new local meshes
174  ffddmbuildDmeshAug(R,Th,mpiCommWorld)
175  ffddmbuildDfespace(R,R,real,def,init,Pk)
176  // interpolate the solution at the previous time step on the new local mesh
177  hi = h;
178
179//  plot(Th,h,cmm="h ",value=1,wait=1);
180  }
181 // plot(hhh,wait=1,cmm="h < 0");
182
183  plot(h,wait=0,cmm=scmm,viso=viso,value=1);
184//  plot(h,cmm="h ",value=1);
185  hin = hi;
186  temps += dt;
187}
188
189Rwritesummary
190