-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathstarmap_2d_ex_linesource.m
137 lines (127 loc) · 5.52 KB
/
starmap_2d_ex_linesource.m
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
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
function starmap_2d_ex_linesource
%STARMAP_2D_EX_LINESOURCE
% Example case for STARMAP_SOLVER, a second order staggered
% grid finite difference solver for linear hyperbolic moment
% approximations to radiative transfer in 1D-3D geometry.
%
% Example: Line source benchmark test case, as defined in
% [Ganapol, JCP 210 (2005) 386-399].
% This test case is adapted to using a narrow Gaussian as
% initial condition.
%
% Remarks:
% 1) Since the numerical scheme does not possess slope
% limiters, the initial Gaussian must not be chosen too
% narrow (or else oscillations will arise).
% 2) For this problem, the SPN approximation is equivalent
% to the PN approximation, see
% [McClarren, Trans. Theory Stat. Phys. 39 (2011) 73-109],
% [Olbrant, Larsen, Frank, Seibold, JCP 238 (2013) 315-336].
% Since SPN contains much fewer moments than PN, it is
% natural to use here.
%
% Version 2.0
% Copyright (c) 06/28/2022 Benjamin Seibold, Martin Frank, and
% Rujeko Chinomona
% http://www.math.temple.edu/~seibold
% https://www.scc.kit.edu/personen/martin.frank.php
% https://rujekoc.github.io/
%
% Contributers: Edgar Olbrant (v1.0), Kerstin Kuepper (v1.5,v2.0).
%
% StaRMAP project website:
% https://github.com/starmap-project
% For license, see files LICENSE.txt or starmap_solver.m, as published on
% https://github.com/starmap-project/starmap
%========================================================================
% Problem Parameters
%========================================================================
prob = struct(...
'name','Line Source Test',... % name of example
'closure','SP',... % type of closure (can be 'P' or 'SP')
'n_mom',39,... % order of moment approximation
'sigma_a',@sigma_a,... % absorption coefficient (defined below)
'sigma_s0',@sigma_s0,... % isotropic scattering coefficient (def. below)
'ic',@initial,... % initial conditions (defined below)
'ax',[-1 1 -1 1]*0.6,... % coordinates of computational domain
'n',[1 1]*150,... % numbers of grid cells in each coordinate direction
'bc',[1 1],... % type of boundary cond. (0 = periodic, 1 = extrapolation)
't_plot',linspace(0,0.5,21),... % output times
'output',@output... % problem-specific output routine (defined below)
);
%========================================================================
% Moment System Setup and Solver Execution
%========================================================================
par = starmap_init(prob); % Configure data structures for starmap solver
starmap_solver(par) % Run solver.
%========================================================================
% Problem Specific Functions
%========================================================================
function f = initial(x,y,k)
% Initial conditions (for (k-1)-st moment).
t = 3.2e-4; % Pseudo-time for smoothing initial Dirac delta.
f = 1/(4*pi*t)*exp(-(x.^2+y.^2)/(4*t))*(k==1);
function f = sigma_a(x,y)
% Absorption coefficient.
f = 0;
function f = sigma_s0(x,y)
% Isotropic scattering coefficient.
f = 1;
function output(par,x,y,U,step)
% Plotting routine, showing the 2D results and a 1D cross-section.
t = par.t_plot(step); % Time of output.
[r0,phi0] = line_source_solution(t); % Compute true solution.
phi_max = max(phi0(1)+eps,max(max(U))/4); % Axis/color scaling.
clf, subplot(1,3,1:2), imagesc(x,y,U') % 2D plot of approximation.
axis xy equal tight, caxis([-1 1]*phi_max*1.6)
title(sprintf('%s with %s%d at t = %0.2f',par.name,par.closure,par.n_mom,t))
colormap jet(255); colorbar
subplot(1,3,3)
r = linspace(0,max(x),100); phi = interp2(x,y,U,r,0); % 1D cross section.
plot(r0,phi0,'-','linewidth',2,'color',[.6 .6 1]) % Plot true solution.
hold on, plot(r,phi,'r-','linewidth',1), hold off % Plot approximation.
axis([0 max(x) [-0.75 2.5]*phi_max])
legend('true solution',sprintf('%s%d model',par.closure,par.n_mom), ...
'Location','SouthEast')
title('Cut at x>0, y=0')
drawnow
%========================================================================
% Reference solution of line source problem
%========================================================================
function [rho,f] = line_source_solution(t)
% Evaluates the true solution of the line source test case,
% using Ganapol's quadrature formulas.
rho = [(1-linspace(1,0,25).^2)*t 1];
f = phi_l(rho,t);
function f = phi_l(rho,t)
eta = rho/t;
ind = find(eta<1);
f = rho*0; phi_l0 = f;
for k = ind % Compute integral.
f(k) = quad(@(w) phi_pt(t*sqrt(eta(k).^2+w.^2),t),...
0,sqrt(1-eta(k).^2),1e-5);
end
phi_l0(ind) = exp(-t)/(2*pi*t^2)./sqrt(1-eta(ind).^2);
f = phi_l0+(2*t)*f;
function f = phi_pt(r,t)
r = max(r,1e-10); % Move zero value into small positive regime.
eta = r/t;
ind = find(eta<1);
g = r*0;
for k = ind % Compute integral.
g(k) = quad(@(u) sec(u/2).^2.*real(...
(eta(k)+1i*tan(u/2)).*...
xi(eta(k),u).^3.*exp(t/2*(1-eta(k).^2).*xi(eta(k),u)) ...
),0,pi,1e-2);
end
f = 1/(2*pi)*exp(-t)./(4*pi*r.*t.^2)*(t/2)^2.*(1-eta.^2).*g; % Assemble
f(ind) = f(ind)+phi_pt1(r(ind),t); % final result.
function f = xi(eta,u)
q = (1+eta)./(1-eta);
q = min(max(q,-1e80),1e80); % Remove singularity at eta = 1.
f = (log(q)+1i*u)./(eta+1i*tan(u/2));
function y = phi_pt1(r,t)
eta = r/t;
q = (1+eta)./(1-eta);
q = min(max(q,-1e80),1e80); % Remove singularity at eta = 1.
y = exp(-t)./(4*pi*r*t^2)*t.*log(q);