Mercurial > repos > public > sbplib
changeset 1251:6424745e1b58 feature/volcano
Merge in latest changes from default
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Wed, 20 Nov 2019 14:26:57 -0800 |
| parents | 10881b234f77 (current diff) 8ec777fb473e (diff) |
| children | |
| files | |
| diffstat | 5 files changed, 757 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
diff -r 10881b234f77 -r 6424745e1b58 +multiblock/joinGrids.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/joinGrids.m Wed Nov 20 14:26:57 2019 -0800 @@ -0,0 +1,90 @@ +% Connects several multiblock grids into one grid +% gs -- Cell array of multiblock grids +% gConnections -- Upper-triangular cell matrix +% gConnections{i,j} specifies all connections between grid i and grid j +% Example: +% gConnections{i,j} = { {{1,'e'},{2,'w'}}, {{5,'s'},{2,'n'}},... }; +% names -- (Optional) cell array of strings, used for boundary groups in new grid +% default: names = {'g1', 'g2', ..., 'gN'}; +% Boundary groups from grid i are contained in g.boundaryGroups.(names{i}), etc. +function g = joinGrids(gs, gConnections, names) + + nGrids = numel(gs); + + % Default names are {'g1', 'g2', ... 'gN'}. + defaultNames = cell(nGrids, 1); + for i = 1:nGrids + defaultNames{i} = sprintf('g%d',i); + end + default_arg('names', defaultNames); + + nBlocks = 0; + for i = 1:nGrids + nBlocks = nBlocks + gs{i}.nBlocks(); + end + + % Create vector of cumulative sum of number of blocks per grid + startIndex = zeros(1, nGrids); + for i = 2:nGrids + startIndex(i) = startIndex(i-1) + gs{i-1}.nBlocks(); + end + + % Create cell array of all grids + grids = cell(nBlocks, 1); + for i = 1:nGrids + for j = 1:gs{i}.nBlocks(); + grids{startIndex(i)+j} = gs{i}.grids{j}; + end + end + + % Create cell matrix of connections + connections = cell(nBlocks, nBlocks); + + % Connections within grids + for i = 1:nGrids + for j = 1:gs{i}.nBlocks() + for k = 1:gs{i}.nBlocks() + connections{startIndex(i)+j,startIndex(i)+k} = gs{i}.connections{j,k}; + end + end + end + + % Connections between grids + for i = 1:nGrids + for j = 1:nGrids + for k = 1:numel(gConnections{i,j}) + b1 = gConnections{i,j}{k}{1}; + id1 = b1{1}; + str1 = b1{2}; + + b2 = gConnections{i,j}{k}{2}; + id2 = b2{1}; + str2 = b2{2}; + + connections{startIndex(i)+id1, startIndex(j)+id2} = {str1, str2}; + end + end + end + + % Boundary groups + boundaryGroups = struct; + for i = 1:nGrids + bgs = gs{i}.boundaryGroups; + bgNames = fieldnames(bgs); + for j = 1:numel(bgNames) + bg = bgs.(bgNames{j}); + + % Shift block id:s in boundary groups + for k = 1:length(bg) + bg{k}{1} = bg{k}{1} + startIndex(i); + end + + bgs.(bgNames{j}) = bg; + end + boundaryGroups.(names{i}) = bgs; + end + + % Create grid object + g = multiblock.Grid(grids, connections, boundaryGroups); + +end \ No newline at end of file
diff -r 10881b234f77 -r 6424745e1b58 +sbp/D2VariableCompatible.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/D2VariableCompatible.m Wed Nov 20 14:26:57 2019 -0800 @@ -0,0 +1,83 @@ +classdef D2VariableCompatible < sbp.OpSet + properties + D1 % SBP operator approximating first derivative + H % Norm matrix + HI % H^-1 + Q % Skew-symmetric matrix + e_l % Left boundary operator + e_r % Right boundary operator + D2 % SBP operator for second derivative + M % Norm matrix, second derivative + d1_l % Left boundary first derivative + d1_r % Right boundary first derivative + m % Number of grid points. + h % Step size + x % grid + borrowing % Struct with borrowing limits for different norm matrices + end + + methods + function obj = D2VariableCompatible(m,lim,order) + + x_l = lim{1}; + x_r = lim{2}; + L = x_r-x_l; + obj.h = L/(m-1); + obj.x = linspace(x_l,x_r,m)'; + + switch order + + case 6 + + [obj.H, obj.HI, obj.D1, D2, ... + ~, obj.e_l, obj.e_r, ~, ~, ~, ~, ~,... + d1_l, d1_r] = ... + sbp.implementations.d4_variable_6(m, obj.h); + + case 4 + [obj.H, obj.HI, obj.D1, D2, obj.e_l,... + obj.e_r, d1_l, d1_r] = ... + sbp.implementations.d2_variable_4(m,obj.h); + case 2 + [obj.H, obj.HI, obj.D1, D2, obj.e_l,... + obj.e_r, d1_l, d1_r] = ... + sbp.implementations.d2_variable_2(m,obj.h); + + otherwise + error('Invalid operator order %d.',order); + end + obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h, + obj.borrowing.M.d1 = obj.H(1,1)/obj.h; % First element in H/h is borrowing also for M + obj.borrowing.R.delta_D = inf; + obj.m = m; + obj.M = []; + + + D1 = obj.D1; + e_r = obj.e_r; + e_l = obj.e_l; + + % D2 = Hinv * (-M + br*er*d1r^T - bl*el*d1l^T); + % Replace d1' by e'*D1 in D2. + correction_l = obj.HI*(e_l*d1_l' - e_l*e_l'*D1); + correction_r = - obj.HI*(e_r*d1_r' - e_r*e_r'*D1); + + D2_compatible = @(b) D2(b) + b(1)*correction_l + b(m)*correction_r; + + obj.D2 = D2_compatible; + obj.d1_l = (e_l'*D1)'; + obj.d1_r = (e_r'*D1)'; + + end + function str = string(obj) + str = [class(obj) '_' num2str(obj.order)]; + end + end + + +end + + + + +
diff -r 10881b234f77 -r 6424745e1b58 diracDiscr.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diracDiscr.m Wed Nov 20 14:26:57 2019 -0800 @@ -0,0 +1,106 @@ +% n-dimensional delta function +% g: cartesian grid +% x_s: source point coordinate vector, e.g. [x; y] or [x; y; z]. +% m_order: Number of moment conditions +% s_order: Number of smoothness conditions +% H: cell array of 1D norm matrices +function d = diracDiscr(g, x_s, m_order, s_order, H) + assertType(g, 'grid.Cartesian'); + + dim = g.d; + d_1D = cell(dim,1); + + % Allow for non-cell input in 1D + if dim == 1 + H = {H}; + end + % Create 1D dirac discr for each coordinate dir. + for i = 1:dim + d_1D{i} = diracDiscr1D(x_s(i), g.x{i}, m_order, s_order, H{i}); + end + + d = d_1D{dim}; + for i = dim-1: -1: 1 + % Perform outer product, transpose, and then turn into column vector + d = (d_1D{i}*d')'; + d = d(:); + end +end + + +% Helper function for 1D delta functions +function ret = diracDiscr1D(x_s, x, m_order, s_order, H) + + % Return zeros if x0 is outside grid + if x_s < x(1) || x_s > x(end) + ret = zeros(size(x)); + return + end + + tot_order = m_order+s_order; %This is equiv. to the number of equations solved for + S = []; + M = []; + + % Get interior grid spacing + middle = floor(length(x)/2); + h = x(middle+1) - x(middle); % Use middle point to allow for staggered grids. + + index = sourceIndices(x_s, x, tot_order, h); + + polynomial = (x(index)-x(index(1)))/(x(index(end))-x(index(1))); + x_0 = (x_s-x(index(1)))/(x(index(end))-x(index(1))); + + quadrature = diag(H); + quadrature_weights = quadrature(index)/h; + + h_polynomial = polynomial(2)-polynomial(1); + b = zeros(tot_order,1); + + for i = 1:m_order + b(i,1) = x_0^(i-1); + end + + for i = 1:tot_order + for j = 1:m_order + M(j,i) = polynomial(i)^(j-1)*h_polynomial*quadrature_weights(i); + end + end + + for i = 1:tot_order + for j = 1:s_order + S(j,i) = (-1)^(i-1)*polynomial(i)^(j-1); + end + end + + A = [M;S]; + + d = A\b; + ret = x*0; + ret(index) = d/h*h_polynomial; +end + + +function I = sourceIndices(x_s, x, tot_order, h) + % Find the indices that are within range of of the point source location + I = find(tot_order*h/2 >= abs(x-x_s)); + + if length(I) > tot_order + if length(I) == tot_order + 2 + I = I(2:end-1); + elseif length(I) == tot_order + 1 + I = I(1:end-1); + end + elseif length(I) < tot_order + if x_s < x(1) + ceil(tot_order/2)*h + I = 1:tot_order; + elseif x_s > x(end) - ceil(tot_order/2)*h + I = length(x)-tot_order+1:length(x); + else + if I(end) < length(x) + I = [I; I(end)+1]; + else + I = [I(1)-1; I]; + end + end + end +end
diff -r 10881b234f77 -r 6424745e1b58 diracDiscrCurve.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diracDiscrCurve.m Wed Nov 20 14:26:57 2019 -0800 @@ -0,0 +1,63 @@ +% 2-dimensional delta function for single-block curvilinear grid +% x_s: source point coordinate vector, e.g. [x; y] or [x; y; z]. +% g: single-block grid containing the source +% m_order: Number of moment conditions +% s_order: Number of smoothness conditions +% order: Order of SBP derivative approximations +% opSet: Cell array of function handle to opSet generator +function d = diracDiscrCurve(x_s, g, m_order, s_order, order, opSet) + default_arg('order', m_order); + default_arg('opSet', {@sbp.D2Variable, @sbp.D2Variable}); + + dim = length(x_s); + assert(dim == 2, 'diracDiscrCurve: Only implemented for 2d.'); + assertType(g, 'grid.Curvilinear'); + + % Compute Jacobian + J = jacobian(g, opSet, order); + + % Find approximate logical coordinates of point source + X = g.points(); + U = g.logic.points(); + U_interp = scatteredInterpolant(X, U(:,1)); + V_interp = scatteredInterpolant(X, U(:,2)); + uS = U_interp(x_s); + vS = V_interp(x_s); + + % Get quadrature matrices for moment conditions + m = g.size(); + ops_u = opSet{1}(m(1), {0, 1}, order); + ops_v = opSet{2}(m(2), {0, 1}, order); + H_u = ops_u.H; + H_v = ops_v.H; + + % Get delta function for logical grid and scale by Jacobian + d = (1./J).*diracDiscr(g, [uS; vS], m_order, s_order, {H_u, H_v}); +end + +function J = jacobian(g, opSet, order) + m = g.size(); + m_u = m(1); + m_v = m(2); + ops_u = opSet{1}(m_u, {0, 1}, order); + ops_v = opSet{2}(m_v, {0, 1}, order); + I_u = speye(m_u); + I_v = speye(m_v); + + D1_u = ops_u.D1; + D1_v = ops_v.D1; + + Du = kr(D1_u,I_v); + Dv = kr(I_u,D1_v); + + coords = g.points(); + x = coords(:,1); + y = coords(:,2); + + x_u = Du*x; + x_v = Dv*x; + y_u = Du*y; + y_v = Dv*y; + + J = x_u.*y_v - x_v.*y_u; +end
diff -r 10881b234f77 -r 6424745e1b58 diracDiscrTest.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diracDiscrTest.m Wed Nov 20 14:26:57 2019 -0800 @@ -0,0 +1,415 @@ +function tests = diracDiscrTest() + tests = functiontests(localfunctions); +end + +%TODO: +% 1. Test discretizing with smoothness conditions. +% Only discretization with moment conditions currently tested. +% 2. Test using other types of grids. Only equidistant grids currently used. + +function testLeftRandom(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + rng(1) % Set seed + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [g, H, fs] = setup1D(order, mom_cond); + xl = g.lim{1}{1}; + h = g.scaling(); + x = g.x{1}; + + % Test random points near left boundary + x0s = xl + 2*h*rand(1,10); + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(g, x0, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testRightRandom(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + rng(1) % Set seed + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [g, H, fs] = setup1D(order, mom_cond); + xr = g.lim{1}{2}; + h = g.scaling(); + x = g.x{1}; + + % Test random points near right boundary + x0s = xr - 2*h*rand(1,10); + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(g, x0, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testInteriorRandom(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + rng(1) % Set seed + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [g, H, fs] = setup1D(order, mom_cond); + xl = g.lim{1}{1}; + xr = g.lim{1}{2}; + h = g.scaling(); + x = g.x{1}; + + % Test random points in interior + x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20); + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(g, x0, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +% x0 outside grid should yield 0 integral! +function testX0OutsideGrid(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [g, H, fs] = setup1D(order, mom_cond); + xl = g.lim{1}{1}; + xr = g.lim{1}{2}; + h = g.scaling(); + x = g.x{1}; + + % Test points outisde grid + x0s = [xl-1.1*h, xr+1.1*h]; + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(g, x0, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - 0); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testAllGP(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [g, H, fs] = setup1D(order, mom_cond); + x = g.x{1}; + + % Test all grid points + x0s = x; + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(g, x0, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testHalfGP(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [g, H, fs] = setup1D(order, mom_cond); + x = g.x{1}; + + % Test halfway between all grid points + x0s = 1/2*(x(2:end)+x(1:end-1)); + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(g, x0, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +%=============== 2D tests ============================== +function testAllGP2D(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [g, H, fs] = setup2D(order, mom_cond); + H_global = kron(H{1}, H{2}); + X = g.points(); + % Test all grid points + x0s = X; + + for j = 1:length(fs) + f = fs{j}; + fx = f(X(:,1), X(:,2)); + for i = 1:length(x0s) + x0 = x0s(i,:); + delta = diracDiscr(g, x0, mom_cond, 0, H); + integral = delta'*H_global*fx; + err = abs(integral - f(x0(1), x0(2))); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testAllRandom2D(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + rng(1) % Set seed + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [g, H, fs] = setup2D(order, mom_cond); + H_global = kron(H{1}, H{2}); + X = g.points(); + xl = g.lim{1}{1}; + xr = g.lim{1}{2}; + yl = g.lim{2}{1}; + yr = g.lim{2}{2}; + h = g.scaling(); + + + % Test random points, even outside grid + Npoints = 100; + x0s = [(xl-3*h(1)) + (xr-xl+6*h(1))*rand(Npoints,1), ... + (yl-3*h(2)) + (yr-yl+6*h(2))*rand(Npoints,1) ]; + + for j = 1:length(fs) + f = fs{j}; + fx = f(X(:,1), X(:,2)); + for i = 1:length(x0s) + x0 = x0s(i,:); + delta = diracDiscr(g, x0, mom_cond, 0, H); + integral = delta'*H_global*fx; + + % Integral should be 0 if point is outside grid + if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr + err = abs(integral - 0); + else + err = abs(integral - f(x0(1), x0(2))); + end + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +%=============== 3D tests ============================== +function testAllGP3D(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [g, H, fs] = setup3D(order, mom_cond); + H_global = kron(kron(H{1}, H{2}), H{3}); + X = g.points(); + % Test all grid points + x0s = X; + + for j = 1:length(fs) + f = fs{j}; + fx = f(X(:,1), X(:,2), X(:,3)); + for i = 1:length(x0s) + x0 = x0s(i,:); + delta = diracDiscr(g, x0, mom_cond, 0, H); + integral = delta'*H_global*fx; + err = abs(integral - f(x0(1), x0(2), x0(3))); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testAllRandom3D(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + rng(1) % Set seed + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [g, H, fs] = setup3D(order, mom_cond); + H_global = kron(kron(H{1}, H{2}), H{3}); + X = g.points(); + xl = g.lim{1}{1}; + xr = g.lim{1}{2}; + yl = g.lim{2}{1}; + yr = g.lim{2}{2}; + zl = g.lim{3}{1}; + zr = g.lim{3}{2}; + h = g.scaling(); + + % Test random points, even outside grid + Npoints = 200; + x0s = [(xl-3*h(1)) + (xr-xl+6*h(1))*rand(Npoints,1), ... + (yl-3*h(2)) + (yr-yl+6*h(2))*rand(Npoints,1), ... + (zl-3*h(3)) + (zr-zl+6*h(3))*rand(Npoints,1) ]; + + for j = 1:length(fs) + f = fs{j}; + fx = f(X(:,1), X(:,2), X(:,3)); + for i = 1:length(x0s) + x0 = x0s(i,:); + delta = diracDiscr(g, x0, mom_cond, 0, H); + integral = delta'*H_global*fx; + + % Integral should be 0 if point is outside grid + if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr || x0(3) < zl || x0(3) > zr + err = abs(integral - 0); + else + err = abs(integral - f(x0(1), x0(2), x0(3))); + end + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + + +% ====================================================== +% ============== Setup functions ======================= +% ====================================================== +function [g, H, fs] = setup1D(order, mom_cond) + + % Grid + xl = -3; + xr = 900; + L = xr-xl; + m = 101; + g = grid.equidistant(m, {xl, xr}); + + % Quadrature + ops = sbp.D2Standard(m, {xl, xr}, order); + H = ops.H; + + % Moment conditions + fs = cell(mom_cond,1); + for p = 0:mom_cond-1 + fs{p+1} = @(x) (x/L).^p; + end + +end + +function [g, H, fs] = setup2D(order, mom_cond) + + % Grid + xlims = {-3, 20}; + ylims = {-11,5}; + Lx = xlims{2} - xlims{1}; + Ly = ylims{2} - ylims{1}; + m = [15, 16]; + g = grid.equidistant(m, xlims, ylims); + + % Quadrature + opsx = sbp.D2Standard(m(1), xlims, order); + opsy = sbp.D2Standard(m(2), ylims, order); + Hx = opsx.H; + Hy = opsy.H; + H = {Hx, Hy}; + + % Moment conditions + fs = cell(mom_cond,1); + for p = 0:mom_cond-1 + fs{p+1} = @(x,y) (x/Lx + y/Ly).^p; + end +end + +function [g, H, fs] = setup3D(order, mom_cond) + + % Grid + xlims = {-3, 20}; + ylims = {-11,5}; + zlims = {2,4}; + Lx = xlims{2} - xlims{1}; + Ly = ylims{2} - ylims{1}; + Lz = zlims{2} - zlims{1}; + m = [13, 14, 15]; + g = grid.equidistant(m, xlims, ylims, zlims); + + % Quadrature + opsx = sbp.D2Standard(m(1), xlims, order); + opsy = sbp.D2Standard(m(2), ylims, order); + opsz = sbp.D2Standard(m(3), zlims, order); + Hx = opsx.H; + Hy = opsy.H; + Hz = opsz.H; + H = {Hx, Hy, Hz}; + + % Moment conditions + fs = cell(mom_cond,1); + for p = 0:mom_cond-1 + fs{p+1} = @(x,y,z) (x/Lx + y/Ly + z/Lz).^p; + end +end