kernelpack-matlab/examples/diffusion_solver_example.m at main · ShankarLab/kernelpack-matlab · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
function diffusion_solver_example()
%DIFFUSION_SOLVER_EXAMPLE Advance a fixed-domain diffusion problem with BDF steps.

t = linspace(0, 2*pi, 80).';
t(end) = [];
curve = [cos(t), 0.8 * sin(t)];

% Build a smooth closed boundary and its implicit model.
surface = kp.geometry.EmbeddedSurface();
surface.setDataSites(curve);
surface.buildClosedGeometricModelPS(2, 0.06, size(curve, 1));
surface.buildLevelSetFromGeometricModel([]);

% Generate the geometry-backed domain used by the diffusion stepper.
generator = kp.nodes.DomainNodeGenerator();
domain = generator.buildDomainDescriptorFromGeometry(surface, 0.1, ...
    'Seed', 23, ...
    'StripCount', 5, ...
    'DoOuterRefinement', true, ...
    'OuterFractionOfh', 0.5, ...
    'OuterRefinementZoneSizeAsMultipleOfh', 2.0);

% Initialize the diffusion solver once, then step it with a BDF update.
solver = kp.solvers.DiffusionSolver( ...
    'LapAssembler', 'fd', ...
    'BCAssembler', 'fd', ...
    'LapStencil', 'wls', ...
    'BCStencil', 'wls');
nu = 0.25;
dt = 0.02;
solver.init(domain, 3, dt, nu);

uExact = @(time, X) exp(-time) .* (X(:, 1).^2 + X(:, 2).^2);
forcing = @(nuValue, time, X) -exp(-time) .* (X(:, 1).^2 + X(:, 2).^2) - 4 * nuValue * exp(-time);
neuCoeff = @(Xb) zeros(size(Xb, 1), 1);
dirCoeff = @(Xb) ones(size(Xb, 1), 1);
bc = @(NeuCoeffs, DirCoeffs, nr, time, Xb) uExact(time, Xb); %#ok<INUSD>

solver.setInitialState(uExact(0, domain.getIntBdryNodes()));
u1 = solver.bdf1Step(dt, forcing, neuCoeff, dirCoeff, bc);
u2 = solver.bdf2Step(2 * dt, forcing, neuCoeff, dirCoeff, bc);
u3 = solver.bdf3Step(3 * dt, forcing, neuCoeff, dirCoeff, bc);

Xphys = domain.getIntBdryNodes();

figure('Color', 'w');
tiledlayout(1, 3);

nexttile;
scatter(Xphys(:, 1), Xphys(:, 2), 22, u1, 'filled');
axis equal;
grid on;
title('BDF1 state');
colorbar;

nexttile;
scatter(Xphys(:, 1), Xphys(:, 2), 22, u2, 'filled');
axis equal;
grid on;
title('BDF2 state');
colorbar;

nexttile;
scatter(Xphys(:, 1), Xphys(:, 2), 22, abs(u3 - uExact(3 * dt, Xphys)), 'filled');
axis equal;
grid on;
title('BDF3 absolute error');
colorbar;
end