MATLAB-Transport-Pack/test_assembly.m at master · robertsj/MATLAB-Transport-Pack · 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
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
%> @file  test_assembly.m
%> @brief Tests the eigenvalue solver on a 2D assembly problem.
%
% That's about 2 minutes per sweep, wish seems reasonable enough.  However,
% if we could limit that to maybe 10 sweeps via CMFD, we're in
% business.  Note, the GMRES

clear class

% ==============================================================================
% INPUT
% ==============================================================================
input = Input();
put(input, 'number_groups',         2);
put(input, 'number_groups',         2);
put(input, 'eigen_tolerance',       1e-4);
put(input, 'eigen_max_iters',       100);
put(input, 'inner_tolerance',       0.0001);
put(input, 'inner_solver',          'SI');
put(input, 'livolant_free_iters',   3);
put(input, 'livolant_accel_iters',  3);
put(input, 'bc_left',               'reflect');
put(input, 'bc_right',              'reflect');
put(input, 'bc_top',                'reflect');
put(input, 'bc_bottom',             'reflect');
input.number_groups = 2;


% ==============================================================================
% MATERIALS (Test two group data)
% ==============================================================================
mat = test_materials(2);

% ==============================================================================
% MESH (Assembly of pins)
% ==============================================================================

% Shared pin properties
pitch   = 1.26;
radii   = 0.54;
number  = 12;

% Pin 1
matid  = [2 1]; % IN to OUT
pin1   = PinCell(pitch, radii, matid);
meshify(pin1, number);
figure(1)
plot_mesh(pin1);

% Pin 2
matid  = [4 1];
pin2   = PinCell(pitch, radii, matid);
meshify(pin2, number);
figure(2)
plot_mesh(pin2);

% Pin 3
matid  = [1];
pin3   = PinCell(pitch, [], matid);
meshify(pin3, number);
figure(3)
plot_mesh(pin3);

% Make the assembly.
pin_map = [ 1 1 1 1 1 1 1 1 1
            1 1 1 2 1 1 1 1 1
            1 1 1 1 1 1 2 1 1
            1 2 1 3 3 1 1 1 1
            1 1 1 3 3 3 1 1 1
            1 1 1 1 3 3 1 2 1
            1 1 2 1 1 1 1 1 1
            1 1 1 1 1 2 1 1 1
            1 1 1 1 1 1 1 1 1 ];
mesh = Assembly({pin1, pin2, pin3}, pin_map);
meshify(mesh);
figure(4)
plot_mesh_map(mesh, 'MATERIAL')


% ==============================================================================
% SETUP
% ==============================================================================
state       = State(input, mesh);
quadrature  = QuadrupleRange(8);
boundary    = BoundaryMesh(input, mesh, quadrature);
q_e         = Source(mesh, 2);                  % Not initialized = not used.
q_f         = FissionSource(state, mesh, mat);  % Inititalized = used.
initialize(q_f);

% ==============================================================================
% SOLVE
% ==============================================================================
solver = Eigensolver(input,         ...
                     state,         ...
                     boundary,      ...
                     mesh,          ...
                     mat,           ...
                     quadrature,	...
                     q_e,           ...
                     q_f);

tic
out = solve(solver);
toc

% ==============================================================================
% POSTPROCESS
% ==============================================================================
VTK.savevtk(state, 2, mesh, 'nine_by_nine_pin.vtk')

save('one_pin.mat')

%figure(5)
%subplot(2, 1, 1)
%f1 = flux(state, 1);
%plot_flux(mesh, f1)
%axis equal
%subplot(2, 1, 2)
%f2 = flux(state, 2);
%plot_flux(mesh, f2)
%axis equal