vgpinit

 init for vector GP

0001 function p=vgpinit(p) 
0002 % init for vector GP
0003 p=stanparam(p); p.neq=4; p.f=@vgpf; p.jac=@vgpjac; p.outfu=@vgpbra; 
0004 pre=sprintf('%s',inputname(1)); p=setfn(p,pre);
0005 lx=5; ly=1*lx; p.tol=1e-6; [p.geo,p.bc]=recnbc2(lx,ly); 
0006 sc=100; qd=[sc, sc, sc, sc]; qm=diag(qd,0); % stiff spring Dirichlet
0007 p.bc=gnbc(4,4,qm,zeros(p.neq,1)); p.bcf=@(p,u,lam) p.bc;
0008 nx=30;ny=nx;p=stanmesh(p,nx,ny);p=setbmesh(p); 
0009 p.dsmin=0.001; p.dsmax=0.05; p.dlammax=0.2; p.nsteps=20; p.lammax=1.5; 
0010 p.jsw=0; p.bifchecksw=0; p.spcalcsw=0;p.pcmp=1; 
0011 p.imax=10; p.bpcmp=3; p.parasw=0; p.smod=5;p.pstyle=3; p.cm=cool; 
0012 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% starting point %%%%%%%%%%%%%%%%%%%%%%%%%%%%
0013 p.lam=0; p.ds=0.1; p.pa=1; p.xi=1/p.np; p.mu1=2; p.mu2=2.2; 
0014 x=p.points(1,:)'; y=p.points(2,:)'; p.r=x.^2+y.^2; p.phi=angle(x+1i*y); 
0015 m=1; n=0; % first dipole (i.e. dipole in u1,v1)
0016 am=pi*factorial(m+1)/2; gm=am; bm=pi*factorial(2*m)/2^(2*m+5); 
0017 dm=pi*factorial(m)/2; a2=(sqrt(dm^2*p.mu1^2+12*am*gm)+dm*p.mu1)/(6*gm); 
0018 a=sqrt(a2);A2=(2*a2*gm-dm*p.mu1)/(bm*3); A=sqrt(A2); 
0019 u=A*tfu(p.r/a,m,n).*cos(m*p.phi); v=0.0*u; u0=[u v]; 
0020 m=1; n=0; % second dipole (i.e. dipole in u2,v2)
0021 am=pi*factorial(m+1)/2; gm=am; bm=pi*factorial(2*m)/2^(2*m+5); 
0022 dm=pi*factorial(m)/2; a2=(sqrt(dm^2*p.mu2^2+12*am*gm)+dm*p.mu2)/(6*gm); 
0023 a=sqrt(a2); A2=(2*a2*gm-dm*p.mu2)/(bm*3); A=sqrt(A2); 
0024 u=A*tfu(p.r/a,m,n).*cos(m*p.phi); v=0.0*u; 
0025 u0=[u0 u v]; p.u=reshape(u0,p.neq*p.np,1); % append to first and reshape
0026 plotsol(p,1,1,3); p.tau=zeros(p.neq*p.np+1,1); p.tau(1)=1; 
0027 p.ngen=10;p.adafac=1.5;p.maxt=p.adafac*p.nt; p.xi=1/p.np; 
0028 res=norm(resi(p,p.u,p.lam), p.normsw); fprintf('initial res=%g\n',res);
0029 [p.u,res]=nlooppde(p,p.u,p.lam); fprintf('first res=%g\n',res); 
0030 plotsol(p,1,1,3); p=meshref(p,'maxt',6000,'eb',0); 
0031 p=setbmesh(p); p.maxt=9000; % set maxt for later adaption
0032 plotsol(p,p.pfig,p.pcmp,p.pstyle); p.mst=10;
0033 [p.u,res]=nlooppde(p,p.u,p.lam); fprintf('first res=%g\n',res); 
0034 
0035