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