Mercurial > repos > public > sbplib
changeset 570:7d5d74940987 feature/grids
Merge laplace refactor. Refactor class to more closely resemble SBP notation. Fix bc bugs
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Fri, 01 Sep 2017 11:07:15 +0200 |
| parents | 12ee11893453 (current diff) f1a01a48779c (diff) |
| children | f8072bb8d1d1 |
| files | |
| diffstat | 2 files changed, 133 insertions(+), 142 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/LaplaceCurvilinear.m Fri Sep 01 10:52:56 2017 +0200 +++ b/+scheme/LaplaceCurvilinear.m Fri Sep 01 11:07:15 2017 +0200 @@ -7,32 +7,37 @@ order % Order accuracy for the approximation - D % non-stabalized scheme operator - M % Derivative norm - a,b + a,b % Parameters of the operator + + + % Inner products and operators for physical coordinates + D % Laplace operator + H, Hi % Inner product + e_w, e_e, e_s, e_n + d_w, d_e, d_s, d_n % Normal derivatives at the boundary + H_w, H_e, H_s, H_n % Boundary inner products + Dx, Dy % Physical derivatives + M % Gradient inner product + + % Metric coefficients J, Ji a11, a12, a22 + x_u + x_v + y_u + y_v - H % Discrete norm - Hi + % Inner product and operators for logical coordinates H_u, H_v % Norms in the x and y directions + Hi_u, Hi_v Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. - Hi_u, Hi_v Hiu, Hiv - e_w, e_e, e_s, e_n du_w, dv_w du_e, dv_e du_s, dv_s du_n, dv_n gamm_u, gamm_v lambda - - Dx, Dy % Physical derivatives - - x_u - x_v - y_u - y_v end methods @@ -44,7 +49,6 @@ default_arg('a', 1); default_arg('b', 1); - if b ~=1 error('Not implemented yet') end @@ -60,7 +64,8 @@ h_u = h(1); h_v = h(2); - % Operators + + % 1D operators ops_u = opSet(m_u, {0, 1}, order); ops_v = opSet(m_v, {0, 1}, order); @@ -85,10 +90,30 @@ d1_l_v = ops_v.d1_l; d1_r_v = ops_v.d1_r; + + % Logical operators Du = kr(D1_u,I_v); Dv = kr(I_u,D1_v); + obj.Hu = kr(H_u,I_v); + obj.Hv = kr(I_u,H_v); + obj.Hiu = kr(Hi_u,I_v); + obj.Hiv = kr(I_u,Hi_v); - % Metric derivatives + e_w = kr(e_l_u,I_v); + e_e = kr(e_r_u,I_v); + e_s = kr(I_u,e_l_v); + e_n = kr(I_u,e_r_v); + obj.du_w = kr(d1_l_u,I_v); + obj.dv_w = (e_w'*Dv)'; + obj.du_e = kr(d1_r_u,I_v); + obj.dv_e = (e_e'*Dv)'; + obj.du_s = (e_s'*Du)'; + obj.dv_s = kr(I_u,d1_l_v); + obj.du_n = (e_n'*Du)'; + obj.dv_n = kr(I_u,d1_r_v); + + + % Metric coefficients coords = g.points(); x = coords(:,1); y = coords(:,2); @@ -104,8 +129,14 @@ a22 = 1./J .* (x_u.^2 + y_u.^2); lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); + obj.x_u = x_u; + obj.x_v = x_v; + obj.y_u = y_u; + obj.y_v = y_v; + + % Assemble full operators - L_12 = spdiags(a12, 0, m_tot, m_tot); + L_12 = spdiag(a12); Duv = Du*L_12*Dv; Dvu = Dv*L_12*Du; @@ -125,31 +156,55 @@ Dvv(p,p) = D; end - obj.H = kr(H_u,H_v); - obj.Hi = kr(Hi_u,Hi_v); - obj.Hu = kr(H_u,I_v); - obj.Hv = kr(I_u,H_v); - obj.Hiu = kr(Hi_u,I_v); - obj.Hiv = kr(I_u,Hi_v); + + % Physical operators + obj.J = spdiag(J); + obj.Ji = spdiag(1./J); + + obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv); + obj.H = obj.J*kr(H_u,H_v); + obj.Hi = obj.Ji*kr(Hi_u,Hi_v); + + obj.e_w = e_w; + obj.e_e = e_e; + obj.e_s = e_s; + obj.e_n = e_n; + + %% normal derivatives + I_w = ind(1,:); + I_e = ind(end,:); + I_s = ind(:,1); + I_n = ind(:,end); - obj.e_w = kr(e_l_u,I_v); - obj.e_e = kr(e_r_u,I_v); - obj.e_s = kr(I_u,e_l_v); - obj.e_n = kr(I_u,e_r_v); - obj.du_w = kr(d1_l_u,I_v); - obj.dv_w = (obj.e_w'*Dv)'; - obj.du_e = kr(d1_r_u,I_v); - obj.dv_e = (obj.e_e'*Dv)'; - obj.du_s = (obj.e_s'*Du)'; - obj.dv_s = kr(I_u,d1_l_v); - obj.du_n = (obj.e_n'*Du)'; - obj.dv_n = kr(I_u,d1_r_v); + a11_w = spdiag(a11(I_w)); + a12_w = spdiag(a12(I_w)); + a11_e = spdiag(a11(I_e)); + a12_e = spdiag(a12(I_e)); + a22_s = spdiag(a22(I_s)); + a12_s = spdiag(a12(I_s)); + a22_n = spdiag(a22(I_n)); + a12_n = spdiag(a12(I_n)); + + s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2); + s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2); + s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2); + s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2); - obj.x_u = x_u; - obj.x_v = x_v; - obj.y_u = y_u; - obj.y_v = y_v; + obj.d_w = -1*(spdiag(1./s_w)*(a11_w*obj.du_w' + a12_w*obj.dv_w'))'; + obj.d_e = (spdiag(1./s_e)*(a11_e*obj.du_e' + a12_e*obj.dv_e'))'; + obj.d_s = -1*(spdiag(1./s_s)*(a22_s*obj.dv_s' + a12_s*obj.du_s'))'; + obj.d_n = (spdiag(1./s_n)*(a22_n*obj.dv_n' + a12_n*obj.du_n'))'; + + obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; + obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; + %% Boundary inner products + obj.H_w = H_v*spdiag(s_w); + obj.H_e = H_v*spdiag(s_e); + obj.H_s = H_u*spdiag(s_s); + obj.H_n = H_u*spdiag(s_n); + + % Misc. obj.m = m; obj.h = [h_u h_v]; obj.order = order; @@ -157,17 +212,11 @@ obj.a = a; obj.b = b; - obj.J = spdiags(J, 0, m_tot, m_tot); - obj.Ji = spdiags(1./J, 0, m_tot, m_tot); obj.a11 = a11; obj.a12 = a12; obj.a22 = a22; - obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv); obj.lambda = lambda; - obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; - obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; - obj.gamm_u = h_u*ops_u.borrowing.M.d1; obj.gamm_v = h_v*ops_v.borrowing.M.d1; end @@ -184,62 +233,34 @@ default_arg('type','neumann'); default_arg('parameter', []); - [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv , ~, ~, ~, scale_factor] = obj.get_boundary_ops(boundary); + [e, d, gamm, H_b, ~] = obj.get_boundary_ops(boundary); switch type % Dirichlet boundary condition case {'D','d','dirichlet'} - % 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; + b1 = gamm*obj.lambda./obj.a11.^2; + b2 = gamm*obj.lambda./obj.a22.^2; - b1 = gamm*u.lambda./u.a11.^2; - b2 = gamm*u.lambda./u.a22.^2; + tau1 = tuning * spdiag(-1./b1 - 1./b2); + tau2 = 1; - tau = -1./b1 - 1./b2; - tau = tuning * spdiag(tau); - sig1 = 1; + tau = (tau1*e + tau2*d)*H_b; - penalty_parameter_1 = halfnorm_inv_n*(tau + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e; - - closure = obj.Ji*obj.a * penalty_parameter_1*e'; - penalty = -obj.Ji*obj.a * penalty_parameter_1; + closure = obj.a*obj.Hi*tau*e'; + penalty = -obj.a*obj.Hi*tau; % Neumann boundary condition case {'N','n','neumann'} - a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); - a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); - d = (a_n * d_n' + a_t*d_t')'; - - tau1 = -s; + tau1 = -1; tau2 = 0; - tau = obj.a * obj.Ji*(tau1*e + tau2*d); - - closure = halfnorm_inv*tau*d'; - penalty = -halfnorm_inv*tau; + tau = (tau1*e + tau2*d)*H_b; - % Characteristic boundary condition - case {'characteristic', 'char', 'c'} - default_arg('parameter', 1); - beta = parameter; + closure = obj.a*obj.Hi*tau*d'; + penalty = -obj.a*obj.Hi*tau; - a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); - a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); - d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative - - tau = -obj.a * 1/beta*obj.Ji*e; - - closure{1} = halfnorm_inv*tau*spdiag(scale_factor)*e'; - closure{2} = halfnorm_inv*tau*beta*d'; - penalty = -halfnorm_inv*tau; % Unknown, boundary condition otherwise @@ -252,16 +273,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, 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); - a_n_v = spdiag(coeff_n_v); - a_t_v = spdiag(coeff_t_v); - - F_u = (s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u')'; - F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')'; + [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); + [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); u = obj; v = neighbour_scheme; @@ -271,24 +284,25 @@ 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); - tau = tuning * spdiag(tau); - sig1 = 1/2; - sig2 = -1/2; + tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); + tau1 = tuning * spdiag(tau1); + tau2 = 1/2; - 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; + sig1 = -1/2; + sig2 = 0; + tau = (e_u*tau1 + tau2*d_u)*H_b_u; + sig = (sig1*e_u + sig2*d_u)*H_b_u; - closure = obj.Ji*obj.a * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u'); - penalty = obj.Ji*obj.a * (-penalty_parameter_1*e_v' + penalty_parameter_2*F_v'); + closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u'); + penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v'); end % Ruturns the boundary ops and sign for the boundary specified by the string boundary. % The right boundary is considered the positive 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, scale_factor] = get_boundary_ops(obj, boundary) + function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary) % gridMatrix = zeros(obj.m(2),obj.m(1)); % gridMatrix(:) = 1:numel(gridMatrix); @@ -298,58 +312,32 @@ switch boundary case 'w' e = obj.e_w; - d_n = obj.du_w; - d_t = obj.dv_w; - s = -1; - + d = obj.d_w; + H_b = obj.H_w; I = ind(1,:); - coeff_n = obj.a11(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); case 'e' e = obj.e_e; - d_n = obj.du_e; - d_t = obj.dv_e; - s = 1; - + d = obj.d_e; + H_b = obj.H_e; I = ind(end,:); - coeff_n = obj.a11(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); case 's' e = obj.e_s; - d_n = obj.dv_s; - d_t = obj.du_s; - s = -1; - + d = obj.d_s; + H_b = obj.H_s; I = ind(:,1)'; - coeff_n = obj.a22(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); case 'n' e = obj.e_n; - d_n = obj.dv_n; - d_t = obj.du_n; - s = 1; - + d = obj.d_n; + H_b = obj.H_n; I = ind(:,end)'; - coeff_n = obj.a22(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); otherwise error('No such boundary: boundary = %s',boundary); end switch boundary case {'w','e'} - halfnorm_inv_n = obj.Hiu; - halfnorm_inv_t = obj.Hiv; - halfnorm_t = obj.Hv; gamm = obj.gamm_u; case {'s','n'} - halfnorm_inv_n = obj.Hiv; - halfnorm_inv_t = obj.Hiu; - halfnorm_t = obj.Hu; gamm = obj.gamm_v; end end @@ -357,7 +345,5 @@ function N = size(obj) N = prod(obj.m); end - - end -end \ No newline at end of file +end