Mercurial > repos > public > sbplib
changeset 905:459eeb99130f feature/utux2D
Include type as (optional) input parameter in the interface method of all schemes.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Thu, 22 Nov 2018 22:03:44 -0800 |
| parents | 14b093a344eb |
| children | 0499239496cf |
| files | +scheme/Beam.m +scheme/Beam2d.m +scheme/Elastic2dVariable.m +scheme/Euler1d.m +scheme/Heat2dVariable.m +scheme/Hypsyst2d.m +scheme/Hypsyst2dCurve.m +scheme/Hypsyst3d.m +scheme/Hypsyst3dCurve.m +scheme/Scheme.m +scheme/Schrodinger.m +scheme/Schrodinger2d.m +scheme/Utux.m +scheme/Utux2D.m +scheme/Wave.m +scheme/Wave2d.m +scheme/Wave2dCurve.m |
| diffstat | 17 files changed, 325 insertions(+), 320 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Beam.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Beam.m Thu Nov 22 22:03:44 2018 -0800 @@ -170,7 +170,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain [e_u,d1_u,d2_u,d3_u,s_u] = obj.get_boundary_ops(boundary);
--- a/+scheme/Beam2d.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Beam2d.m Thu Nov 22 22:03:44 2018 -0800 @@ -161,7 +161,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain [e_u,d1_u,d2_u,d3_u,s_u,gamm_u,delt_u, halfnorm_inv] = obj.get_boundary_ops(boundary);
--- a/+scheme/Elastic2dVariable.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Elastic2dVariable.m Thu Nov 22 22:03:44 2018 -0800 @@ -1,7 +1,7 @@ classdef Elastic2dVariable < scheme.Scheme % Discretizes the elastic wave equation: -% rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i +% rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i % opSet should be cell array of opSets, one per dimension. This % is useful if we have periodic BC in one direction. @@ -37,7 +37,7 @@ e_l, e_r d1_l, d1_r % Normal derivatives at the boundary E % E{i}^T picks out component i - + H_boundary % Boundary inner products % Kroneckered norms and coefficients @@ -223,14 +223,14 @@ tau_l{j}{i} = sparse(m_tot,dim*m_tot); tau_r{j}{i} = sparse(m_tot,dim*m_tot); for k = 1:dim - T_l{j}{i,k} = ... + T_l{j}{i,k} = ... -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})... - -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})... + -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})... -d(i,k)*MU*e_l{j}*d1_l{j}'; - T_r{j}{i,k} = ... + T_r{j}{i,k} = ... d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})... - +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})... + +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})... +d(i,k)*MU*e_r{j}*d1_r{j}'; tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; @@ -271,7 +271,7 @@ default_arg('parameter', []); % j is the coordinate direction of the boundary - % nj: outward unit normal component. + % nj: outward unit normal component. % nj = -1 for west, south, bottom boundaries % nj = 1 for east, north, top boundaries [j, nj] = obj.get_boundary_number(boundary); @@ -329,20 +329,20 @@ db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... + d(i,j)* a_mu_i*MU ... - + db(i,j)*a_mu_ij*MU ); + + db(i,j)*a_mu_ij*MU ); % Loop over components that Dirichlet penalties end up on for i = 1:dim C = T{k,i}; A = -d(i,k)*alpha(i,j); B = A + C; - closure = closure + E{i}*RHOi*Hi*B'*e{j}*H_gamma*(e{j}'*E{k}' ); + closure = closure + E{i}*RHOi*Hi*B'*e{j}*H_gamma*(e{j}'*E{k}' ); penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e{j}*H_gamma; - end + end % Free boundary condition case {'F','f','Free','free','traction','Traction','t','T'} - closure = closure - E{k}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{k} ); + closure = closure - E{k}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{k} ); penalty{k} = penalty{k} + E{k}*RHOi*Hi*e{j}*H_gamma; % Unknown boundary condition @@ -352,7 +352,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain tuning = 1.2;
--- a/+scheme/Euler1d.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Euler1d.m Thu Nov 22 22:03:44 2018 -0800 @@ -446,7 +446,7 @@ closure = @closure_fun; end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) error('NOT DONE') % u denotes the solution in the own domain % v denotes the solution in the neighbour domain
--- a/+scheme/Heat2dVariable.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Heat2dVariable.m Thu Nov 22 22:03:44 2018 -0800 @@ -1,9 +1,9 @@ classdef Heat2dVariable < scheme.Scheme % Discretizes the Laplacian with variable coefficent, -% In the Heat equation way (i.e., the discretization matrix is not necessarily +% In the Heat equation way (i.e., the discretization matrix is not necessarily % symmetric) -% u_t = div * (kappa * grad u ) +% u_t = div * (kappa * grad u ) % opSet should be cell array of opSets, one per dimension. This % is useful if we have periodic BC in one direction. @@ -28,7 +28,7 @@ H, Hi % Inner products e_l, e_r d1_l, d1_r % Normal derivatives at the boundary - + H_boundary % Boundary inner products end @@ -160,7 +160,7 @@ default_arg('parameter', []); % j is the coordinate direction of the boundary - % nj: outward unit normal component. + % nj: outward unit normal component. % nj = -1 for west, south, bottom boundaries % nj = 1 for east, north, top boundaries [j, nj] = obj.get_boundary_number(boundary); @@ -176,19 +176,19 @@ Hi = obj.Hi; H_gamma = obj.H_boundary{j}; KAPPA = obj.KAPPA; - kappa_gamma = e{j}'*KAPPA*e{j}; + kappa_gamma = e{j}'*KAPPA*e{j}; switch type % Dirichlet boundary condition case {'D','d','dirichlet','Dirichlet'} - closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' ); + closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' ); penalty = nj*Hi*d{j}*kappa_gamma*H_gamma; % Free boundary condition case {'N','n','neumann','Neumann'} - closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' ); - penalty = nj*Hi*e{j}*kappa_gamma*H_gamma; + closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' ); + penalty = nj*Hi*e{j}*kappa_gamma*H_gamma; % Unknown boundary condition otherwise @@ -196,7 +196,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain error('Interface not implemented');
--- a/+scheme/Hypsyst2d.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Hypsyst2d.m Thu Nov 22 22:03:44 2018 -0800 @@ -6,10 +6,10 @@ x,y % Grid X,Y % Values of x and y for each grid point order % Order accuracy for the approximation - + D % non-stabalized scheme operator A, B, E %Coefficient matrices - + H % Discrete norm % Norms in the x and y directions Hxi,Hyi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. @@ -17,65 +17,65 @@ e_w, e_e, e_s, e_n params %parameters for the coeficient matrice end - + methods %Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu function obj = Hypsyst2d(m, lim, order, A, B, E, params) default_arg('E', []) xlim = lim{1}; ylim = lim{2}; - + if length(m) == 1 m = [m m]; end - + obj.A=A; obj.B=B; obj.E=E; - + m_x = m(1); m_y = m(2); obj.params = params; - + ops_x = sbp.D2Standard(m_x,xlim,order); ops_y = sbp.D2Standard(m_y,ylim,order); - + obj.x = ops_x.x; obj.y = ops_y.x; - + obj.X = kr(obj.x,ones(m_y,1)); obj.Y = kr(ones(m_x,1),obj.y); - + Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y); Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y); Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y); - + obj.n = length(A(obj.params,0,0)); - + I_n = eye(obj.n);I_x = speye(m_x); obj.I_x = I_x; I_y = speye(m_y); obj.I_y = I_y; - - + + D1_x = kr(I_n, ops_x.D1, I_y); obj.Hxi = kr(I_n, ops_x.HI, I_y); D1_y = kr(I_n, I_x, ops_y.D1); obj.Hyi = kr(I_n, I_x, ops_y.HI); - + obj.e_w = kr(I_n, ops_x.e_l, I_y); obj.e_e = kr(I_n, ops_x.e_r, I_y); obj.e_s = kr(I_n, I_x, ops_y.e_l); obj.e_n = kr(I_n, I_x, ops_y.e_r); - + obj.m = m; obj.h = [ops_x.h ops_y.h]; obj.order = order; - + obj.D = -Aevaluated*D1_x-Bevaluated*D1_y-Eevaluated; - + end - + % Closure functions return the opertors applied to the own doamin to close the boundary % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. @@ -92,18 +92,18 @@ error('No such boundary condition') end end - - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) error('An interface function does not exist yet'); end - + function N = size(obj) N = obj.m; end - + function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y) params = obj.params; - + if isa(mat,'function_handle') [rows,cols] = size(mat(params,0,0)); matVec = mat(params,X',Y'); @@ -116,7 +116,7 @@ cols = cols/side; end ret = cell(rows,cols); - + for ii = 1:rows for jj=1:cols ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); @@ -124,13 +124,13 @@ end ret = cell2mat(ret); end - + %Characteristic boundary conditions function [closure, penalty] = boundary_condition_char(obj,boundary) params = obj.params; x = obj.x; y = obj.y; - + switch boundary case {'w','W','west'} e_ = obj.e_w; @@ -164,7 +164,7 @@ pos = signVec(1); zeroval = signVec(2); neg = signVec(3); - + switch boundPos case {'l'} tau = sparse(obj.n*side,pos); @@ -180,13 +180,13 @@ penalty = -Hi*e_*V*tau*Vi_minus; end end - + % General boundary condition in the form Lu=g(x) function [closure,penalty] = boundary_condition_general(obj,boundary,L) params = obj.params; x = obj.x; y = obj.y; - + switch boundary case {'w','W','west'} e_ = obj.e_w; @@ -218,14 +218,14 @@ boundPos = 'r'; Hi = obj.Hyi; [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end)); - L = obj.evaluateCoefficientMatrix(L,x,y(end)); + L = obj.evaluateCoefficientMatrix(L,x,y(end)); side = max(length(x)); end - + pos = signVec(1); zeroval = signVec(2); neg = signVec(3); - + switch boundPos case {'l'} tau = sparse(obj.n*side,pos); @@ -233,7 +233,7 @@ Vi_minus = Vi(pos+zeroval+1:obj.n*side,:); V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); - + tau(1:pos,:) = -abs(D(1:pos,1:pos)); R = -inv(L*V_plus)*(L*V_minus); closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; @@ -243,7 +243,7 @@ tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); Vi_plus = Vi(1:pos,:); Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); - + V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); R = -inv(L*V_minus)*(L*V_plus); @@ -251,13 +251,13 @@ penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; end end - + % Function that diagonalizes a symbolic matrix A as A=V*D*Vi % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign % [d+ ] % D = [ d0 ] % [ d-] - % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D + % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y) params = obj.params; syms xs ys @@ -265,12 +265,12 @@ Vi = inv(V); xs = x; ys = y; - + side = max(length(x),length(y)); Dret = zeros(obj.n,side*obj.n); Vret = zeros(obj.n,side*obj.n); Viret = zeros(obj.n,side*obj.n); - + for ii = 1:obj.n for jj = 1:obj.n Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); @@ -278,7 +278,7 @@ Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); end end - + D = sparse(Dret); V = sparse(Vret); Vi = sparse(Viret); @@ -286,16 +286,16 @@ Vi = obj.evaluateCoefficientMatrix(Vi,x,y); D = obj.evaluateCoefficientMatrix(D,x,y); DD = diag(D); - + poseig = (DD>0); zeroeig = (DD==0); negeig = (DD<0); - + D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; end - + end end \ No newline at end of file
--- a/+scheme/Hypsyst2dCurve.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Hypsyst2dCurve.m Thu Nov 22 22:03:44 2018 -0800 @@ -4,19 +4,19 @@ n % size of system h % Grid spacing X,Y % Values of x and y for each grid point - + J, Ji % Jacobaian and inverse Jacobian xi,eta Xi,Eta - + A,B X_eta, Y_eta X_xi,Y_xi order % Order accuracy for the approximation - + D % non-stabalized scheme operator Ahat, Bhat, E - + H % Discrete norm Hxii,Hetai % Kroneckerd norms in xi and eta. I_xi,I_eta, I_N, onesN @@ -24,93 +24,93 @@ index_w, index_e,index_s,index_n params % Parameters for the coeficient matrice end - - + + methods % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu function obj = Hypsyst2dCurve(m, order, A, B, E, params, ti) default_arg('E', []) xilim = {0 1}; etalim = {0 1}; - + if length(m) == 1 m = [m m]; end obj.params = params; obj.A=A; obj.B=B; - + obj.Ahat=@(params,x,y,x_eta,y_eta)(A(params,x,y).*y_eta-B(params,x,y).*x_eta); obj.Bhat=@(params,x,y,x_xi,y_xi)(B(params,x,y).*x_xi-A(params,x,y).*y_xi); obj.E=@(params,x,y,~,~)E(params,x,y); - + m_xi = m(1); m_eta = m(2); m_tot=m_xi*m_eta; - + ops_xi = sbp.D2Standard(m_xi,xilim,order); ops_eta = sbp.D2Standard(m_eta,etalim,order); - + obj.xi = ops_xi.x; obj.eta = ops_eta.x; - + obj.Xi = kr(obj.xi,ones(m_eta,1)); obj.Eta = kr(ones(m_xi,1),obj.eta); - + obj.n = length(A(obj.params,0,0)); obj.onesN=ones(obj.n); - + obj.index_w=1:m_eta; - obj.index_e=(m_tot-m_e - + obj.index_e=(m_tot-m_e + metric_termsta+1):m_tot; obj.index_s=1:m_eta:(m_tot-m_eta+1); obj.index_n=(m_eta):m_eta:m_tot; - + I_n = eye(obj.n); I_xi = speye(m_xi); obj.I_xi = I_xi; I_eta = speye(m_eta); obj.I_eta = I_eta; - + D1_xi = kr(I_n, ops_xi.D1, I_eta); obj.Hxii = kr(I_n, ops_xi.HI, I_eta); D1_eta = kr(I_n, I_xi, ops_eta.D1); obj.Hetai = kr(I_n, I_xi, ops_eta.HI); - + obj.e_w = kr(I_n, ops_xi.e_l, I_eta); obj.e_e = kr(I_n, ops_xi.e_r, I_eta); obj.e_s = kr(I_n, I_xi, ops_eta.e_l); - obj.e_n = kr(I_n, I_xi, - + obj.e_n = kr(I_n, I_xi, + metric_termsops_eta.e_r); - + [X,Y] = ti.map(obj.xi,obj.eta); - + [x_xi,x_eta] = gridDerivatives(X,ops_xi.D1,ops_eta.D1); [y_xi,y_eta] = gridDerivatives(Y,ops_xi.D1, ops_eta.D1); - + obj.X = reshape(X,m_tot,1); obj.Y = reshape(Y,m_tot,1); obj.X_xi = reshape(x_xi,m_tot,1); obj.Y_xi = reshape(y_xi,m_tot,1); obj.X_eta = reshape(x_eta,m_tot,1); obj.Y_eta = reshape(y_eta,m_tot,1); - + Ahat_evaluated = obj.evaluateCoefficientMatrix(obj.Ahat, obj.X, obj.Y,obj.X_eta,obj.Y_eta); Bhat_evaluated = obj.evaluateCoefficientMatrix(obj.Bhat, obj.X, obj.Y,obj.X_xi,obj.Y_xi); E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]); - + obj.m = m; obj.h = [ops_xi.h ops_eta.h]; obj.order = order; obj.J = obj.X_xi.*obj.Y_eta - obj.X_eta.*obj.Y_xi; obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); - + obj.D = obj.Ji*(-Ahat_evaluated*D1_xi-Bhat_evaluated*D1_eta)-E_evaluated; - + end - + % Closure functions return the opertors applied to the own doamin to close the boundary % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w',General boundary conditions'n','s'. @@ -127,18 +127,18 @@ error('No such boundary condition') end end - - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundaryGeneral boundary conditions) + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) error('An interface function does not exist yet'); end - + function N = size(obj) N = obj.m; end - + function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_) params = obj.params; - + if isa(mat,'function_handle') [rows,cols] = size(mat(params,0,0,0,0)); x_ = kr(obj.onesN,x_); @@ -152,7 +152,7 @@ side = max(length(X),length(Y)); cols = cols/side; end - + ret = cell(rows,cols); for ii = 1:rows for jj = 1:cols @@ -161,7 +161,7 @@ end ret = cell2mat(ret); end - + %Characteristic boundary conditions function [closure, penalty] = boundary_condition_char(obj,boundary) params = obj.params; @@ -169,7 +169,7 @@ Y = obj.Y; xi = obj.xi; eta = obj.eta; - + switch boundary case {'w','W','west'} e_ = obj.e_w; @@ -200,11 +200,11 @@ [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); side = max(length(xi)); end - + pos = signVec(1); zeroval = signVec(2); neg = signVec(3); - + switch boundPos case {'l'} tau = sparse(obj.n*side,pos); @@ -218,10 +218,10 @@ Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); closure = Hi*e_*V*tau*Vi_minus*e_'; penalty = -Hi*e_*V*tau*Vi_minus; - end + end end - - + + % General boundary condition in the form Lu=g(x) function [closure,penalty] = boundary_condition_general(obj,boundary,L) params = obj.params; @@ -229,7 +229,7 @@ Y = obj.Y; xi = obj.xi; eta = obj.eta; - + switch boundary case {'w','W','west'} e_ = obj.e_w; @@ -237,7 +237,7 @@ boundPos = 'l'; Hi = obj.Hxii; [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w)); - + Ji_vec = diag(obj.Ji); Ji = diag(Ji_vec(obj.index_w)); xi_x = Ji*obj.Y_eta(obj.index_w); @@ -250,7 +250,7 @@ boundPos = 'r'; Hi = obj.Hxii; [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); - + Ji_vec = diag(obj.Ji); Ji = diag(Ji_vec(obj.index_e)); xi_x = Ji*obj.Y_eta(obj.index_e); @@ -263,7 +263,7 @@ boundPos = 'l'; Hi = obj.Hetai; [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s)); - + Ji_vec = diag(obj.Ji); Ji = diag(Ji_vec(obj.index_s)); eta_x = Ji*obj.Y_xi(obj.index_s); @@ -276,7 +276,7 @@ boundPos = 'r'; Hi = obj.Hetai; [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); - + Ji_vec = diag(obj.Ji); Ji = diag(Ji_vec(obj.index_n)); eta_x = Ji*obj.Y_xi(obj.index_n); @@ -284,11 +284,11 @@ L = obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]); side = max(length(xi)); end - + pos = signVec(1); zeroval = signVec(2); neg = signVec(3); - + switch boundPos case {'l'} tau = sparse(obj.n*side,pos); @@ -296,7 +296,7 @@ Vi_minus = Vi(pos+1:obj.n*side,:); V_plus = V(:,1:pos); V_minus = V(:,(pos)+1:obj.n*side); - + tau(1:pos,:) = -abs(D(1:pos,1:pos)); R = -inv(L*V_plus)*(L*V_minus); closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; @@ -306,7 +306,7 @@ tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); Vi_plus = Vi(1:pos,:); Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); - + V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); R = -inv(L*V_minus)*(L*V_plus); @@ -314,7 +314,7 @@ penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; end end - + % Function that diagonalizes a symbolic matrix A as A=V*D*Vi % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign % [d+ ] @@ -329,22 +329,22 @@ else xs_ = 0; end - + if(sum(abs(y_))~= 0) syms ys_; else ys_ = 0; end - + [V, D] = eig(mat(params,xs,ys,xs_,ys_)); Vi = inv(V); syms xs ys xs_ ys_ - + xs = x; ys = y; xs_ = x_; ys_ = y_; - + side = max(length(x),length(y)); Dret = zeros(obj.n,side*obj.n); Vret = zeros(obj.n,side*obj.n); @@ -356,7 +356,7 @@ Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); end end - + D = sparse(Dret); V = sparse(Vret); Vi = sparse(Viret); @@ -364,11 +364,11 @@ D = obj.evaluateCoefficientMatrix(D,x,y,x_,y_); Vi = obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_); DD = diag(D); - + poseig = (DD>0); zeroeig = (DD==0); negeig = (DD<0); - + D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
--- a/+scheme/Hypsyst3d.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Hypsyst3d.m Thu Nov 22 22:03:44 2018 -0800 @@ -7,11 +7,11 @@ X, Y, Z% Values of x and y for each grid point Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces order % Order accuracy for the approximation - + D % non-stabalized scheme operator A, B, C, E % Symbolic coefficient matrices Aevaluated,Bevaluated,Cevaluated, Eevaluated - + H % Discrete norm Hx, Hy, Hz % Norms in the x, y and z directions Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. @@ -19,8 +19,8 @@ e_w, e_e, e_s, e_n, e_b, e_t params % Parameters for the coeficient matrice end - - + + methods % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Cu_z-Eu function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator) @@ -28,11 +28,11 @@ xlim = lim{1}; ylim = lim{2}; zlim = lim{3}; - + if length(m) == 1 m = [m m m]; end - + obj.A = A; obj.B = B; obj.C = C; @@ -41,7 +41,7 @@ m_y = m(2); m_z = m(3); obj.params = params; - + switch operator case 'upwind' ops_x = sbp.D1Upwind(m_x,xlim,order); @@ -52,29 +52,29 @@ ops_y = sbp.D2Standard(m_y,ylim,order); ops_z = sbp.D2Standard(m_z,zlim,order); end - + obj.x = ops_x.x; obj.y = ops_y.x; obj.z = ops_z.x; - + obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1)); obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1)); obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z); - + obj.Yx = kr(obj.y,ones(m_z,1)); obj.Zx = kr(ones(m_y,1),obj.z); obj.Xy = kr(obj.x,ones(m_z,1)); obj.Zy = kr(ones(m_x,1),obj.z); obj.Xz = kr(obj.x,ones(m_y,1)); obj.Yz = kr(ones(m_z,1),obj.y); - + obj.Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z); obj.Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z); obj.Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z); obj.Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z); - + obj.n = length(A(obj.params,0,0,0)); - + I_n = speye(obj.n); I_x = speye(m_x); obj.I_x = I_x; @@ -83,31 +83,31 @@ I_z = speye(m_z); obj.I_z = I_z; I_N = kr(I_n,I_x,I_y,I_z); - + obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z); obj.Hx = ops_x.H; obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z); obj.Hy = ops_y.H; obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI); obj.Hz = ops_z.H; - + obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z); obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z); obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z); obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z); obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l); obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r); - + obj.m = m; obj.h = [ops_x.h ops_y.h ops_x.h]; obj.order = order; - + switch operator case 'upwind' alphaA = max(abs(eig(A(params,obj.x(end),obj.y(end),obj.z(end))))); alphaB = max(abs(eig(B(params,obj.x(end),obj.y(end),obj.z(end))))); alphaC = max(abs(eig(C(params,obj.x(end),obj.y(end),obj.z(end))))); - + Ap = (obj.Aevaluated+alphaA*I_N)/2; Am = (obj.Aevaluated-alphaA*I_N)/2; Dpx = kr(I_n, ops_x.Dp, I_y,I_z); @@ -116,7 +116,7 @@ temp = Ap*Dmx; obj.D = obj.D-temp; clear Ap Am Dpx Dmx - + Bp = (obj.Bevaluated+alphaB*I_N)/2; Bm = (obj.Bevaluated-alphaB*I_N)/2; Dpy = kr(I_n, I_x, ops_y.Dp,I_z); @@ -126,20 +126,20 @@ temp = Bp*Dmy; obj.D = obj.D-temp; clear Bp Bm Dpy Dmy - - + + Cp = (obj.Cevaluated+alphaC*I_N)/2; Cm = (obj.Cevaluated-alphaC*I_N)/2; Dpz = kr(I_n, I_x, I_y,ops_z.Dp); Dmz = kr(I_n, I_x, I_y,ops_z.Dm); - + temp = Cm*Dpz; obj.D = obj.D-temp; temp = Cp*Dmz; obj.D = obj.D-temp; clear Cp Cm Dpz Dmz obj.D = obj.D-obj.Eevaluated; - + case 'standard' D1_x = kr(I_n, ops_x.D1, I_y,I_z); D1_y = kr(I_n, I_x, ops_y.D1,I_z); @@ -149,7 +149,7 @@ error('Opperator not supported'); end end - + % Closure functions return the opertors applied to the own doamin to close the boundary % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. @@ -167,15 +167,15 @@ error('No such boundary condition') end end - - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) error('An interface function does not exist yet'); end - + function N = size(obj) N = obj.m; end - + function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z) params = obj.params; side = max(length(X),length(Y)); @@ -189,7 +189,7 @@ side = max(length(X),length(Y)); cols = cols/side; end - + ret = cell(rows,cols); for ii = 1:rows for jj = 1:cols @@ -198,10 +198,10 @@ end ret = cell2mat(ret); end - + function [BM] = boundary_matrices(obj,boundary) params = obj.params; - + switch boundary case {'w','W','west'} BM.e_ = obj.e_w; @@ -248,7 +248,7 @@ end BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); end - + % Characteristic bouyndary consitions function [closure, penalty]=boundary_condition_char(obj,BM) side = BM.side; @@ -260,7 +260,7 @@ Hi = BM.Hi; D = BM.D; e_ = BM.e_; - + switch BM.boundpos case {'l'} tau = sparse(obj.n*side,pos); @@ -276,9 +276,9 @@ penalty = -Hi*e_*V*tau*Vi_minus; end end - + % General boundary condition in the form Lu=g(x) - function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L) + function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L) side = BM.side; pos = BM.pos; neg = BM.neg; @@ -288,7 +288,7 @@ Hi = BM.Hi; D = BM.D; e_ = BM.e_; - + switch boundary case {'w','W','west'} L = obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx); @@ -303,7 +303,7 @@ case {'t','T','top'} L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end)); end - + switch BM.boundpos case {'l'} tau = sparse(obj.n*side,pos); @@ -311,7 +311,7 @@ Vi_minus = Vi(pos+zeroval+1:obj.n*side,:); V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); - + tau(1:pos,:) = -abs(D(1:pos,1:pos)); R = -inv(L*V_plus)*(L*V_minus); closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; @@ -321,7 +321,7 @@ tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); Vi_plus = Vi(1:pos,:); Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); - + V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); R = -inv(L*V_minus)*(L*V_plus); @@ -329,7 +329,7 @@ penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; end end - + % Function that diagonalizes a symbolic matrix A as A=V*D*Vi % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign % [d+ ] @@ -344,13 +344,13 @@ xs = x; ys = y; zs = z; - - + + side = max(length(x),length(y)); Dret = zeros(obj.n,side*obj.n); Vret = zeros(obj.n,side*obj.n); Viret= zeros(obj.n,side*obj.n); - + for ii=1:obj.n for jj=1:obj.n Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); @@ -358,7 +358,7 @@ Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); end end - + D = sparse(Dret); V = sparse(Vret); Vi = sparse(Viret); @@ -366,11 +366,11 @@ Vi= obj.evaluateCoefficientMatrix(Vi,x,y,z); D = obj.evaluateCoefficientMatrix(D,x,y,z); DD = diag(D); - + poseig = (DD>0); zeroeig = (DD==0); negeig = (DD<0); - + D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; Vi= [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
--- a/+scheme/Hypsyst3dCurve.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Hypsyst3dCurve.m Thu Nov 22 22:03:44 2018 -0800 @@ -5,22 +5,22 @@ h % Grid spacing X, Y, Z% Values of x and y for each grid point Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces - + xi,eta,zeta Xi, Eta, Zeta - + Eta_xi, Zeta_xi, Xi_eta, Zeta_eta, Xi_zeta, Eta_zeta % Metric terms X_xi, X_eta, X_zeta,Y_xi,Y_eta,Y_zeta,Z_xi,Z_eta,Z_zeta % Metric terms - + order % Order accuracy for the approximation - + D % non-stabalized scheme operator Aevaluated, Bevaluated, Cevaluated, Eevaluated % Numeric Coeffiecient matrices Ahat, Bhat, Chat % Symbolic Transformed Coefficient matrices A, B, C, E % Symbolic coeffiecient matrices - + J, Ji % JAcobian and inverse Jacobian - + H % Discrete norm % Norms in the x, y and z directions Hxii,Hetai,Hzetai, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. @@ -30,14 +30,14 @@ index_w, index_e,index_s,index_n, index_b, index_t params %parameters for the coeficient matrice end - - + + methods function obj = Hypsyst3dCurve(m, order, A, B,C, E, params,ti,operator) xilim ={0 1}; etalim = {0 1}; zetalim = {0 1}; - + if length(m) == 1 m = [m m m]; end @@ -47,11 +47,11 @@ m_tot = m_xi*m_eta*m_zeta; obj.params = params; obj.n = length(A(obj,0,0,0)); - + obj.m = m; obj.order = order; obj.onesN = ones(obj.n); - + switch operator case 'upwind' ops_xi = sbp.D1Upwind(m_xi,xilim,order); @@ -64,21 +64,21 @@ otherwise error('Operator not available') end - + obj.xi = ops_xi.x; obj.eta = ops_eta.x; obj.zeta = ops_zeta.x; - + obj.Xi = kr(obj.xi,ones(m_eta,1),ones(m_zeta,1)); obj.Eta = kr(ones(m_xi,1),obj.eta,ones(m_zeta,1)); obj.Zeta = kr(ones(m_xi,1),ones(m_eta,1),obj.zeta); - - + + [X,Y,Z] = ti.map(obj.Xi,obj.Eta,obj.Zeta); obj.X = X; obj.Y = Y; obj.Z = Z; - + I_n = eye(obj.n); I_xi = speye(m_xi); obj.I_xi = I_xi; @@ -86,19 +86,19 @@ obj.I_eta = I_eta; I_zeta = speye(m_zeta); obj.I_zeta = I_zeta; - + I_N=kr(I_n,I_xi,I_eta,I_zeta); - + O_xi = ones(m_xi,1); O_eta = ones(m_eta,1); O_zeta = ones(m_zeta,1); - - + + obj.Hxi = ops_xi.H; obj.Heta = ops_eta.H; obj.Hzeta = ops_zeta.H; obj.h = [ops_xi.h ops_eta.h ops_zeta.h]; - + switch operator case 'upwind' D1_xi = kr((ops_xi.Dp+ops_xi.Dm)/2, I_eta,I_zeta); @@ -109,11 +109,11 @@ D1_eta = kr(I_xi, ops_eta.D1,I_zeta); D1_zeta = kr(I_xi, I_eta,ops_zeta.D1); end - + obj.A = A; obj.B = B; obj.C = C; - + obj.X_xi = D1_xi*X; obj.X_eta = D1_eta*X; obj.X_zeta = D1_zeta*X; @@ -123,55 +123,55 @@ obj.Z_xi = D1_xi*Z; obj.Z_eta = D1_eta*Z; obj.Z_zeta = D1_zeta*Z; - + obj.Ahat = @transform_coefficient_matrix; obj.Bhat = @transform_coefficient_matrix; obj.Chat = @transform_coefficient_matrix; obj.E = @(obj,x,y,z,~,~,~,~,~,~)E(obj,x,y,z); - + obj.Aevaluated = obj.evaluateCoefficientMatrix(obj.Ahat,obj.X, obj.Y,obj.Z, obj.X_eta,obj.X_zeta,obj.Y_eta,obj.Y_zeta,obj.Z_eta,obj.Z_zeta); obj.Bevaluated = obj.evaluateCoefficientMatrix(obj.Bhat,obj.X, obj.Y,obj.Z, obj.X_zeta,obj.X_xi,obj.Y_zeta,obj.Y_xi,obj.Z_zeta,obj.Z_xi); obj.Cevaluated = obj.evaluateCoefficientMatrix(obj.Chat,obj.X,obj.Y,obj.Z, obj.X_xi,obj.X_eta,obj.Y_xi,obj.Y_eta,obj.Z_xi,obj.Z_eta); - + switch operator case 'upwind' clear D1_xi D1_eta D1_zeta alphaA = max(abs(eig(obj.Ahat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_eta(end),obj.X_zeta(end),obj.Y_eta(end),obj.Y_zeta(end),obj.Z_eta(end),obj.Z_zeta(end))))); alphaB = max(abs(eig(obj.Bhat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_zeta(end),obj.X_xi(end),obj.Y_zeta(end),obj.Y_xi(end),obj.Z_zeta(end),obj.Z_xi(end))))); alphaC = max(abs(eig(obj.Chat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_xi(end),obj.X_eta(end),obj.Y_xi(end),obj.Y_eta(end),obj.Z_xi(end),obj.Z_eta(end))))); - + Ap = (obj.Aevaluated+alphaA*I_N)/2; Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta); diffSum = -Ap*Dmxi; clear Ap Dmxi - + Am = (obj.Aevaluated-alphaA*I_N)/2; - + obj.Aevaluated = []; Dpxi = kr(I_n, ops_xi.Dp, I_eta,I_zeta); temp = Am*Dpxi; diffSum = diffSum-temp; clear Am Dpxi - + Bp = (obj.Bevaluated+alphaB*I_N)/2; Dmeta = kr(I_n, I_xi, ops_eta.Dm,I_zeta); temp = Bp*Dmeta; diffSum = diffSum-temp; clear Bp Dmeta - + Bm = (obj.Bevaluated-alphaB*I_N)/2; obj.Bevaluated = []; Dpeta = kr(I_n, I_xi, ops_eta.Dp,I_zeta); temp = Bm*Dpeta; diffSum = diffSum-temp; clear Bm Dpeta - + Cp = (obj.Cevaluated+alphaC*I_N)/2; Dmzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dm); temp = Cp*Dmzeta; diffSum = diffSum-temp; clear Cp Dmzeta - + Cm = (obj.Cevaluated-alphaC*I_N)/2; clear I_N obj.Cevaluated = []; @@ -179,72 +179,72 @@ temp = Cm*Dpzeta; diffSum = diffSum-temp; clear Cm Dpzeta temp - + obj.J = obj.X_xi.*obj.Y_eta.*obj.Z_zeta... +obj.X_zeta.*obj.Y_xi.*obj.Z_eta... +obj.X_eta.*obj.Y_zeta.*obj.Z_xi... -obj.X_xi.*obj.Y_zeta.*obj.Z_eta... -obj.X_eta.*obj.Y_xi.*obj.Z_zeta... -obj.X_zeta.*obj.Y_eta.*obj.Z_xi; - + obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); - + obj.D = obj.Ji*diffSum-obj.Eevaluated; - + case 'standard' D1_xi = kr(I_n,D1_xi); D1_eta = kr(I_n,D1_eta); D1_zeta = kr(I_n,D1_zeta); - + obj.J = obj.X_xi.*obj.Y_eta.*obj.Z_zeta... +obj.X_zeta.*obj.Y_xi.*obj.Z_eta... +obj.X_eta.*obj.Y_zeta.*obj.Z_xi... -obj.X_xi.*obj.Y_zeta.*obj.Z_eta... -obj.X_eta.*obj.Y_xi.*obj.Z_zeta... -obj.X_zeta.*obj.Y_eta.*obj.Z_xi; - + obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); - + obj.D = obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated; otherwise error('Operator not supported') end - + obj.Hxii = kr(I_n, ops_xi.HI, I_eta,I_zeta); obj.Hetai = kr(I_n, I_xi, ops_eta.HI,I_zeta); obj.Hzetai = kr(I_n, I_xi,I_eta, ops_zeta.HI); - + obj.index_w = (kr(ops_xi.e_l, O_eta,O_zeta)==1); obj.index_e = (kr(ops_xi.e_r, O_eta,O_zeta)==1); obj.index_s = (kr(O_xi, ops_eta.e_l,O_zeta)==1); obj.index_n = (kr(O_xi, ops_eta.e_r,O_zeta)==1); obj.index_b = (kr(O_xi, O_eta, ops_zeta.e_l)==1); obj.index_t = (kr(O_xi, O_eta, ops_zeta.e_r)==1); - + obj.e_w = kr(I_n, ops_xi.e_l, I_eta,I_zeta); obj.e_e = kr(I_n, ops_xi.e_r, I_eta,I_zeta); obj.e_s = kr(I_n, I_xi, ops_eta.e_l,I_zeta); obj.e_n = kr(I_n, I_xi, ops_eta.e_r,I_zeta); obj.e_b = kr(I_n, I_xi, I_eta, ops_zeta.e_l); obj.e_t = kr(I_n, I_xi, I_eta, ops_zeta.e_r); - + obj.Eta_xi = kr(obj.eta,ones(m_xi,1)); obj.Zeta_xi = kr(ones(m_eta,1),obj.zeta); obj.Xi_eta = kr(obj.xi,ones(m_zeta,1)); obj.Zeta_eta = kr(ones(m_xi,1),obj.zeta); obj.Xi_zeta = kr(obj.xi,ones(m_eta,1)); - obj.Eta_zeta = kr(ones(m_zeta,1),obj.eta); + obj.Eta_zeta = kr(ones(m_zeta,1),obj.eta); end - + function [ret] = transform_coefficient_matrix(obj,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2) ret = obj.A(obj,x,y,z).*(y_1.*z_2-z_1.*y_2); ret = ret+obj.B(obj,x,y,z).*(x_2.*z_1-x_1.*z_2); ret = ret+obj.C(obj,x,y,z).*(x_1.*y_2-x_2.*y_1); end - - + + % Closure functions return the opertors applied to the own doamin to close the boundary % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. @@ -253,7 +253,7 @@ function [closure, penalty] = boundary_condition(obj,boundary,type,L) default_arg('type','char'); BM = boundary_matrices(obj,boundary); - + switch type case{'c','char'} [closure,penalty] = boundary_condition_char(obj,BM); @@ -263,15 +263,15 @@ error('No such boundary condition') end end - - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) error('An interface function does not exist yet'); end - + function N = size(obj) N = obj.m; end - + % Evaluates the symbolic Coeffiecient matrix mat function [ret] = evaluateCoefficientMatrix(obj,mat, X, Y, Z , x_1 , x_2 , y_1 , y_2 , z_1 , z_2) params = obj.params; @@ -294,7 +294,7 @@ end matVec(abs(matVec)<10^(-10)) = 0; ret = cell(rows,cols); - + for ii = 1:rows for jj = 1:cols ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); @@ -302,7 +302,7 @@ end ret = cell2mat(ret); end - + function [BM] = boundary_matrices(obj,boundary) params = obj.params; BM.boundary = boundary; @@ -385,7 +385,7 @@ BM.side = sum(BM.index); BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); end - + % Characteristic boundary condition function [closure, penalty] = boundary_condition_char(obj,BM) side = BM.side; @@ -397,7 +397,7 @@ Hi = BM.Hi; D = BM.D; e_ = BM.e_; - + switch BM.boundpos case {'l'} tau = sparse(obj.n*side,pos); @@ -413,7 +413,7 @@ penalty = -Hi*e_*V*tau*Vi_minus; end end - + % General boundary condition in the form Lu=g(x) function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L) side = BM.side; @@ -426,7 +426,7 @@ D = BM.D; e_ = BM.e_; index = BM.index; - + switch BM.boundary case{'b','B','bottom'} Ji_vec = diag(obj.Ji); @@ -434,10 +434,10 @@ Zeta_x = Ji*(obj.Y_xi(index).*obj.Z_eta(index)-obj.Z_xi(index).*obj.Y_eta(index)); Zeta_y = Ji*(obj.X_eta(index).*obj.Z_xi(index)-obj.X_xi(index).*obj.Z_eta(index)); Zeta_z = Ji*(obj.X_xi(index).*obj.Y_eta(index)-obj.Y_xi(index).*obj.X_eta(index)); - + L = obj.evaluateCoefficientMatrix(L,Zeta_x,Zeta_y,Zeta_z,[],[],[],[],[],[]); end - + switch BM.boundpos case {'l'} tau = sparse(obj.n*side,pos); @@ -445,7 +445,7 @@ Vi_minus = Vi(pos+zeroval+1:obj.n*side,:); V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); - + tau(1:pos,:) = -abs(D(1:pos,1:pos)); R = -inv(L*V_plus)*(L*V_minus); closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; @@ -455,7 +455,7 @@ tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); Vi_plus = Vi(1:pos,:); Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); - + V_plus = V(:,1:pos); V_minus = V(:,(pos+zeroval)+1:obj.n*side); R = -inv(L*V_minus)*(L*V_plus); @@ -463,7 +463,7 @@ penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; end end - + % Function that diagonalizes a symbolic matrix A as A=V*D*Vi % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign % [d+ ] @@ -478,38 +478,38 @@ else x_1s = 0; end - + if(sum(abs(x_2))>eps) syms x_2s; else x_2s = 0; end - - + + if(sum(abs(y_1))>eps) syms y_1s else y_1s = 0; end - + if(sum(abs(y_2))>eps) syms y_2s; else y_2s = 0; end - + if(sum(abs(z_1))>eps) syms z_1s else z_1s = 0; end - + if(sum(abs(z_2))>eps) syms z_2s; else z_2s = 0; end - + syms xs ys zs [V, D] = eig(mat(obj,xs,ys,zs,x_1s,x_2s,y_1s,y_2s,z_1s,z_2s)); Vi = inv(V); @@ -522,12 +522,12 @@ y_2s = y_2; z_1s = z_1; z_2s = z_2; - + side = max(length(x),length(y)); Dret = zeros(obj.n,side*obj.n); Vret = zeros(obj.n,side*obj.n); Viret = zeros(obj.n,side*obj.n); - + for ii=1:obj.n for jj=1:obj.n Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); @@ -535,7 +535,7 @@ Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); end end - + D = sparse(Dret); V = sparse(Vret); Vi = sparse(Viret); @@ -543,11 +543,11 @@ D = obj.evaluateCoefficientMatrix(D,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); Vi = obj.evaluateCoefficientMatrix(Vi,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); DD = diag(D); - + poseig = (DD>0); zeroeig = (DD==0); negeig = (DD<0); - + D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
--- a/+scheme/Scheme.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Scheme.m Thu Nov 22 22:03:44 2018 -0800 @@ -26,7 +26,12 @@ % interface to. % penalty may be a cell array if there are several penalties with different weights [closure, penalty] = boundary_condition(obj,boundary,type) % TODO: Change name to boundaryCondition - [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + % type Specifies the type of interface coupling. + % The format of type is different for every scheme. + % Some schemes may only have one type implemented, in which case + % the input argument is a dummy. + [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) % TODO: op = getBoundaryOperator()?? % makes sense to have it available through a method instead of random properties
--- a/+scheme/Schrodinger.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Schrodinger.m Thu Nov 22 22:03:44 2018 -0800 @@ -90,7 +90,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain [e_u,d_u,s_u] = obj.get_boundary_ops(boundary);
--- a/+scheme/Schrodinger2d.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Schrodinger2d.m Thu Nov 22 22:03:44 2018 -0800 @@ -1,9 +1,9 @@ classdef Schrodinger2d < scheme.Scheme % Discretizes the Laplacian with constant coefficent, -% in the Schrödinger equation way (i.e., the discretization matrix is not necessarily +% in the Schrödinger equation way (i.e., the discretization matrix is not necessarily % definite) -% u_t = a*i*Laplace u +% u_t = a*i*Laplace u % opSet should be cell array of opSets, one per dimension. This % is useful if we have periodic BC in one direction. @@ -30,7 +30,7 @@ d1_l, d1_r % Normal derivatives at the boundary e_w, e_e, e_s, e_n d_w, d_e, d_s, d_n - + H_boundary % Boundary inner products interpolation_type % MC or AWW @@ -167,7 +167,7 @@ default_arg('parameter', []); % j is the coordinate direction of the boundary - % nj: outward unit normal component. + % nj: outward unit normal component. % nj = -1 for west, south, bottom boundaries % nj = 1 for east, north, top boundaries [j, nj] = obj.get_boundary_number(boundary); @@ -188,13 +188,13 @@ % Dirichlet boundary condition case {'D','d','dirichlet','Dirichlet'} - closure = nj*Hi*d{j}*a*1i*H_gamma*(e{j}' ); + closure = nj*Hi*d{j}*a*1i*H_gamma*(e{j}' ); penalty = -nj*Hi*d{j}*a*1i*H_gamma; % Free boundary condition case {'N','n','neumann','Neumann'} - closure = -nj*Hi*e{j}*a*1i*H_gamma*(d{j}' ); - penalty = nj*Hi*e{j}*a*1i*H_gamma; + closure = -nj*Hi*e{j}*a*1i*H_gamma*(d{j}' ); + penalty = nj*Hi*e{j}*a*1i*H_gamma; % Unknown boundary condition otherwise @@ -202,7 +202,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain % Get neighbour boundary operator @@ -251,7 +251,7 @@ d = obj.d1_l; end e = e{coord}; - d = d{coord}; + d = d{coord}; Hi = obj.Hi; sigma = -n*1i*a/2; @@ -281,15 +281,15 @@ case 'AWW' %String 'C2F' indicates that ICF2 is more accurate. interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); - interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); - if grid_ratio < 1 + interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); + if grid_ratio < 1 % Local is coarser than neighbour I_neighbour2local_e = interpOpSetF2C.IF2C; I_neighbour2local_d = interpOpSetC2F.IF2C; I_local2neighbour_e = interpOpSetC2F.IC2F; I_local2neighbour_d = interpOpSetF2C.IC2F; elseif grid_ratio > 1 - % Local is finer than neighbour + % Local is finer than neighbour I_neighbour2local_e = interpOpSetC2F.IC2F; I_neighbour2local_d = interpOpSetF2C.IC2F; I_local2neighbour_e = interpOpSetF2C.IF2C; @@ -300,7 +300,7 @@ ' is not available.' ]); end - else + else % No interpolation required I_neighbour2local_e = speye(m,m); I_neighbour2local_d = speye(m,m); @@ -310,8 +310,8 @@ closure = tau*Hi*d*H_gamma*e' + sigma*Hi*e*H_gamma*d'; penalty = -tau*Hi*d*H_gamma*I_neighbour2local_e*e_neighbour' ... - -sigma*Hi*e*H_gamma*I_neighbour2local_d*d_neighbour'; - + -sigma*Hi*e*H_gamma*I_neighbour2local_d*d_neighbour'; + end % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
--- a/+scheme/Utux.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Utux.m Thu Nov 22 22:03:44 2018 -0800 @@ -16,15 +16,15 @@ end - methods + methods function obj = Utux(g ,order, operator) default_arg('operator','Standard'); - + m = g.size(); xl = g.getBoundary('l'); xr = g.getBoundary('r'); xlim = {xl, xr}; - + switch operator % case 'NonEquidistant' % ops = sbp.D1Nonequidistant(m,xlim,order); @@ -38,12 +38,12 @@ otherwise error('Unvalid operator') end - + obj.grid = g; obj.H = ops.H; obj.Hi = ops.HI; - + obj.e_l = ops.e_l; obj.e_r = ops.e_r; obj.D = -obj.D1; @@ -62,27 +62,27 @@ % neighbour_boundary is a string specifying which boundary to interface to. function [closure, penalty] = boundary_condition(obj,boundary,type) default_arg('type','dirichlet'); - tau =-1*obj.e_l; - closure = obj.Hi*tau*obj.e_l'; + tau =-1*obj.e_l; + closure = obj.Hi*tau*obj.e_l'; penalty = -obj.Hi*tau; - + end - - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) switch boundary % Upwind coupling case {'l','left'} tau = -1*obj.e_l; - closure = obj.Hi*tau*obj.e_l'; + closure = obj.Hi*tau*obj.e_l'; penalty = -obj.Hi*tau*neighbour_scheme.e_r'; case {'r','right'} tau = 0*obj.e_r; - closure = obj.Hi*tau*obj.e_r'; + closure = obj.Hi*tau*obj.e_r'; penalty = -obj.Hi*tau*neighbour_scheme.e_l'; end - + end - + function N = size(obj) N = obj.m; end
--- a/+scheme/Utux2D.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Utux2D.m Thu Nov 22 22:03:44 2018 -0800 @@ -5,7 +5,7 @@ grid % Grid order % Order accuracy for the approximation v0 % Initial data - + a % Wave speed a = [a1, a2]; % Can either be a constant vector or a cell array of function handles. @@ -15,36 +15,36 @@ % Derivatives Dx, Dy - + % Boundary operators e_w, e_e, e_s, e_n - + D % Total discrete operator % String, type of interface coupling % Default: 'upwind' % Other: 'centered' - coupling_type + coupling_type % String, type of interpolation operators % Default: 'AWW' (Almquist Wang Werpers) % Other: 'MC' (Mattsson Carpenter) interpolation_type - + % Cell array, damping on upwstream and downstream sides. interpolation_damping end - methods + methods function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type, interpolation_damping) - + default_arg('interpolation_damping',{0,0}); - default_arg('interpolation_type','AWW'); - default_arg('coupling_type','upwind'); - default_arg('a',1/sqrt(2)*[1, 1]); + default_arg('interpolation_type','AWW'); + default_arg('coupling_type','upwind'); + default_arg('a',1/sqrt(2)*[1, 1]); default_arg('opSet',@sbp.D2Standard); assert(isa(g, 'grid.Cartesian')) @@ -55,7 +55,7 @@ else a = {a(1), a(2)}; end - + m = g.size(); m_x = m(1); m_y = m(2); @@ -70,13 +70,13 @@ ops_y = opSet(m_y, ylim, order); Ix = speye(m_x); Iy = speye(m_y); - + % Norms Hx = ops_x.H; Hy = ops_y.H; Hxi = ops_x.HI; Hyi = ops_y.HI; - + obj.H_x = Hx; obj.H_y = Hy; obj.H = kron(Hx,Hy); @@ -85,13 +85,13 @@ obj.Hy = kron(Ix,Hy); obj.Hxi = kron(Hxi,Iy); obj.Hyi = kron(Ix,Hyi); - + % Derivatives Dx = ops_x.D1; Dy = ops_y.D1; obj.Dx = kron(Dx,Iy); obj.Dy = kron(Ix,Dy); - + % Boundary operators obj.e_w = kr(ops_x.e_l, Iy); obj.e_e = kr(ops_x.e_r, Iy); @@ -117,23 +117,23 @@ % neighbour_boundary is a string specifying which boundary to interface to. function [closure, penalty] = boundary_condition(obj,boundary,type) default_arg('type','dirichlet'); - + sigma = -1; % Scalar penalty parameter switch boundary case {'w','W','west','West'} tau = sigma*obj.a{1}*obj.e_w*obj.H_y; closure = obj.Hi*tau*obj.e_w'; - + case {'s','S','south','South'} tau = sigma*obj.a{2}*obj.e_s*obj.H_x; closure = obj.Hi*tau*obj.e_s'; - end + end penalty = -obj.Hi*tau; - + end - - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) + % Get neighbour boundary operator switch neighbour_boundary case {'e','E','east','East'} @@ -149,9 +149,9 @@ e_neighbour = neighbour_scheme.e_s; m_neighbour = neighbour_scheme.m(1); end - + switch obj.coupling_type - + % Upwind coupling (energy dissipation) case 'upwind' sigma_ds = -1; %"Downstream" penalty @@ -169,7 +169,7 @@ % Check grid ratio for interpolation switch boundary case {'w','W','west','West','e','E','east','East'} - m = obj.m(2); + m = obj.m(2); case {'s','S','south','South','n','N','north','North'} m = obj.m(1); end @@ -197,15 +197,15 @@ case 'AWW' %String 'C2F' indicates that ICF2 is more accurate. interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); - interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); - if grid_ratio < 1 + interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); + if grid_ratio < 1 % Local is coarser than neighbour I_neighbour2local_us = interpOpSetC2F.IF2C; I_neighbour2local_ds = interpOpSetF2C.IF2C; I_local2neighbour_us = interpOpSetC2F.IC2F; I_local2neighbour_ds = interpOpSetF2C.IC2F; elseif grid_ratio > 1 - % Local is finer than neighbour + % Local is finer than neighbour I_neighbour2local_us = interpOpSetF2C.IC2F; I_neighbour2local_ds = interpOpSetC2F.IC2F; I_local2neighbour_us = interpOpSetF2C.IF2C; @@ -216,12 +216,12 @@ ' is not available.' ]); end - else + else % No interpolation required I_neighbour2local_us = speye(m,m); I_neighbour2local_ds = speye(m,m); - end - + end + int_damp_us = obj.interpolation_damping{1}; int_damp_ds = obj.interpolation_damping{2}; @@ -238,7 +238,7 @@ beta = int_damp_ds*obj.a{1}... *obj.e_w*obj.H_y; - closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w'; + closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w'; case {'e','E','east','East'} tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y; closure = obj.Hi*tau*obj.e_e'; @@ -246,10 +246,10 @@ beta = int_damp_us*obj.a{1}... *obj.e_e*obj.H_y; - closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e'; + closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e'; case {'s','S','south','South'} tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x; - closure = obj.Hi*tau*obj.e_s'; + closure = obj.Hi*tau*obj.e_s'; penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; beta = int_damp_ds*obj.a{2}... @@ -262,12 +262,12 @@ beta = int_damp_us*obj.a{2}... *obj.e_n*obj.H_x; - closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n'; + closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n'; end - - + + end - + function N = size(obj) N = obj.m; end
--- a/+scheme/Wave.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Wave.m Thu Nov 22 22:03:44 2018 -0800 @@ -112,7 +112,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain [e_u,d_u,s_u] = obj.get_boundary_ops(boundary);
--- a/+scheme/Wave2d.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Wave2d.m Thu Nov 22 22:03:44 2018 -0800 @@ -158,7 +158,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain [e_u,d_u,s_u,gamm_u, halfnorm_inv] = obj.get_boundary_ops(boundary);
--- a/+scheme/Wave2dCurve.m Thu Nov 22 22:03:06 2018 -0800 +++ b/+scheme/Wave2dCurve.m Thu Nov 22 22:03:44 2018 -0800 @@ -243,7 +243,7 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain tuning = 1.2;