rbconvdemo

 script for Rayleigh-Benard
 note that this script already incorporates some knowledge 
 like approximate location of Bif. points and occurrence 
 of imperfect pitchfork, which has been found before (by toying 
 around with pde2path)

0001 % script for Rayleigh-Benard
0002 % note that this script already incorporates some knowledge
0003 % like approximate location of Bif. points and occurrence
0004 % of imperfect pitchfork, which has been found before (by toying
0005 % around with pde2path)
0006 close all; clear all; 
0007 mypause('Init p for Rayleigh-Benard convection for no-slip BC.\n'); 
0008 ns=[];ns=rbconvinit(ns);
0009 mypause('Trivial branch continuation for bifurcation detection.'); 
0010 ns=cont(ns); 
0011 mypause('To speed up demo go directly to point later on branch near 2nd BP.'); 
0012 ns.lam=771; ns.ds=0.1;ns=cont(ns);
0013 % ----------------------------------------------------
0014 mypause('Switch branch at first bifurcation and continue with pmcont'); 
0015 ns1=swibra('ns','bp1','ns1');ns1.ds=2;ns1=pmcont(ns1);
0016 mypause('Switch branch at 2nd bifurcation and continue with pmcont'); 
0017 ns2=swibra('ns','bp2','ns2');ns2.ds=2;ns2=pmcont(ns2); 
0018 % ----------------------------------------------------
0019 mypause('Plot bifurcation diagram:');
0020 figure(3); clf;cmp=1;
0021 plotbra(ns,3,cmp,'lw',6,'ms',10,'cl','b');
0022 plotbra(ns1,3,cmp,'ms',10,'cl','r','lab',10);
0023 plotbra(ns2,3,cmp,'ms',10,'cl','k','lab',8);
0024 axis([739 940 0 3]); xlabel('R'); ylabel('max\psi');
0025 mypause('Quiver plot of last no-slip solution');
0026 arrowplot(ns2,8,1); 
0027 % ----------------------------------------------------
0028 mypause('Initialize with stress-free boundary conditions');
0029 sf=[];sf=rbconvinit_stressfree(sf);
0030 mypause('Seek first bifurcation');
0031 sf.lam=665.5;clf(2);sf=cont(sf);
0032 mypause('Seek 2nd bifurcation');
0033 sf.lam=761;sf=cont(sf); 
0034 mypause('Switch branch at 1st bifurcation and continue with pmcont');
0035 sf1=swibra('sf','bp1','sf1');sf1.ds=2;df1.nsteps=10; sf1=pmcont(sf1); 
0036 % ----------------------------------------------------
0037 fprintf('Switch branch at 2nd bifurcation and continue with pmcont\n');
0038 mypause('but only to close to the imperfect pitchfork');
0039 sf2=swibra('sf','bp2','sf2');sf2.ds=2;sf2.lammax=845;sf2=pmcont(sf2); 
0040 mypause('Extend branch by standard cont');
0041 sf2.ds=1;sf2.dsmax=2;sf2.lammax=900;sf2.nsteps=50;sf2=cont(sf2);
0042 % ----------------------------------------------------
0043 mypause('Switch to disconnected branch by time-integration');
0044 p=sf2;p=resetc(p); r=resi(p,p.u,p.lam); Gu=getGu(p,p.u,p.lam,r); 
0045 [K,MId]=assema(p.points,p.tria,0,1,zeros(p.neq,1)); 
0046 [phi1,mu1]=eigs(Gu,MId,1,0,p.evopts);
0047 p1=p;p1.u=p.u+0.01*phi1/norm(phi1);
0048 p1=tint(p1,0.1,400,50); 
0049 mypause('Now continue from p1'); 
0050 sf3=p1; sf3.restart=1; sf3.ds=-sf3.ds; sf3.lammax=sf3.lammax+2; 
0051 sf3.bifchecksw=0;sf3.dlammax=10;sf3.nsteps=100;sf3=cont(sf3);
0052 % ----------------------------------------------------
0053 mypause('Plot bifurcation diagram'); 
0054 clf(3);cmp=1;
0055 plotbra(sf,3,cmp,'lw',5,'ms',10,'cl','b');
0056 plotbra(sf1,3,cmp,'ms',10,'cl','r','lab',40,'inc',40); 
0057 plotbra(sf2,3,cmp,'ms',10,'cl','k','lab',[30 60]);
0058 plotbra(sf3,3,cmp,'ms',10,'cl','b','lab',65,'inc',65);
0059 axis([660 920 0 3]);xlabel('R'); ylabel('max\psi');
0060 mypause('Plot some solutions'); 
0061 p1=loadp('sf1','p40','p1'); p2=loadp('sf2','p20','p2'); 
0062 p3=loadp('sf2','p60','p3'); p4=loadp('sf3','p65','p4'); 
0063 arrowplot(p1,8,1); arrowplot(p2,9,1); 
0064 arrowplot(p3,10,1); arrowplot(p4,11,1);