Reza Malek-Madani

stommel1

tau0=10000000000; a=6000; b=4000; beta=.00000002;rho0=10^15;
tau0x=-tau0 Cos[Pi a y/b];
k=1/beta D[tau0x/rho0,y];eps=.025;
psi=k (Exp[-x/eps]+x-1);
ContourPlot[psi,{x,0,1},{y,0,b/a},PlotPoints->200,
  PlotLabel -> "Eps =0.025"];
NumOfPoints=100;     (*Number of points on the boundary of parcels*)
radius = 0.031;      (*Radius of each initial parcel*)
(* Initial Centers of individual parcels*) 
centers={{0.8, 0.3}, {0.8, 0.4}, {0.8, 0.5}, {0.5, 0.3}, 
    {0.5, 0.4}, {0.5, 0.5}}; 
tfin=0.0005;             (*Total Time*)
interval = 0.00001;       (*Time between snapshots*)
diffeqn[tfinal_,aa_,bb_]:=NDSolve[{x'[t]==f[x[t],y[t]],
                                 y'[t]==g[x[t],y[t]],
                                 x[0]==aa, y[0]==bb},
                                 {x,y}, {t,0,tfinal}, MaxSteps->10000];
f[x_,y_] = D[psi,y]; g[x_, y_] = -D[psi, x];

Do[initdata[i]=N[Table[{centers[[i,1]]+radius*Cos[w], 
                        centers[[i,2]]+radius*Sin[w]}, 
        {w, 0, 2*Pi, 2*Pi/NumOfPoints}]];
     oldsolution[i]=Table[diffeqn[tfin,initdata[i][[j,1]],
                           initdata[i][[j,2]]], {j, Length[initdata[i]]}];
     data[i]=Table[{x[t], y[t]} /. oldsolution[i], {t, 0, tfin, interval}];
     ColoredSnapshots[i]=Table[Graphics[{RGBColor[1-0.2*i, 0.2*i, 0],
     Map[Polygon, Flatten[data[i][[j]],1], {0}]}],  {j, Length[data[i]]}],
{i,Length[centers]}];
graph1=ContourPlot[psi, {x, 0, 1}, {y, 0, b/a}, PlotPoints->50, 
    Contours->20];
Do[
    giffile=Show[graph1,Table[ColoredSnapshots[k][[i]],
    {k,Length[centers]}],
    PlotRange->{{0,1}, {0, b/a}},
    AspectRatio->Automatic];
    llabel=StringJoin[ToString[i],".gif"];
    Display[llabel,giffile,"GIF"],
{i, Length[ColoredSnapshots[1]]}];
go to Top