addpath('../../SOURCES_MATLAB/'); SF_Start('verbosity',4);
NOTICE - Initialization already done WARNING - Detected verbosity = 4 in Publish mode ; not recommended ! switching to 2 instead
Generate a mesh and a starting flow
SF_core_setopt('verbosity',2); ffmesh = SF_Mesh('Mesh_ConvectionCell.edp','Params',[4 20 5],'problemtype','2dboussinesq');
WARNING - Your program uses depreciated syntax 'Params' to pass parameters WARNING - It is advised to switch to new method using 'Options', and modify your Freefem solvers to detect parameters using 'getARGV'
compute a "base flow" for Ra = 2000
bf = SF_BaseFlow(ffmesh,'Ra',2000,'Pr',0.03,'Qmag',0);
WARNING - Folder ./WORK/MESHBF not found : creating it !
look for eigenmodes
[ev,em] = SF_Stability(bf,'k',0,'shift',5,'nev',10,'sort','lr'); % plot the leading eigenmode figure; SF_Plot(em(1),'T'); hold on; SF_Plot(em(1),{'ux','uy'},'fgridparam',[40 10],'xlim',[0 4],'title','Leading eigenmode for Ra = 2000') pause(0.1);
WARNING - It is now advised to generate a file Spectrum.ff2m in your stability solvers (see examples, e.g. Stab_2D.edp) ans = Quiver with properties: Color: [0 0 1] LineStyle: '-' LineWidth: 0.5000 XData: [400x1 double] YData: [400x1 double] ZData: [] UData: [400x1 double] VData: [400x1 double] WData: [] Use GET to show all properties
We construct a "guess" as a superposition of the previous bf + the unstable eigenmode
amp = 5; % amplitude (max. velocity) of the expected pattern. NB this coefficient has to be tuned by trial-error guess = SF_Add({bf,em(1)},'coefs',[1 amp]);
WARNING - in SF_Add : could not process with auxiliary field vort WARNING - in SF_Add : could not process with auxiliary field psi
compute a steady state using this "guess" as a starting point
bfNL = SF_BaseFlow(guess,'Ra',2000,'Pr',0.03,'Qmag',0);
plot this nonlinear state
figure; SF_Plot(bfNL,'T'); hold on; SF_Plot(bfNL,{'ux','uy'},'fgridparam',[40 10],'xlim',[0 4],'title','Nonlinear state for Ra = 2000') pause(0.1);
ans = Quiver with properties: Color: [0 0 1] LineStyle: '-' LineWidth: 0.5000 XData: [400x1 double] YData: [400x1 double] ZData: [] UData: [400x1 double] VData: [400x1 double] WData: [] Use GET to show all properties
If you want the value of the normalized heat flux and the max velocity (norm) it is available in this way
bfNL.HeatFlux bfNL.Umax
ans = 1.0672 ans = 2.6961
look for 3D eigenmodes
[ev,em] = SF_Stability(bfNL,'k',pi,'shift',10+10i,'nev',20,'sort','lr','plotspectrum','yes'); figure; SF_Plot(em(1),'T','title','Leading 3D eigenmode of 2D-convection state for Ra = 2000') hold on; SF_Plot(em(1),'T.im','contour','only')
WARNING - It is now advised to generate a file Spectrum.ff2m in your stability solvers (see examples, e.g. Stab_2D.edp) ans = 2x1 graphics array: ColorBar Patch ans = 2x1 Contour array: Contour Contour
Appendix
if(1==0)
do a loop over Ra to investigate the stability of this nonlinear state with respect to 2D perturbations
RaTAB = [2000 :200 : 4000] for j=1:length(RaTAB) Ra = RaTAB(j); bfNL = SF_BaseFlow(bfNL,'Ra',Ra); ev = SF_Stability(bfNL,'k',0,'shift',10+10i,'nev',10,'sort','lr'); evTAB(j) = ev(1) end figure; plot(RaTAB,evTAB);
end % [[PUBLISH]] % fclose(fopen('SCRIPT_RayleighBenard.success','w+'));