// run with MPI: ff-mpirun -np 4 script.edp // NBPROC 4 // // code borrowed from examples/ffddm/heat-torus-3d-surf.edp /* @software{mjm_cpp_$PROJ, author = {Michael J Marchywka}, title = {$PROJ}, abstract=(), institution={}, license={Knowledge sir should be free to all }, publisher={Mike Marchywka}, email={marchywka@hotmail.com}, authorid={orcid.org/0000-0001-9237-455X}, filename = {$PROJ}, url = {}, version = {0.0.0}, date-started = {$DATE} } */ // finally had to remove all that junk //load "hpddm" // HPDDM plugin load "medit" // HPDDM plugin //macro dimension()2// EOM // 2D, 3D, or 3S include "macro_ddm.idp" // additional DDM functions //macro def(i)i// EOM // scalar field definition // syntax errors due to forgotten macro are a huge problem // and not clear from eror message !!!!!!!!!!!!!!!!!!!! //macro init(i)i// EOM // scalar field initialization macro Grad(u)[dx(u), dy(u)]// EOM // two-dimensional gradient func Pk = P2; // finite element space int s = getARGV("-split", 1); // refinement factor int sz = getARGV("-sz", 50); // refinement factor int sz2 = getARGV("-sz2", 10); // refinement factor real dt = getARGV("-dt", 1e-4); // time step int iMax = getARGV("-iMax", 50); // number of iterations int nx =getARGV("-nx", 30); int ny = getARGV("-ny", 30); cout<<" sz "< mat; // local operator real f1z= getARGV("-f1z", 100); // real f1az= getARGV("-f1az", 0); // real f2z= getARGV("-f2z", 100); // real f2az= getARGV("-f2az", 0); // real fgnz= getARGV("-fgnz", 100); // real fbz= getARGV("-fbz", 0); // real tfz= getARGV("-tfz", 20); // real c2= getARGV("-c2", .001); // real cfgn= getARGV("-cfgn", .001); // real df1= getARGV("-df1", 1); // real df1a= getARGV("-df1a", .001); // real df2= getARGV("-df2", 1); // real df2a= getARGV("-df2a", .0001); // real dfgn= getARGV("-dfgn", 1); // real dfb= getARGV("-dfb", .001); // real dzed= getARGV("-dzed",6000); // real edfb= getARGV("-edfb",100); // real cdfb= getARGV("-cdfb",.1); // func real tfx( real xx, real yy) { //return 1000; //cout<= -sz2)) if ((yy<=sz2)&&(yy>= -sz2)) return tfz; return 0; } Vh tf= tfx(x,y); Wh [f1,f1a,f2,f2a,fgn,fb]=[f1z,f1az,f2z,f2az,fgnz,fbz]; Wh [f1w,f1aw,f2w,f2aw,fgnw,fbw]=[f1z,f1az,f2z,f2az,fgnz,fbz]; Wh [vf1,vf1a,vf2,vf2a,vfgn,vfb]; Vh D= dzed/(1+cdfb*exp(fb*edfb)) ; problem vDiff([f1,f1a,f2,f2a,fgn,fb], [vf1,vf1a,vf2,vf2a,vfgn,vfb],solver=GMRES,init=1) = //problem vDiff(f1,f1a,f2,f2a,fgn,fb, vf1,vf1a,vf2,vf2a,vfgn,vfb,solver=CG,init=j) = int2d(Th) ( +dt*(df1*WVDIFF(D,f1,vf1)+df1a*WVDIFF(D,f1a,vf1a)+df2*WVDIFF(D,f2,vf2) +df2a*WVDIFF(D,f2a,vf2a)+dfgn*WVDIFF(D,fgn,vfgn)+dfb*WVDIFF(D,fb,vfb) ) +(f1*vf1+f1a*vf1a+f2*vf2+f2a*vf2a+fgn*vfgn+fb*vfb) ) -int2d(Th)(f1w*vf1+f1aw*vf1a+f2w*vf2+f2aw*vf2a+fgnw*vfgn+fbw*vfb) +int2d(Th)(dt*f1w*vf1*tf) +int2d(Th)(-dt*f1w*vf1a*tf) +int2d(Th)(dt*f2w*vf2*f1aw*c2) +int2d(Th)(-dt*f2w*vf2a*f1aw*c2) +int2d(Th)(dt*fgnw*vfgn*f2aw*cfgn) +int2d(Th)(-dt*fgnw*vfb*f2aw*cfgn) //+on(1,2,3,4,[f1,f1a,f2,f2a,fgn,fb]=[f1z,f1az,f2z,f2az,fgnz,fbz]) +on(1,2,3,4,f1=f1z) +on(1,2,3,4,f1a=f1az) +on(1,2,3,4,f2=f2z) +on(1,2,3,4,f2a=f2az) +on(1,2,3,4,fgn=fgnz) +on(1,2,3,4,fb=fbz) ; for(int i=0; i<10000; ++i) { vDiff; [f1w,f1aw,f2w,f2aw,fgnw,fbw]=[f1,f1a,f2,f2a,fgn,fb]; //medit("asdf",Th,fb); macro myplot()cmm = "Global solution at iteration " + i, fill = 1, value = 1// plot(fb,myplot); if ((i%1000)==9){ medit("f1",Th,f1); medit("f1a",Th,f1a); medit("f2",Th,f2); medit("f2a",Th,f2a); medit("fgn",Th,fgn); medit("fb",Th,fb); medit("D",Th,D); } D= dzed/(1+cdfb*exp(fb*edfb)) ; } // i