added README_changes.txt

[wrffire.git] / wrfv2_fire / phys / prop_test.m

function prop_test\r
% a jm/mkim version of propagate\r
% follows Osher-Fedkiw and Mitchell toolbox\r
\r
% defaults\r
% grid\r
%xs=15; % domain size\r
xs=10;\r
ys=10;\r
%m=121 % domain mesh\r
m=101\r
n=101\r
time0=0\r
time1=10\r
r=0\r
speed=2  % wind speed\r
wind='c'\r
alpha=0.75*pi\r
dt=0.2;\r
plotSteps=50\r
normal_spread_c=1  % speed = r + wind^e * c\r
normal_spread_e=1\r
\r
example='p'\r
\r
switch example\r
    case {'circlespin','c'}\r
    wind='c'\r
    ignition='c'\r
    case {'gostraight','g'}\r
    wind='s'\r
    ignition='c'\r
    case {'physical','p'}\r
    wind='s';\r
    ignition='c';\r
    alpha=pi;\r
    r=0.02;\r
    m=401\r
    n=401\r
    %m=151\r
    %n=151\r
    speed=10;\r
    xs=6*(m-1);\r
    ys=6*(n-1);\r
    time1=200;\r
    dt=1;\r
    plotSteps=400\r
    normal_spread_c=0.185060861\r
    normal_spread_e=1.310758329\r
    %normal_spread_c=1\r
    %normal_spread_e=1\r
    case {'windmill','w'}\r
    time1=3;\r
    r=0.1;\r
    wind='c'\r
    ignition='w'    \r
    case {'spinstardemo','s'},\r
    % same as spinStarDemo('low',1)\r
    xs=2; % domain size\r
    ys=2;\r
    m=101 % domain mesh\r
    n=101\r
    % normal propagation speed\r
    r=0.2;\r
    % max wind speed\r
    w=1;\r
    wind='c'\r
    ignition='s'\r
    time0=0\r
    time1=1\r
    case {'o','only growth'}\r
    speed=0\r
    wind='s'\r
    ignition='c'\r
    r=0.1\r
end\r
\r
phi=zeros(m,n); % allocate level function\r
\r
dx=xs/(m-1);\r
dy=ys/(n-1);\r
x=dx*[0:m-1];\r
y=dy*[0:n-1];\r
\r
switch wind\r
    case {'straight','s'}\r
        vx=speed*cos(alpha)*ones(m,n);\r
        vy=speed*sin(alpha)*ones(m,n);\r
    case {'circular','c'}\r
        [xx,yy]=ndgrid(x-xs/2,y-ys/2);\r
        rotation=-0.75 * pi;\r
        rotation=0.3 * pi;\r
        vx=yy*rotation;\r
        vy=-xx*rotation;\r
    otherwise\r
        error('no wind')\r
end\r
\r
switch ignition\r
    case {'star','s'}\r
        points = 7;\r
        shift = 2.5;\r
        scale = 0.20;\r
        [xx,yy]=ndgrid(x-xs/2,y-ys/2);\r
        [ theta, rad ] = cart2pol(xx, yy);\r
        phi = rad - scale * (cos(points * theta) + shift);\r
    case {'circle','c'}\r
        c=[0.8*xs,0.5*ys]; % center\r
        d=3*max(dx,dy);      % radius\r
        for i=1:m\r
            for j=1:n\r
%               phi(i,j)=sign(norm([i*dx,j*dy]-c)-d);\r
                phi(i,j)=(norm([i*dx,j*dy]-c)-d);\r
            end\r
        end\r
    case {'widmill','w'}\r
        phi=ones(m,n);\r
        phi(floor(m*0.20):ceil(m*0.45),floor(n*0.49):ceil(n*0.51))=-1;\r
        phi(floor(m*0.55):ceil(m*0.80),floor(n*0.49):ceil(n*0.51))=-1;\r
        phi(floor(m*0.49):ceil(m*0.51),floor(n*0.20):ceil(n*0.45))=-1;\r
        phi(floor(m*0.49):ceil(m*0.51),floor(n*0.55):ceil(n*0.80))=-1;\r
    otherwise\r
        error('no ignition')\r
