matlab画图基础程序

%% Brenner_scheme_2D_for_Euler %% computational domain omega=[0,1]^2 %======================================================================= %======================================================================= function nf = MC_KH() maxNumCompThreads(90) alpha = 0.8; beta = 0.0; mesh = [64 128 256 512 1024 2048]; Nmesh = 1; Nsample = 100; parfor isample = 1:100 brenner_kh(2048, alpha, beta, isample); end nf=1; end %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ function nf = brenner_kh(nx, alpha, beta, label) %clear all; clc; close all; format long; %----------------------------------------------------------------------- % basic parameters lx = 1.0; ly = 1.0; t_max = 2.0; gamma = 1.4; CFL = 0.1; epsilon = 0.01; %nx=32; ny = nx; x = linspace(0, lx, nx + 1); y = linspace(0, ly, ny + 1); h = lx / nx; cfl_hx = CFL * h; p_flag = 0; MESH = num2str(nx); %alpha=1.2; %% alpha should be (0,2) %beta=0.2; if (alpha == -1 && beta == -1) ha = 0.0; mu = 0.0; DIR = strcat('./DataKH/upwind_', MESH, '_'); else ha = h ^ (alpha - 1.0); mu = h ^ beta; S1 = num2str(alpha); S2 = num2str(beta);S3=num2str(label); DIR = strcat('./DataKH/grid',MESH,'/', S1, '_', S2, '_', MESH, '_', S3,'_'); end t = 0.0; nt = 0; %% Initialization [rho, mx, my, U, V, ene, pre, tem] = initial_kelvin_helmholtz_rand(nx, ny, h, gamma, epsilon); % subplot(2,2,1); surf(rho','LineStyle','none'); xlabel('x');ylabel('y');title('\rho');colorbar; view(0,90); % subplot(2,2,2); surf(U','LineStyle','none'); xlabel('x');ylabel('y');title('u');colorbar;view(0,90); % subplot(2,2,3); surf(V','LineStyle','none'); xlabel('x');ylabel('y');title('v');colorbar;view(0,90); % subplot(2,2,4); surf(pre','LineStyle','none'); xlabel('x');ylabel('y');title('p');colorbar;view(0,90); % pause; dtmin = 1; fprintf('nx=%d\n', nx); while (t < t_max) nt = nt + 1; if (max(max(tem)) < 0.0) error('fatal error, negative temperature,nt=%d,t=%.2e', nt, t); end dt = get_time_step(gamma, cfl_hx, U, V, tem, ha, mu, h); % % 2D good! % dtmin=min(dt,dtmin); % dt=min(dt,dt0); t = t + dt; rhoc = rho(2:end - 1, 2:end - 1); mxc = mx(2:end - 1, 2:end - 1); myc = my(2:end - 1, 2:end - 1); prec = pre(2:end - 1, 2:end - 1); enec = ene(2:end - 1, 2:end - 1); Uc = U(2:end - 1, 2:end - 1); Vc = V(2:end - 1, 2:end - 1); Ua = avg(U(:, 2:end - 1)); Va = avg(V(2:end - 1, :)')'; Px = diff(avg(pre(:, 2:end - 1))) / h; Py = diff(avg(pre(2:end - 1, :)'))' / h; divu = discrete_divergence(h, Ua, Va); ue = 0.5 * (U .* U + V .* V); [flux_rho, flux_mx, flux_my, flux_ene] = get_upwindflux(h, rho, mx, my, ene, Ua, Va); % % the fluxes has the size (nx,ny) % rhoc = rhoc + dt * (- flux_rho + arti_diff(mu, h, rho)); mxc = mxc + dt * (-flux_mx - Px + arti_diff(ha, h, U) + arti_diff(mu, h, mx)); myc = myc + dt * (-flux_my - Py + arti_diff(ha, h, V) + arti_diff(mu, h, my)); %enec = enec + dt * (-flux_ene - Px .* Uc -Py .* Vc - divu .* prec + arti_diff(ha, h, ue) + arti_diff(mu, h, ene)); % rho = apply_boundary_condition(rhoc); mx = apply_boundary_condition(mxc); my = apply_boundary_condition(myc); ene = apply_boundary_condition(enec); [U, V, pre, tem] = update_other_values(gamma, rho, mx, my, ene); if floor(25 * t / t_max) > floor(25 * (t - dt) / t_max), fprintf('.'), end end fprintf('Minimum time step due to CFL condition: dtmin=%f\n', dt); fprintf('Computation finished for n=%d\n', nx); nf = 1; ent = rho .* (2.5 * log(pre) - 3.5 * log(rho)); filename = strcat(DIR, 'rho.txt'); dlmwrite(filename, rho, 'precision', 16); %filename=strcat(DIR,'U.txt'); dlmwrite(filename, U,'precision',16); %filename=strcat(DIR,'V.txt'); dlmwrite(filename, V,'precision',16); %filename=strcat(DIR,'P.txt'); dlmwrite(filename, pre,'precision',16); filename = strcat(DIR, 'm1.txt'); dlmwrite(filename, mx, 'precision', 16); filename = strcat(DIR, 'm2.txt'); dlmwrite(filename, my, 'precision', 16); %filename = strcat(DIR, 'E.txt'); dlmwrite(filename, ene, 'precision', 16); filename = strcat(DIR, 'S.txt'); dlmwrite(filename, ent, 'precision', 16); if (p_flag == 1) filename = strcat(DIR, 'rho.png'); fa = avg(avg(rho)')'; myplot(x, y, fa, filename); filename = strcat(DIR, 'P.png'); fa = avg(avg(pre)')'; myplot(x, y, fa, filename); filename = strcat(DIR, 'U.png'); fa = avg(avg(U)')'; myplot(x, y, fa, filename); filename = strcat(DIR, 'V.png'); fa = avg(avg(V)')'; myplot(x, y, fa, filename); end % % rhoc=rho(2:end-1,2:end-1); % figure;contour(avg(x),avg(y),rhoc',29,'LineWidth',2); reset_fig; axis equal;%colorbar; % filename=strcat('./Data2D/rho'); % saveas(gcf,filename, 'png'); % % figure;contour(avg(x),avg(y),U(2:end-1,2:end-1)',30,'LineWidth',2); reset_fig; axis equal; %colorbar; % filename=strcat('./Data2D/u'); % saveas(gcf,filename, 'png'); % % figure;contour(avg(x),avg(y),V(2:end-1,2:end-1)',30,'LineWidth',2); reset_fig; %colorbar; % axis equal; % filename=strcat('./Data2D/v'); % saveas(gcf,filename, 'png'); % % figure;contour(avg(x),avg(y),pre(2:end-1,2:end-1)',30,'LineWidth',2); reset_fig;axis equal;%colorbar; % filename=strcat('./Data2D/p'); % saveas(gcf,filename, 'png'); % % figure;contour(avg(x),avg(y),pre(2:end-1,2:end-1)',29,'LineWidth',2); reset_fig;axis equal;%colorbar; % filename=strcat('./Data2D/p2'); % saveas(gcf,filename, 'png'); end %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ function dt = get_time_step(gamma, cfl_hx, U, V, tem, ha, mu, h) % % 2D good! um = sqrt(U .* U + V .* V); vel_max = max(max(um)); c_max = sqrt(gamma * max(max(tem))); umax = vel_max + c_max; dt = cfl_hx / umax; dt0 = 0.24 * h / max(ha, (mu + 0.5 * vel_max)); dt = min(dt, dt0); end %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ function B = apply_boundary_condition(A) [nx, ny] = size(A); B = zeros(nx + 2, ny + 2); B(2:end - 1, 2:end - 1) = A; B = periodic_2d(B); end %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ function rnd_data = get_rnd_data(data_flag) if (data_flag == 1) fin = strcat('./rndata/rnd1.txt'); else fin = strcat('./rndata/rnd2.txt'); end rnd_data = dlmread(fin); end %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ function [line_low, line_high] = getP_position_KH(nx, h, epsilon) rndl = get_rnd_data(1); rndh = get_rnd_data(2); nlength = nx + 2; line_low = zeros(nlength, 1); line_high = zeros(nlength, 1); K = 10; for i = 1:nlength xi = (i - 0.5) * h; PY_low = 0; PY_high = 0; for n = 1:K PY_low = PY_low + rndl(n) * cos(rndl(n + 10) + 2 * n * pi * xi); PY_high = PY_high + rndh(n) * cos(rndh(n + 10) + 2 * n * pi * xi); end line_low(i) = 0.25 + epsilon * PY_low; line_high(i) = 0.75 + epsilon * PY_high; end end %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ function [line_low, line_high] = getP_position_KH_rand(nx, h, epsilon) K = 10; randla = rand(K, 1); randlb = -pi + (pi + pi) * rand(K, 1); randha = rand(K, 1); randhb = -pi + (pi + pi) * rand(K, 1); suml = sum(randla); randla=randla ./ suml; sumh = sum(randha); randha=randha ./ sumh; nlength = nx + 2; line_low = zeros(nlength, 1);