Mercurial > repos > public > sbplib
changeset 97:33dba20b4b9d feature/arclen-param
Merged with deafult.
| author | Martin Almquist <martin.almquist@it.uu.se> |
|---|---|
| date | Wed, 02 Dec 2015 16:29:54 +0100 |
| parents | 7b092f0fd620 (current diff) 19d0c9325a3e (diff) |
| children | b30f3d8845f4 |
| files | +multiblock/equal_step_size.m |
| diffstat | 8 files changed, 141 insertions(+), 51 deletions(-) [+] |
line wrap: on
line diff
diff -r 7b092f0fd620 -r 33dba20b4b9d +anim/make_movie.sh --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+anim/make_movie.sh Wed Dec 02 16:29:54 2015 +0100 @@ -0,0 +1,15 @@ +#!/bin/bash + +# $1 is the folder where the PNGS are found + +#ffmpeg -qp 0 -framerate 30 -i test%04d.png -c:v mpeg4 -r 30 $1_pos.mp4 +# ffmpeg -framerate 30 -i $1/test_pos1%04d.png -q:v 0 -c:v mpeg4 -r 30 $1_pos1.mp4 +# ffmpeg -framerate 30 -i $1/test_pos2%04d.png -q:v 0 -c:v mpeg4 -r 30 $1_pos2.mp4 +ffmpeg -framerate 30 -i $1/fig%04d.png -q:v 0 -c:v mpeg4 -r 30 $1.mp4 +# rm -f test*.png + + + +# -b 5000k to set bitrate of file. +# -q:v 0 to set quality? This is the shit!!! +# -crf set quality 23 is default 18 is almost lossless 51 is worst \ No newline at end of file
diff -r 7b092f0fd620 -r 33dba20b4b9d +anim/setup_fig_mov.m --- a/+anim/setup_fig_mov.m Mon Nov 30 12:13:22 2015 +0100 +++ b/+anim/setup_fig_mov.m Wed Dec 02 16:29:54 2015 +0100 @@ -2,6 +2,11 @@ fresh_dir(dirname) + figure_handle.Units = 'centimeters'; + figure_handle.PaperUnits = 'centimeters'; + figure_handle.PaperPosition(3:4) = figure_handle.Position(3:4); + + frame_nr = 0; function save_frame_F() saveas(figure_handle,sprintf('%s/%s%04d.png',dirname,'fig',frame_nr));
diff -r 7b092f0fd620 -r 33dba20b4b9d +anim/setup_time_quantity_plot.m --- a/+anim/setup_time_quantity_plot.m Mon Nov 30 12:13:22 2015 +0100 +++ b/+anim/setup_time_quantity_plot.m Wed Dec 02 16:29:54 2015 +0100 @@ -11,8 +11,6 @@ end axis_handle = gca; - legend() - function update(t_now,varargin) if ishandle(axis_handle)
diff -r 7b092f0fd620 -r 33dba20b4b9d +multiblock/equal_step_size.m --- a/+multiblock/equal_step_size.m Mon Nov 30 12:13:22 2015 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -% Choose the number of points on m, M, so that the step size is equal to that for N points on n. -function [M] = equal_step_size(n,N,m) - M = round((m*(N-1)+n)/n); -end \ No newline at end of file
diff -r 7b092f0fd620 -r 33dba20b4b9d +scheme/Euler1d.m --- a/+scheme/Euler1d.m Mon Nov 30 12:13:22 2015 +0100 +++ b/+scheme/Euler1d.m Wed Dec 02 16:29:54 2015 +0100 @@ -17,6 +17,14 @@ end + properties (Constant) + SUBSONIC_INFLOW = 1; + SUBSONIC_OUTFLOW = -1; + NO_FLOW = 0; + SUPERSONIC_INFLOW = 2; + SUPERSONIC_OUTFLOW = -2; + end + methods function obj = Euler1d(m,xlim,order,gama,opsGen,do_upwind) default_arg('opsGen',@sbp.Ordinary); @@ -155,6 +163,40 @@ % Devide columns by stuff to get rid of extra comp? end + function fs = flowStateL(obj, q) + q_l = obj.e_L'*q; + c = obj.c(q); + v = q_l(2,:)/q_l(1,:); + + if v > c + fs = scheme.Euler1d.SUPERSONIC_INFLOW; + elseif v > 0 + fs = scheme.Euler1d.SUBSONIC_INFLOW; + elseif v > -c + fs = scheme.Euler1d.SUBSONIC_OUTFLOW; + else + fs = scheme.Euler1d.SUPERSONIC_INFLOW; + end + end + + % returns positiv values for inlfow, negative for outflow. + % +-1 for subsonic + function fs = flowStateR(obj, q) + q_r = obj.e_R'*q; + c = obj.c(q); + v = q_r(2,:)/q_r(1,:); + + if v < -c + fs = scheme.Euler1d.SUPERSONIC_INFLOW; + elseif v < 0 + fs = scheme.Euler1d.SUBSONIC_INFLOW; + elseif v < c + fs = scheme.Euler1d.SUBSONIC_OUTFLOW; + else + fs = scheme.Euler1d.SUPERSONIC_INFLOW; + end + end + % Enforces the boundary conditions % w+ = R*w- + g(t) function closure = boundary_condition(obj,boundary, type, varargin) @@ -170,7 +212,9 @@ closure = obj.boundary_condition_inflow(boundary,varargin{:}); case 'outflow' closure = obj.boundary_condition_outflow(boundary,varargin{:}); - case 'outflow_rho' + case 'inflow_rho' + closure = obj.boundary_condition_inflow_rho(boundary,varargin{:}); + case 'outflow_rho' closure = obj.boundary_condition_outflow_rho(boundary,varargin{:}); case 'char' closure = obj.boundary_condition_char(boundary,varargin{:}); @@ -207,9 +251,7 @@ Tin = T(:,p_in); Tot = T(:,p_ot); - % Convert bc from ordinary form to characteristic form. - % Lq = g => w_in = Rw_ot + g_tilde - + % Calculate eigen value matrix Lambda = obj.Lambda(q_s); % Setup the penalty parameter @@ -217,18 +259,11 @@ tau2 = zeros(length(p_ot),length(p_in)); % Penalty only on ingoing char. tauHat = [tau1; tau2]; - tau = -s*e_S*sparse(T*tauHat(pt,:)); + tau = e_S*sparse(T*tauHat(pt,:)); L = L_fun(rho,u,e); g = g_fun(t); - % printExpr('s') - % penalty = tauHat(pt,:)*inv(L*Tin)*(L*q_s - g); - % tauHatPt = tauHat(pt,:); - % display(tauHatPt); - % display(penalty); - % pause - o = 1/2*obj.Hi * tau * inv(L*Tin)*(L*q_s - g); end closure = @closure_fun; @@ -317,6 +352,29 @@ end + function closure = boundary_condition_inflow_rho(obj, boundary, rho_data, v_data) + [~,~,s] = obj.get_boundary_ops(boundary); + + switch s + case -1 + p_in = [1 2]; + case 1 + p_in = [1 3]; + end + + a = obj.gamma - 1; + L = @(rho,~,~)[ + 0 1/rho 0; + 1 0 0; + ]; + g = @(t)[ + v_data(t); + rho_data(t); + ]; + + closure = boundary_condition_L(obj, boundary, L, g, p_in); + end + function closure = boundary_condition_outflow_rho(obj, boundary, rho_data) [~,~,s] = obj.get_boundary_ops(boundary);
diff -r 7b092f0fd620 -r 33dba20b4b9d +scheme/Wave2dCurve.m --- a/+scheme/Wave2dCurve.m Mon Nov 30 12:13:22 2015 +0100 +++ b/+scheme/Wave2dCurve.m Wed Dec 02 16:29:54 2015 +0100 @@ -175,19 +175,29 @@ switch type % Dirichlet boundary condition case {'D','d','dirichlet'} - error('not implemented') - alpha = obj.alpha; + % v denotes the solution in the neighbour domain + tuning = 1.2; + % tuning = 20.2; + [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary); + + a_n = spdiag(coeff_n); + a_t = spdiag(coeff_t); + + F = (s * a_n * d_n' + s * a_t*d_t')'; + + u = obj; - % tau1 < -alpha^2/gamma - tuning = 1.1; - tau1 = -tuning*alpha/gamm; - tau2 = s*alpha; + b1 = gamm*u.lambda./u.a11.^2; + b2 = gamm*u.lambda./u.a22.^2; - p = tau1*e + tau2*d; + tau = -1./b1 - 1./b2; + tau = tuning * spdiag(tau(:)); + sig1 = 1/2; - closure = halfnorm_inv*p*e'; + penalty_parameter_1 = halfnorm_inv_n*(tau + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e; - pp = halfnorm_inv*p; + closure = obj.Ji*obj.c^2 * penalty_parameter_1*e'; + pp = -obj.Ji*obj.c^2 * penalty_parameter_1; switch class(data) case 'double' penalty = pp*data; @@ -234,8 +244,8 @@ % v denotes the solution in the neighbour domain tuning = 1.2; % tuning = 20.2; - [e_u, d_n_u, d_t_u, coeff_n_u, coeff_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t] = obj.get_boundary_ops(boundary); - [e_v, d_n_v, d_t_v, coeff_n_v, coeff_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t] = neighbour_scheme.get_boundary_ops(boundary); + [e_u, d_n_u, d_t_u, coeff_n_u, coeff_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t, I_u] = obj.get_boundary_ops(boundary); + [e_v, d_n_v, d_t_v, coeff_n_v, coeff_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); a_n_u = spdiag(coeff_n_u); a_t_u = spdiag(coeff_t_u); @@ -248,29 +258,33 @@ u = obj; v = neighbour_scheme; - b1_u = gamm_u*u.lambda./u.a11.^2; - b2_u = gamm_u*u.lambda./u.a22.^2; - b1_v = gamm_v*v.lambda./v.a11.^2; - b2_v = gamm_v*v.lambda./v.a22.^2; + b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; + b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; + b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; + b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; - - tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); - m_tot = obj.m(1)*obj.m(2); + tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); tau = tuning * spdiag(tau(:)); sig1 = 1/2; - sig2 = -1/2*s_u; + sig2 = -1/2; - penalty_parameter_1 = halfnorm_inv_u_n*(tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t)*e_u; + penalty_parameter_1 = halfnorm_inv_u_n*(e_u*tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u); penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u; closure = obj.Ji*obj.c^2 * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u'); - penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' - penalty_parameter_2*F_v'); + penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' + penalty_parameter_2*F_v'); end % Ruturns the boundary ops and sign for the boundary specified by the string boundary. % The right boundary is considered the positive boundary - function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,a_n, a_t] = get_boundary_ops(obj,boundary) + % + % I -- the indecies of the boundary points in the grid matrix + function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I] = get_boundary_ops(obj,boundary) + + gridMatrix = zeros(obj.m(2),obj.m(1)); + gridMatrix(:) = 1:numel(gridMatrix); + switch boundary case 'w' e = obj.e_w; @@ -278,32 +292,36 @@ d_t = obj.dv_w; s = -1; - coeff_n = obj.a11(:,1); - coeff_t = obj.a12(:,1); + I = gridMatrix(:,1); + coeff_n = obj.a11(I); + coeff_t = obj.a12(I); case 'e' e = obj.e_e; d_n = obj.du_e; d_t = obj.dv_e; s = 1; - coeff_n = obj.a11(:,end); - coeff_t = obj.a12(:,end); + I = gridMatrix(:,end); + coeff_n = obj.a11(I); + coeff_t = obj.a12(I); case 's' e = obj.e_s; d_n = obj.dv_s; d_t = obj.du_s; s = -1; - coeff_n = obj.a22(1,:)'; - coeff_t = obj.a12(1,:)'; + I = gridMatrix(1,:)'; + coeff_n = obj.a22(I); + coeff_t = obj.a12(I); case 'n' e = obj.e_n; d_n = obj.dv_n; d_t = obj.du_n; s = 1; - coeff_n = obj.a22(end,:)'; - coeff_t = obj.a12(end,:)'; + I = gridMatrix(end,:)'; + coeff_n = obj.a22(I); + coeff_t = obj.a12(I); otherwise error('No such boundary: boundary = %s',boundary); end
diff -r 7b092f0fd620 -r 33dba20b4b9d printSize.m --- a/printSize.m Mon Nov 30 12:13:22 2015 +0100 +++ b/printSize.m Wed Dec 02 16:29:54 2015 +0100 @@ -1,5 +1,4 @@ function printSize(A) - warning('Deprecated! Use printExpr() instead!'); - s = size(A); - fprintf('%8s has size: [%d, %d]\n',inputname(1),s(1),s(2)); + result = size(A); + fprintf('size(%s) => %s\n', inputname(1), toString(result)); end \ No newline at end of file