end\r
\r
\r
%r=ones(m,n).*r + 0.04*sqrt(vx.*vx+vy.*vy); % shrinkage correction\r
r=ones(m,n).*r;  % make sure it is array\r
\r
data2=phi;\r
data3=phi;\r
tNow=time0;\r
dofort=1;\r
domat=0;\r
for i=1:plotSteps\r
    tNext=tNow+dt\r
%    % call the fortran implementation\r
    if dofort,\r
        [data3]=prop_test_f(data3,tNow,tNext,...\r
            normal_spread_c,normal_spread_e,vx,vy,r,dx,dy);\r
    end\r
    %call own imp\r
    %[ux,uy]=get_advection(phi,r,vx,vy,dx,dy);\r
    %vis_wind(ux,uy,dx,dy)\r
    if domat,\r
        [dummy,data2]=prop_ls_cir(data2,tNow,tNext,vx,vy,r,dx,dy,@spread_rate);\r
    end\r
    %[tNow,data2]=prop_ls(data2,tNow,tNext,ux,uy,0,dx,dy);\r
    if dofort & domat, \r
        err_fort=norm(data3(:)-data2(:))\r
        if err_fort>1e-10,\r
            i,warning('large difference between matlab and fortran')\r
        end\r
    end\r
    % display\r
    if dofort,\r
        vis(data3,vx,vy,dx,dy,tNow);\r
    else\r
        vis(data2,vx,vy,dx,dy,tNow);\r
    end\r
    tNow = tNext;\r
end\r
%err_fort=norm(data3(:)-data2(:))\r
% call the toolbox routines \r
% for that need to be in toolboxls/Examples/OsherFedkiw\r
if 0\r
data = normal_advection(phi,time0,time1,vx,vy,r,dx,dy,'low');\r
fprintf('\n');\r
err_toolbox=norm(data(:)-data2(:))\r
end \r
\r
end\r
\r
%-------------------------------------------\r
\r
function vis(u,vx,vy,dx,dy,tNow)\r
[m,n]=size(u);\r
x=[0:m-1]*dx;\r
y=[0:n-1]*dy;\r
figure(1);clf\r
hold off\r
vis_type= '2d';\r
drawn=false;\r
switch vis_type\r
    case '3d'\r
    h=surf(y,x,u);\r
    set(h,'edgecolor','none')\r
    alpha(0.3)\r
    hold on\r
    contour3(y,x,u,[0 0],'k')\r
    hold off\r
    drawn=true;\r
    case '2d'\r
    [m,n]=size(u);\r
    contour(y,x,u,[0 0],'k')\r
    hold on\r
    vis_wind(vx,vy,dx,dy)\r
    axis equal\r
    drawn=true;\r
    otherwise\r
        disp('no visualization selected, set vis=type to 2d or 3d')\r
end\r
if drawn,\r
        title(sprintf('t=%g',tNow));\r
    drawnow;\r
    pause(0.3); % for ctrl-c\r
end\r
end\r
\r
function [phi_out]=prop_test_f(phi,ts,te,c,e,vx,vy,r,dx,dy)\r
fmt='%25.15e \n';\r
if any(any(isnan(phi))), error('phi NaN'), end\r
[m,n]=size(phi);\r
fid = fopen('prop_test_in.txt','wt');\r
fprintf(fid,fmt,m,n,ts,te,c,e,r,dx,dy);\r
fprintf(fid,fmt,phi,vx,vy);\r
fclose(fid);\r
! prop_test_prog.exe\r
fid = fopen('prop_test_out.txt','rt');\r
if fid < 0,\r
    error('cannot open prop_test_out.txt')\r
end\r
[in,c]=fscanf(fid,'%g');\r
fclose(fid);\r
if c~=m*n,\r
    fprintf('read %g terms should have %g\n',c,m*n+1)\r
    error('bad number of terms in output file')\r
end\r
%t=in(1);\r
phi_out=reshape(in(1:c),size(phi));\r
if any(any(isnan(phi_out))), error('phi_out NaN'), end\r
end\r
\r
function speed=spread_rate(r,vx,vy,nvx,nvy,scale);\r
    speed = r + max(vx.*nvx + vy.*nvy,0);\r
end\r
\r