Mercurial > repos > public > sbplib
changeset 609:8cbecf22075b feature/utux2D
Merge to get interpolation operators.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sat, 14 Oct 2017 22:36:31 -0700 |
| parents | 0f9d20dbb7ce (diff) c923fe6197ff (current diff) |
| children | b7b3c11fab4d |
| files | |
| diffstat | 7 files changed, 610 insertions(+), 127 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/+domain/Circle.m Sat Oct 14 22:36:31 2017 -0700 @@ -0,0 +1,98 @@ +classdef Circle < multiblock.DefCurvilinear + properties + r, c + + hs + r_arc + omega + end + + methods + function obj = Circle(r, c, hs) + default_arg('r', 1); + default_arg('c', [0; 0]); + default_arg('hs', 0.435); + + + % alpha = 0.75; + % hs = alpha*r/sqrt(2); + + % Square should not be a square, it should be an arc. The arc radius + % is chosen so that the three angles of the meshes are all equal. + % This gives that the (half)arc opening angle of should be omega = pi/12 + omega = pi/12; + r_arc = hs*(2*sqrt(2))/(sqrt(3)-1); % = hs* 1/sin(omega) + c_arc = c - [(1/(2-sqrt(3))-1)*hs; 0]; + + cir = parametrization.Curve.circle(c,r,[-pi/4 pi/4]); + + c2 = cir(0); + c3 = cir(1); + + s1 = [-hs; -hs]; + s2 = [ hs; -hs]; + s3 = [ hs; hs]; + s4 = [-hs; hs]; + + sp2 = parametrization.Curve.line(s2,c2); + sp3 = parametrization.Curve.line(s3,c3); + + Se1 = parametrization.Curve.circle(c_arc,r_arc,[-omega, omega]); + Se2 = Se1.rotate(c,pi/2); + Se3 = Se2.rotate(c,pi/2); + Se4 = Se3.rotate(c,pi/2); + + + S = parametrization.Ti(Se1,Se2,Se3,Se4).rotate_edges(-1); + + A = parametrization.Ti(sp2, cir, sp3.reverse, Se1.reverse); + B = A.rotate(c,1*pi/2).rotate_edges(-1); + C = A.rotate(c,2*pi/2).rotate_edges(-1); + D = A.rotate(c,3*pi/2).rotate_edges(0); + + blocks = {S,A,B,C,D}; + blocksNames = {'S','A','B','C','D'}; + + conn = cell(5,5); + conn{1,2} = {'e','w'}; + conn{1,3} = {'n','s'}; + conn{1,4} = {'w','s'}; + conn{1,5} = {'s','w'}; + + conn{2,3} = {'n','e'}; + conn{3,4} = {'w','e'}; + conn{4,5} = {'w','s'}; + conn{5,2} = {'n','s'}; + + boundaryGroups = struct(); + boundaryGroups.E = multiblock.BoundaryGroup({2,'e'}); + boundaryGroups.N = multiblock.BoundaryGroup({3,'n'}); + boundaryGroups.W = multiblock.BoundaryGroup({4,'n'}); + boundaryGroups.S = multiblock.BoundaryGroup({5,'e'}); + boundaryGroups.all = multiblock.BoundaryGroup({{2,'e'},{3,'n'},{4,'n'},{5,'e'}}); + + obj = obj@multiblock.DefCurvilinear(blocks, conn, boundaryGroups, blocksNames); + + obj.r = r; + obj.c = c; + obj.hs = hs; + obj.r_arc = r_arc; + obj.omega = omega; + end + + function ms = getGridSizes(obj, m) + m_S = m; + + % m_Radial + s = 2*obj.hs; + innerArc = obj.r_arc*obj.omega; + outerArc = obj.r*pi/2; + shortSpoke = obj.r-s/sqrt(2); + x = (1/(2-sqrt(3))-1)*obj.hs; + longSpoke = (obj.r+x)-obj.r_arc; + m_R = parametrization.equal_step_size((innerArc+outerArc)/2, m_S, (shortSpoke+longSpoke)/2); + + ms = {[m_S m_S], [m_R m_S], [m_S m_R], [m_S m_R], [m_R m_S]}; + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/+domain/Rectangle.m Sat Oct 14 22:36:31 2017 -0700 @@ -0,0 +1,188 @@ +classdef Rectangle < multiblock.Definition + properties + + blockTi % Transfinite interpolation objects used for plotting + xlims + ylims + blockNames % Cell array of block labels + nBlocks + connections % Cell array specifying connections between blocks + boundaryGroups % Structure of boundaryGroups + + end + + + methods + % Creates a divided rectangle + % x and y are vectors of boundary and interface positions. + % blockNames: cell array of labels. The id is default. + function obj = Rectangle(x,y,blockNames) + default_arg('blockNames',[]); + + n = length(y)-1; % number of blocks in the y direction. + m = length(x)-1; % number of blocks in the x direction. + N = n*m; % number of blocks + + if ~issorted(x) + error('The elements of x seem to be in the wrong order'); + end + if ~issorted(flip(y)) + error('The elements of y seem to be in the wrong order'); + end + + % Dimensions of blocks and number of points + blockTi = cell(N,1); + xlims = cell(N,1); + ylims = cell(N,1); + for i = 1:n + for j = 1:m + p1 = [x(j), y(i+1)]; + p2 = [x(j+1), y(i)]; + I = flat_index(m,j,i); + blockTi{I} = parametrization.Ti.rectangle(p1,p2); + xlims{I} = {x(j), x(j+1)}; + ylims{I} = {y(i+1), y(i)}; + end + end + + % Interface couplings + conn = cell(N,N); + for i = 1:n + for j = 1:m + I = flat_index(m,j,i); + if i < n + J = flat_index(m,j,i+1); + conn{I,J} = {'s','n'}; + end + + if j < m + J = flat_index(m,j+1,i); + conn{I,J} = {'e','w'}; + end + end + end + + % Block names (id number as default) + if isempty(blockNames) + obj.blockNames = cell(1, N); + for i = 1:N + obj.blockNames{i} = sprintf('%d', i); + end + else + assert(length(blockNames) == N); + obj.blockNames = blockNames; + end + nBlocks = N; + + % Boundary groups + boundaryGroups = struct(); + nx = m; + ny = n; + E = cell(1,ny); + W = cell(1,ny); + S = cell(1,nx); + N = cell(1,nx); + for i = 1:ny + E_id = flat_index(m,nx,i); + W_id = flat_index(m,1,i); + E{i} = {E_id,'e'}; + W{i} = {W_id,'w'}; + end + for j = 1:nx + S_id = flat_index(m,j,ny); + N_id = flat_index(m,j,1); + S{j} = {S_id,'s'}; + N{j} = {N_id,'n'}; + end + boundaryGroups.E = multiblock.BoundaryGroup(E); + boundaryGroups.W = multiblock.BoundaryGroup(W); + boundaryGroups.S = multiblock.BoundaryGroup(S); + boundaryGroups.N = multiblock.BoundaryGroup(N); + boundaryGroups.all = multiblock.BoundaryGroup([E,W,S,N]); + boundaryGroups.WS = multiblock.BoundaryGroup([W,S]); + boundaryGroups.WN = multiblock.BoundaryGroup([W,N]); + boundaryGroups.ES = multiblock.BoundaryGroup([E,S]); + boundaryGroups.EN = multiblock.BoundaryGroup([E,N]); + + obj.connections = conn; + obj.nBlocks = nBlocks; + obj.boundaryGroups = boundaryGroups; + obj.blockTi = blockTi; + obj.xlims = xlims; + obj.ylims = ylims; + + end + + + % Returns a multiblock.Grid given some parameters + % ms: cell array of [mx, my] vectors + % For same [mx, my] in every block, just input one vector. + function g = getGrid(obj, ms, varargin) + + default_arg('ms',[21,21]) + + % Extend ms if input is a single vector + if (numel(ms) == 2) && ~iscell(ms) + m = ms; + ms = cell(1,obj.nBlocks); + for i = 1:obj.nBlocks + ms{i} = m; + end + end + + grids = cell(1, obj.nBlocks); + for i = 1:obj.nBlocks + grids{i} = grid.equidistant(ms{i}, obj.xlims{i}, obj.ylims{i}); + end + + g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); + end + + % label is the type of label used for plotting, + % default is block name, 'id' show the index for each block. + function show(obj, label, gridLines, varargin) + default_arg('label', 'name') + default_arg('gridLines', false); + + if isempty('label') && ~gridLines + for i = 1:obj.nBlocks + obj.blockTi{i}.show(2,2); + end + axis equal + return + end + + if gridLines + m = 10; + for i = 1:obj.nBlocks + obj.blockTi{i}.show(m,m); + end + end + + + switch label + case 'name' + labels = obj.blockNames; + case 'id' + labels = {}; + for i = 1:obj.nBlocks + labels{i} = num2str(i); + end + otherwise + axis equal + return + end + + for i = 1:obj.nBlocks + parametrization.Ti.label(obj.blockTi{i}, labels{i}); + end + + axis equal + end + + % Returns the grid size of each block in a cell array + % The input parameters are determined by the subclass + function ms = getGridSizes(obj, varargin) + end + end +end
--- a/+multiblock/Def.m Sat Oct 14 22:32:25 2017 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,101 +0,0 @@ -classdef Def - properties - nBlocks - blockMaps % Maps from logical blocks to physical blocks build from transfinite interpolation - blockNames - connections % Cell array specifying connections between blocks - boundaryGroups % Structure of boundaryGroups - end - - methods - % Defines a multiblock setup for transfinite interpolation blocks - % TODO: How to bring in plotting of points? - function obj = Def(blockMaps, connections, boundaryGroups, blockNames) - default_arg('boundaryGroups', struct()); - default_arg('blockNames',{}); - - nBlocks = length(blockMaps); - - obj.nBlocks = nBlocks; - - obj.blockMaps = blockMaps; - - assert(all(size(connections) == [nBlocks, nBlocks])); - obj.connections = connections; - - - if isempty(blockNames) - obj.blockNames = cell(1, nBlocks); - for i = 1:length(blockMaps) - obj.blockNames{i} = sprintf('%d', i); - end - else - assert(length(blockNames) == nBlocks); - obj.blockNames = blockNames; - end - - obj.boundaryGroups = boundaryGroups; - end - - function g = getGrid(obj, varargin) - ms = obj.getGridSizes(varargin{:}); - - grids = cell(1, obj.nBlocks); - for i = 1:obj.nBlocks - grids{i} = grid.equidistantCurvilinear(obj.blockMaps{i}.S, ms{i}); - end - - g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); - end - - function show(obj, label, gridLines, varargin) - default_arg('label', 'name') - default_arg('gridLines', false); - - if isempty('label') && ~gridLines - for i = 1:obj.nBlocks - obj.blockMaps{i}.show(2,2); - end - axis equal - return - end - - if gridLines - ms = obj.getGridSizes(varargin{:}); - for i = 1:obj.nBlocks - obj.blockMaps{i}.show(ms{i}(1),ms{i}(2)); - end - end - - - switch label - case 'name' - labels = obj.blockNames; - case 'id' - labels = {}; - for i = 1:obj.nBlocks - labels{i} = num2str(i); - end - otherwise - axis equal - return - end - - for i = 1:obj.nBlocks - parametrization.Ti.label(obj.blockMaps{i}, labels{i}); - end - - axis equal - end - end - - methods (Abstract) - % Returns the grid size of each block in a cell array - % The input parameters are determined by the subclass - ms = getGridSizes(obj, varargin) - % end - end - -end - -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DefCurvilinear.m Sat Oct 14 22:36:31 2017 -0700 @@ -0,0 +1,101 @@ +classdef DefCurvilinear < multiblock.Definition + properties + nBlocks + blockMaps % Maps from logical blocks to physical blocks build from transfinite interpolation + blockNames + connections % Cell array specifying connections between blocks + boundaryGroups % Structure of boundaryGroups + end + + methods + % Defines a multiblock setup for transfinite interpolation blocks + % TODO: How to bring in plotting of points? + function obj = DefCurvilinear(blockMaps, connections, boundaryGroups, blockNames) + default_arg('boundaryGroups', struct()); + default_arg('blockNames',{}); + + nBlocks = length(blockMaps); + + obj.nBlocks = nBlocks; + + obj.blockMaps = blockMaps; + + assert(all(size(connections) == [nBlocks, nBlocks])); + obj.connections = connections; + + + if isempty(blockNames) + obj.blockNames = cell(1, nBlocks); + for i = 1:length(blockMaps) + obj.blockNames{i} = sprintf('%d', i); + end + else + assert(length(blockNames) == nBlocks); + obj.blockNames = blockNames; + end + + obj.boundaryGroups = boundaryGroups; + end + + function g = getGrid(obj, varargin) + ms = obj.getGridSizes(varargin{:}); + + grids = cell(1, obj.nBlocks); + for i = 1:obj.nBlocks + grids{i} = grid.equidistantCurvilinear(obj.blockMaps{i}.S, ms{i}); + end + + g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); + end + + function show(obj, label, gridLines, varargin) + default_arg('label', 'name') + default_arg('gridLines', false); + + if isempty('label') && ~gridLines + for i = 1:obj.nBlocks + obj.blockMaps{i}.show(2,2); + end + axis equal + return + end + + if gridLines + ms = obj.getGridSizes(varargin{:}); + for i = 1:obj.nBlocks + obj.blockMaps{i}.show(ms{i}(1),ms{i}(2)); + end + end + + + switch label + case 'name' + labels = obj.blockNames; + case 'id' + labels = {}; + for i = 1:obj.nBlocks + labels{i} = num2str(i); + end + otherwise + axis equal + return + end + + for i = 1:obj.nBlocks + parametrization.Ti.label(obj.blockMaps{i}, labels{i}); + end + + axis equal + end + end + + methods (Abstract) + % Returns the grid size of each block in a cell array + % The input parameters are determined by the subclass + ms = getGridSizes(obj, varargin) + % end + end + +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Definition.m Sat Oct 14 22:36:31 2017 -0700 @@ -0,0 +1,11 @@ +classdef Definition + methods (Abstract) + + % Returns a multiblock.Grid given some parameters + g = getGrid(obj, varargin) + + % label is the type of label used for plotting, + % default is block name, 'id' show the index for each block. + show(obj, label, gridLines, varargin) + end +end
--- a/+scheme/Utux.m Sat Oct 14 22:32:25 2017 -0700 +++ b/+scheme/Utux.m Sat Oct 14 22:36:31 2017 -0700 @@ -2,7 +2,7 @@ properties m % Number of points in each direction, possibly a vector h % Grid spacing - x % Grid + grid % Grid order % Order accuracy for the approximation H % Discrete norm @@ -17,41 +17,40 @@ methods - function obj = Utux(m,xlim,order,operator) - default_arg('a',1); - - %Old operators - % [x, h] = util.get_grid(xlim{:},m); - %ops = sbp.Ordinary(m,h,order); - - + 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); - obj.D1 = ops.D1; +% case 'NonEquidistant' +% ops = sbp.D1Nonequidistant(m,xlim,order); +% obj.D1 = ops.D1; case 'Standard' ops = sbp.D2Standard(m,xlim,order); obj.D1 = ops.D1; - case 'Upwind' - ops = sbp.D1Upwind(m,xlim,order); - obj.D1 = ops.Dm; +% case 'Upwind' +% ops = sbp.D1Upwind(m,xlim,order); +% obj.D1 = ops.Dm; otherwise error('Unvalid operator') end - obj.x=ops.x; + + 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; + obj.D = -obj.D1; obj.m = m; obj.h = ops.h; obj.order = order; - obj.x = ops.x; end % Closure functions return the opertors applied to the own doamin to close the boundary @@ -61,17 +60,27 @@ % data is a function returning the data that should be applied at the boundary. % neighbour_scheme is an instance of Scheme that should be interfaced to. % neighbour_boundary is a string specifying which boundary to interface to. - function [closure, penalty] = boundary_condition(obj,boundary,type,data) - default_arg('type','neumann'); - default_arg('data',0); + function [closure, penalty] = boundary_condition(obj,boundary,type) + default_arg('type','dirichlet'); tau =-1*obj.e_l; closure = obj.Hi*tau*obj.e_l'; - penalty = 0*obj.e_l; + penalty = -obj.Hi*tau; end function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - error('An interface function does not exist yet'); + switch boundary + % Upwind coupling + case {'l','left'} + tau = -1*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'; + penalty = -obj.Hi*tau*neighbour_scheme.e_l'; + end + end function N = size(obj) @@ -81,9 +90,9 @@ end methods(Static) - % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u + % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u % and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Utux2D.m Sat Oct 14 22:36:31 2017 -0700 @@ -0,0 +1,177 @@ +classdef Utux2D < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + grid % Grid + order % Order accuracy for the approximation + v0 % Initial data + + a % Wave speed a = [a1, a2]; + + H % Discrete norm + H_x, H_y % Norms in the x and y directions + Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms + + % 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 + + + end + + + methods + function obj = Utux2D(g ,order, opSet, a, coupling_type) + + default_arg('coupling_type','upwind'); + default_arg('a',1/sqrt(2)*[1, 1]); + default_arg('opSet',@sbp.D2Standard); + assert(isa(g, 'grid.Cartesian')) + + m = g.size(); + m_x = m(1); + m_y = m(2); + m_tot = g.N(); + + xlim = {g.x{1}(1), g.x{1}(end)}; + ylim = {g.x{2}(1), g.x{2}(end)}; + obj.grid = g; + + % Operator sets + ops_x = opSet(m_x, xlim, order); + 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); + obj.Hi = kron(Hxi,Hyi); + obj.Hx = kron(Hx,Iy); + 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); + obj.e_s = kr(Ix, ops_y.e_l); + obj.e_n = kr(Ix, ops_y.e_r); + + obj.m = m; + obj.h = [ops_x.h ops_y.h]; + obj.order = order; + obj.a = a; + obj.coupling_type = coupling_type; + obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy); + + end + % Closure functions return the opertors applied to the own domain 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'. + % type is a string specifying the type of boundary condition if there are several. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % 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 + penalty = -obj.Hi*tau; + + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + % Get neighbour boundary operator + switch neighbour_boundary + case {'e','E','east','East'} + e_neighbour = neighbour_scheme.e_e; + case {'w','W','west','West'} + e_neighbour = neighbour_scheme.e_w; + case {'n','N','north','North'} + e_neighbour = neighbour_scheme.e_n; + case {'s','S','south','South'} + e_neighbour = neighbour_scheme.e_s; + end + + switch obj.coupling_type + + % Upwind coupling (energy dissipation) + case 'upwind' + sigma_ds = -1; %"Downstream" penalty + sigma_us = 0; %"Upstream" penalty + + % Energy-preserving coupling (no energy dissipation) + case 'centered' + sigma_ds = -1/2; %"Downstream" penalty + sigma_us = 1/2; %"Upstream" penalty + + otherwise + error(['Interface coupling type ' coupling_type ' is not available.']) + end + + switch boundary + case {'w','W','west','West'} + tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y; + closure = obj.Hi*tau*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'; + case {'s','S','south','South'} + tau = sigma_ds*obj.a(2)*obj.e_s*obj.H_x; + closure = obj.Hi*tau*obj.e_s'; + case {'n','N','north','North'} + tau = sigma_us*obj.a(2)*obj.e_n*obj.H_x; + closure = obj.Hi*tau*obj.e_n'; + end + penalty = -obj.Hi*tau*e_neighbour'; + + end + + function N = size(obj) + N = obj.m; + end + + end + + methods(Static) + % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') + function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) + [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); + [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); + end + end +end \ No newline at end of file