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+'));