Monday, 2 October 2017

Belonging, purpose, transcendence and storytelling...

Navier-Stokes, good luck making the code shorter...


N  = 192; dt = 0.25; p = zeros(N, N);
u = zeros(N, N); v = zeros(N, N); ustar = zeros(N, N); vstar = zeros(N, N);
k = pi*[0:(N/2-1) (-N/2):(-1)]; [x y]  = meshgrid(k, k); ds = -(x.^2 + y.^2);
ds(1,1) = 1; px = 0.0; py = 0.0; gx = 0.0; gy = 0.0;
for t = 0:1000
    % constant velocity zone
    u(N/2-10:N/2+10, N/2-10:N/2+10) = 0.2;

    % solve intermediate velocity U*
    [ux uy] = gradient(u);
    [vx vy] = gradient(v);
    ustar = u + (u.*ux + v.*uy + del2(u)/20 - px/20 + gy)*dt;
    vstar = v + (u.*vx + v.*vy + del2(v)/20 - py/20 + gx)*dt;

    % solve Poisson equation: del2(p) = f
    fhat = fft2(divergence(ustar, vstar));
    p = N^2*real(ifft2(fhat./ds));

    % enforce del(U) = 0
    [px py] = gradient(p);
    u = ustar - px*dt;
    v = vstar - py*dt;

    cla; surface(sqrt(u.^2 + v.^2), 'edgecolor', 'none'); drawnow;
end

Espoo is full of furry animals.
Just a basic wave equation...


u  = zeros(192, 192);
du = zeros(192, 192);
[x y] = meshgrid(linspace(-1, 1, 192), linspace(-1, 1, 192));
for t = 0:0.002:0.4
    u = u + exp(-1000*((x-0.7*sin(7*t)).^2+(y-0.7*sin(11*t)-0.2).^2));
    u = u + 0.5*exp(-1000*((x+0.7*sin(5*t)).^2+(y+0.7).^2))*sin(3*31*t);

    du = du + del2(u);   
    u = u + du;
       
    cla; imagesc(u); drawnow;
end


...

No comments:

Post a Comment