For obtaining a Schlieren image from an equation for density, I need to calculate the first derivative of density and make a contour plot.
The below code snippet calculates the first derivative w.r.t x and y and makes a contour plot. All is well until this step. hSol is an interpolating function polynomial, TRup is the final time of calculation/rupture time in NDSolve and fac is just a value between 0 and 1.
T = 1 - hSol[x, y, fac*TRup]/(hSol[x, y, fac*TRup] + 1);
a = 9.2/10^4;
b = 4.5/10^7;
r = 0.97/(1 + a*T + b*T^2);
ContourPlot[Evaluate[D[r, x] + D[r, y], {x, 0, L}, {y, 0, L}],
PlotRange -> {{0, L}, {0, L}, {0, 3.5}},
BaseStyle -> {FontWeight -> "Bold", FontSize -> 18},
Contours -> 75]
If I replace the dx dy terms with Grad[r] or Laplacian[r], it errors out:
- Either it locks up the RAM and my computer crashes.
- Or All the code in my notebook turns to BLUE as if it were never evaluated in the first place.
What's that all about?
I tried it with a simpler problem as in Example 1 and i could replace y' in the ParametricPlot of In[3] with Laplacian[] and things worked well. (I do realize that y' is not Laplacian...)
*(MY CODE)*Mathematica code that produces the interpolating function polynomial:
(*"Specifies the number of previous lines of input and output to keep in a Mathematica session."*)
$HistoryLength=0;
(*Loads necessary packages to perform vector analysis*)
Needs["VectorAnalysis`"]
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
(*Clears equation symbols*)
Clear[Eq0,EvapThickFilm,h,Bo,\[Epsilon],K1,\[Delta],Bi,m,r] (*Eq0 and EvapThickFilm are variable names. These can be changed as long as the change is made through the entire document*)
(*'h' is the film thickness variable, 'Bo' is the Bond number (G/S), '\[Epsilon]' is the epsilon from the maximizing wavelength equation, 'K1' is the non-equilibrium coefficient, '\[Delta]' is the delta term from the max wavelength equation, 'Bi' is the Biot number, 'm' is the m term and 'r' is the modified Rayleigh number term*)
Eq0[h_,{Bo_,\[Epsilon]_,K1_,\[Delta]_,Bi_,m_,r_}]:=\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]h\)+Div[-h^3 Bo Grad[h]+h^3 Grad[Laplacian[h]]+(\[Delta] h^3)/(Bi h+K1)^3 Grad[h]+m (h/(K1+Bi h))^2 Grad[h]]+\[Epsilon]/(Bi h+K1) + (r)D[D[(h^2/(K1+Bi h)),x] h^3,x] ==0; (*This is the evolution equation in 3D cartesian coordinates*)
SetCoordinates[Cartesian[x,y,z]];
EvapThickFilm[Bo_,\[Epsilon]_,K1_,\[Delta]_,Bi_,m_,r_]:=Eq0[h[x,y,t],{Bo,\[Epsilon],K1,\[Delta],Bi,m,r}];
TraditionalForm[EvapThickFilm[Bo,\[Epsilon],K1,\[Delta],Bi,m,r]];
(*Domain size in non dimensional units*)L=76.1311; (*Maximum time of iteration*)TMax=5000*100;
Off[NDSolve::mxsst];
Clear[Kvar];
(*Kvar[t_]:= Piecewise[{{1,t<=1},{2,t>1}}]*)
(*Ktemp = Array[0.001+0.001#^2&,13]*)
hSol=h/.NDSolve[{
(*Bo,\[Epsilon],K1,\[Delta],Bi,m,r*)
(*Non dimensional number values are fed into the evolution equation which is called here with 'EvapThickFilm[...]'. The arguments, in order are, Bo, \[Epsilon], K1, \[Delta], m, r. These are the exact values used in Oron et al.*)
EvapThickFilm[0.0,0,1,0,1,0.0544911,0],
(*Periodic boundary conditions at the left and right walls in x and y directions*)
h[0,y,t]==h[L,y,t],
h[x,0,t]==h[x,L,t],
(*h[x,y,0] == 1.1+Cos[x] Sin[2y] *)
(*Initial condition as used in Oron et al.*)
h[x,y,0]==1+(-0.05 Cos[2\[Pi] x/L] -0.05 Sin[2\[Pi] x/L])Cos[2 \[Pi] y/L]
},
h,
{x, 0, L},
{y,0, L},
{t, 0, TMax},
Method->{"BDF","MaxDifferenceOrder"->1},
MaxStepFraction->1/50
][[1]]
fac = 0.99;
hGrid = InterpolatingFunctionGrid[hSol];
{TMin, TRup} = InterpolatingFunctionDomain[hSol][[3]]
Length[hGrid]
{nX, nY, nT} = Drop[Dimensions[hGrid], -1]
