Compare commits

..
40 changed files with 3616 additions and 3323 deletions
-15
View File
@@ -53,16 +53,6 @@ Discretization improvements
ComplexHypreParMatrix::GetSystemMatrix, which previously assumed equal
trial and test spaces.
- Added FiniteElementSpace::GetBoundaryLoopEdgeDofs to extract the edge DOFs on
the perimeter loop of a set of boundary elements, with a ParFiniteElementSpace
overload that reconciles the selection across processor boundaries so the
result is partition invariant. This is useful for imposing boundary conditions
on boundary edge DOFs.
- Added a MaxAbs reduction to GroupCommunicator that selects the signed value of
largest magnitude across a group, keeping its sign. Equal-magnitude ties
resolve deterministically to the positive value.
Meshing improvements
--------------------
- Added support for nonuniform anisotropic mesh refinement on parallel quad/hex
@@ -152,11 +142,6 @@ Miscellaneous
using the new method ApplyDofSigns() in class ParFiniteElementSpace: the
method will return immediately if no sign flips are needed.
- Added support for coefficient-weighted LOR transfer in
L2ProjectionGridTransfer. The transfer conserves the weighted mass, for
example when transferring velocity while conserving density-weighted momentum.
This is illustrated in the lor-transfer and plor-transfer miniapps.
- Added support for saving DataCollection output on the node-local storage,
instead of requiring that the filesystem is shared among all the ranks.
@@ -962,7 +962,7 @@ function(mfem_export_mk_files)
set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} -I${dir}")
endforeach()
# TODO: MFEM_TPLFLAGS: add other TPL flags, in addition to the -I flags.
set(MFEM_INCFLAGS "-I\$(MFEM_INC_DIR) \$(MFEM_TPLFLAGS)")
set(MFEM_INCFLAGS "-I\$(MFEM_INC_DIR) -I\$(MFEM_SOURCE_DIR) \$(MFEM_TPLFLAGS)")
set(MFEM_PICFLAG "")
if (BUILD_SHARED_LIBS)
set(MFEM_PICFLAG "${CMAKE_SHARED_LIBRARY_CXX_FLAGS}")
@@ -1098,6 +1098,7 @@ function(mfem_export_mk_files)
"${INSTALL_LIB_DIR}" "${CMAKE_INSTALL_PREFIX}" MFEM_LIB_DIR)
mfem_path_to_fullpath(
"${INSTALL_SHARE_DIR}/mfem/test.mk" "${CMAKE_INSTALL_PREFIX}" MFEM_TEST_MK)
set(MFEM_INCFLAGS "-I\$(MFEM_INC_DIR) \$(MFEM_TPLFLAGS)")
set(MFEM_CONFIG_EXTRA "")
# Create the install-tree version of 'config.mk'
-131
View File
@@ -1,131 +0,0 @@
// Define the cube sizes
L_outer = 1.0;
L_inner = 0.5;
// Set mesh size and algorithm
mesh_size = 0.4;
Mesh.Algorithm3D = 1; // Delaunay algorithm for 3D mesh
Mesh.CharacteristicLengthFactor = 1.0;
Mesh.MshFileVersion = 2.2;
// Define center point for concentric cubes
cx = 0.5;
cy = 0.5;
cz = 0.5;
// Define the points (vertices of the outer cube)
Point(1) = {cx-L_outer/2, cy-L_outer/2, cz-L_outer/2, mesh_size};
Point(2) = {cx+L_outer/2, cy-L_outer/2, cz-L_outer/2, mesh_size};
Point(3) = {cx+L_outer/2, cy+L_outer/2, cz-L_outer/2, mesh_size};
Point(4) = {cx-L_outer/2, cy+L_outer/2, cz-L_outer/2, mesh_size};
Point(5) = {cx-L_outer/2, cy-L_outer/2, cz+L_outer/2, mesh_size};
Point(6) = {cx+L_outer/2, cy-L_outer/2, cz+L_outer/2, mesh_size};
Point(7) = {cx+L_outer/2, cy+L_outer/2, cz+L_outer/2, mesh_size};
Point(8) = {cx-L_outer/2, cy+L_outer/2, cz+L_outer/2, mesh_size};
// Define the points (vertices of the inner cube)
Point(9) = {cx-L_inner/2, cy-L_inner/2, cz-L_inner/2, mesh_size};
Point(10) = {cx+L_inner/2, cy-L_inner/2, cz-L_inner/2, mesh_size};
Point(11) = {cx+L_inner/2, cy+L_inner/2, cz-L_inner/2, mesh_size};
Point(12) = {cx-L_inner/2, cy+L_inner/2, cz-L_inner/2, mesh_size};
Point(13) = {cx-L_inner/2, cy-L_inner/2, cz+L_inner/2, mesh_size};
Point(14) = {cx+L_inner/2, cy-L_inner/2, cz+L_inner/2, mesh_size};
Point(15) = {cx+L_inner/2, cy+L_inner/2, cz+L_inner/2, mesh_size};
Point(16) = {cx-L_inner/2, cy+L_inner/2, cz+L_inner/2, mesh_size};
// Define the lines (edges of the outer cube)
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};
Line(9) = {1, 5};
Line(10) = {2, 6};
Line(11) = {3, 7};
Line(12) = {4, 8};
// Define the lines (edges of the inner cube)
Line(13) = {9, 10};
Line(14) = {10, 11};
Line(15) = {11, 12};
Line(16) = {12, 9};
Line(17) = {13, 14};
Line(18) = {14, 15};
Line(19) = {15, 16};
Line(20) = {16, 13};
Line(21) = {9, 13};
Line(22) = {10, 14};
Line(23) = {11, 15};
Line(24) = {12, 16};
// Define the surfaces (faces of the outer cube)
Line Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Line Loop(2) = {5, 6, 7, 8};
Plane Surface(2) = {2};
Line Loop(3) = {9, 5, -10, -1};
Plane Surface(3) = {3};
Line Loop(4) = {10, 6, -11, -2};
Plane Surface(4) = {4};
Line Loop(5) = {11, 7, -12, -3};
Plane Surface(5) = {5};
Line Loop(6) = {12, 8, -9, -4};
Plane Surface(6) = {6};
// Define the surfaces (faces of the inner cube)
Line Loop(7) = {13, 14, 15, 16};
Plane Surface(7) = {7};
Line Loop(8) = {17, 18, 19, 20};
Plane Surface(8) = {8};
Line Loop(9) = {21, 17, -22, -13};
Plane Surface(9) = {9};
Line Loop(10) = {22, 18, -23, -14};
Plane Surface(10) = {10};
Line Loop(11) = {23, 19, -24, -15};
Plane Surface(11) = {11};
Line Loop(12) = {24, 20, -21, -16};
Plane Surface(12) = {12};
// Define the volumes
Surface Loop(1) = {1, 2, 3, 4, 5, 6};
Surface Loop(2) = {7, 8, 9, 10, 11, 12};
Volume(1) = {1, 2}; // Outer volume with inner hole
Volume(2) = {2}; // Inner volume
// Assign physical groups
Physical Volume(1) = {1}; // Outer volume
Physical Volume(2) = {2}; // Inner volume
// Outer cube surfaces
Physical Surface(1) = {1}; // Outer bottom
Physical Surface(2) = {2}; // Outer top
Physical Surface(3) = {3}; // Outer front
Physical Surface(4) = {4}; // Outer right
Physical Surface(5) = {5}; // Outer back
Physical Surface(6) = {6}; // Outer left
// Inner cube surfaces
Physical Surface(7) = {7}; // Inner bottom (-xy)
Physical Surface(8) = {8}; // Inner top (+xy)
Physical Surface(9) = {9}; // Inner front (-xz)
Physical Surface(10) = {10}; // Inner right (+yz)
Physical Surface(11) = {11}; // Inner back (+xz)
Physical Surface(12) = {12}; // Inner left (-yz)
// Mesh control
Mesh.OptimizeNetgen = 1;
Mesh.Optimize = 1;
Mesh.ElementOrder = 1;
-907
View File
@@ -1,907 +0,0 @@
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
138
1 0 0 0
2 1 0 0
3 1 1 0
4 0 1 0
5 0 0 1
6 1 0 1
7 1 1 1
8 0 1 1
9 0.25 0.25 0.25
10 0.75 0.25 0.25
11 0.75 0.75 0.25
12 0.25 0.75 0.25
13 0.25 0.25 0.75
14 0.75 0.25 0.75
15 0.75 0.75 0.75
16 0.25 0.75 0.75
17 0.3333333333325025 0 0
18 0.6666666666657889 0 0
19 1 0.3333333333325025 0
20 1 0.6666666666657889 0
21 0.6666666666675911 1 0
22 0.3333333333347203 1 0
23 0 0.6666666666675911 0
24 0 0.3333333333347203 0
25 0.3333333333325025 0 1
26 0.6666666666657889 0 1
27 1 0.3333333333325025 1
28 1 0.6666666666657889 1
29 0.6666666666675911 1 1
30 0.3333333333347203 1 1
31 0 0.6666666666675911 1
32 0 0.3333333333347203 1
33 0 0 0.3333333333325025
34 0 0 0.6666666666657889
35 1 0 0.3333333333325025
36 1 0 0.6666666666657889
37 1 1 0.3333333333325025
38 1 1 0.6666666666657889
39 0 1 0.3333333333325025
40 0 1 0.6666666666657889
41 0.5000000000003468 0.25 0.25
42 0.75 0.5000000000003468 0.25
43 0.5000000000013763 0.75 0.25
44 0.25 0.5000000000013763 0.25
45 0.5000000000003468 0.25 0.75
46 0.75 0.5000000000003468 0.75
47 0.5000000000013763 0.75 0.75
48 0.25 0.5000000000013763 0.75
49 0.25 0.25 0.5000000000003468
50 0.75 0.25 0.5000000000003468
51 0.75 0.75 0.5000000000003468
52 0.25 0.75 0.5000000000003468
53 0.7113248654055673 0.4999999999991457 0
54 0.2886751345942123 0.5000000000011557 0
55 0.5000000000006117 0.7525600817161773 0
56 0.4999999999993867 0.2474399182839603 0
57 0.2423197548524782 0.7576802451481532 0
58 0.757680245147464 0.2423197548520695 0
59 0.2423197548507857 0.2423197548513912 0
60 0.7576802451491019 0.7576802451486099 0
61 0.7113248654055673 0.4999999999991457 1
62 0.2886751345942123 0.5000000000011557 1
63 0.5000000000006117 0.7525600817161773 1
64 0.4999999999993867 0.2474399182839603 1
65 0.2423197548524782 0.7576802451481532 1
66 0.757680245147464 0.2423197548520695 1
67 0.2423197548507857 0.2423197548513912 1
68 0.7576802451491019 0.7576802451486099 1
69 0.4999999999993203 0 0.301447615129799
70 0.4999999999992795 0 0.7028666213189801
71 0.7525600817158393 0 0.5007190394076877
72 0.2474399182836191 0 0.5007190394076877
73 0.7576802451479793 0 0.7576802451479793
74 0.2423197548517962 0 0.7576802451477375
75 0.7576802451484569 0 0.2423197548510767
76 0.2423197548513188 0 0.2423197548513187
77 1 0.4999999999993203 0.301447615129799
78 1 0.4999999999992795 0.7028666213189801
79 1 0.7525600817158394 0.5007190394076877
80 1 0.2474399182836191 0.5007190394076877
81 1 0.7576802451479794 0.7576802451479794
82 1 0.2423197548517962 0.7576802451477376
83 1 0.7576802451484569 0.2423197548510768
84 1 0.2423197548513188 0.2423197548513188
85 0.5000000000008327 1 0.3014476151298047
86 0.500000000000961 1 0.7028666213191928
87 0.2474399182842484 1 0.5007190394077241
88 0.7525600817164384 1 0.5007190394078933
89 0.2423197548520873 1 0.7576802451480517
90 0.7576802451481496 1 0.2423197548518761
91 0.2423197548516099 1 0.2423197548510044
92 0.7576802451486874 1 0.7576802451481952
93 0 0.5000000000008327 0.3014476151298047
94 0 0.500000000000961 0.7028666213191928
95 0 0.2474399182842484 0.5007190394077241
96 0 0.7525600817164384 0.5007190394078933
97 0 0.2423197548520873 0.7576802451480517
98 0 0.7576802451481496 0.2423197548518761
99 0 0.2423197548516099 0.2423197548510044
100 0 0.7576802451486874 0.7576802451481952
101 0.3968750000003409 0.603125000000244 0.25
102 0.4374999999998713 0.4375000000001287 0.25
103 0.5739583333335919 0.5718750000001767 0.25
104 0.6093749999999631 0.3906250000003402 0.25
105 0.3968750000003409 0.603125000000244 0.75
106 0.4374999999998713 0.4375000000001287 0.75
107 0.5739583333335919 0.5718750000001767 0.75
108 0.6093749999999631 0.3906250000003402 0.75
109 0.3806942419826734 0.25 0.3806942419826734
110 0.5625000000001735 0.25 0.4375000000000001
111 0.4254282069971791 0.25 0.5712615403304835
112 0.6093749999998808 0.25 0.6093749999998808
113 0.75 0.3806942419826734 0.3806942419826734
114 0.75 0.5625000000001735 0.4375000000000001
115 0.75 0.4254282069971791 0.5712615403304835
116 0.75 0.6093749999998808 0.6093749999998808
117 0.3968750000004991 0.75 0.3968750000002804
118 0.4375000000000001 0.75 0.5625000000004308
119 0.5739583333336153 0.75 0.4281250000000707
120 0.6093749999998166 0.75 0.6093749999993661
121 0.25 0.3968750000004991 0.3968750000002804
122 0.25 0.4375000000000001 0.5625000000004308
123 0.25 0.5739583333336153 0.4281250000000707
124 0.25 0.6093749999998166 0.6093749999993661
125 0.4962939304035875 0.5214350017087855 0.4925553109323813
126 0.3432581985549767 0.6471275530923523 0.3554206606168006
127 0.6442168181713744 0.5929232373214773 0.3593785632300704
128 0.625174517421737 0.3491579444372839 0.3604636129462503
129 0.6130544111091688 0.6576364353993367 0.5046720648819553
130 0.4281518698369243 0.3632662294430484 0.3548726263205102
131 0.3639383531355198 0.3520949250221278 0.4978389662613994
132 0.629585530087249 0.3489878162230438 0.5124654846339122
133 0.3710853378652663 0.6517586292698121 0.6382643241302075
134 0.5917018263727056 0.6522525456211955 0.6390764961119443
135 0.3530810314228338 0.4582477062424107 0.6430239639637545
136 0.6571010904289998 0.5300774811423468 0.6603970500567977
137 0.6484596018596915 0.361399676127967 0.6360267588157177
138 0.4782020887035478 0.3534611476388013 0.6141275027013793
$EndNodes
$Elements
760
1 2 2 1 1 1 17 59
2 2 2 1 1 24 1 59
3 2 2 1 1 18 2 58
4 2 2 1 1 2 19 58
5 2 2 1 1 20 3 60
6 2 2 1 1 3 21 60
7 2 2 1 1 22 4 57
8 2 2 1 1 4 23 57
9 2 2 1 1 17 18 56
10 2 2 1 1 17 56 59
11 2 2 1 1 56 18 58
12 2 2 1 1 19 20 53
13 2 2 1 1 19 53 58
14 2 2 1 1 53 20 60
15 2 2 1 1 21 22 55
16 2 2 1 1 21 55 60
17 2 2 1 1 55 22 57
18 2 2 1 1 23 24 54
19 2 2 1 1 23 54 57
20 2 2 1 1 54 24 59
21 2 2 1 1 54 53 55
22 2 2 1 1 53 54 56
23 2 2 1 1 55 53 60
24 2 2 1 1 53 56 58
25 2 2 1 1 54 55 57
26 2 2 1 1 56 54 59
27 2 2 2 2 5 25 67
28 2 2 2 2 32 5 67
29 2 2 2 2 26 6 66
30 2 2 2 2 6 27 66
31 2 2 2 2 28 7 68
32 2 2 2 2 7 29 68
33 2 2 2 2 30 8 65
34 2 2 2 2 8 31 65
35 2 2 2 2 25 26 64
36 2 2 2 2 25 64 67
37 2 2 2 2 64 26 66
38 2 2 2 2 27 28 61
39 2 2 2 2 27 61 66
40 2 2 2 2 61 28 68
41 2 2 2 2 29 30 63
42 2 2 2 2 29 63 68
43 2 2 2 2 63 30 65
44 2 2 2 2 31 32 62
45 2 2 2 2 31 62 65
46 2 2 2 2 62 32 67
47 2 2 2 2 62 61 63
48 2 2 2 2 61 62 64
49 2 2 2 2 63 61 68
50 2 2 2 2 61 64 66
51 2 2 2 2 62 63 65
52 2 2 2 2 64 62 67
53 2 2 3 3 17 1 76
54 2 2 3 3 1 33 76
55 2 2 3 3 2 18 75
56 2 2 3 3 35 2 75
57 2 2 3 3 5 25 74
58 2 2 3 3 34 5 74
59 2 2 3 3 26 6 73
60 2 2 3 3 6 36 73
61 2 2 3 3 18 17 69
62 2 2 3 3 69 17 76
63 2 2 3 3 18 69 75
64 2 2 3 3 25 26 70
65 2 2 3 3 25 70 74
66 2 2 3 3 70 26 73
67 2 2 3 3 33 34 72
68 2 2 3 3 33 72 76
69 2 2 3 3 72 34 74
70 2 2 3 3 36 35 71
71 2 2 3 3 71 35 75
72 2 2 3 3 36 71 73
73 2 2 3 3 69 70 71
74 2 2 3 3 70 69 72
75 2 2 3 3 69 71 75
76 2 2 3 3 72 69 76
77 2 2 3 3 71 70 73
78 2 2 3 3 70 72 74
79 2 2 4 4 19 2 84
80 2 2 4 4 2 35 84
81 2 2 4 4 3 20 83
82 2 2 4 4 37 3 83
83 2 2 4 4 6 27 82
84 2 2 4 4 36 6 82
85 2 2 4 4 28 7 81
86 2 2 4 4 7 38 81
87 2 2 4 4 20 19 77
88 2 2 4 4 77 19 84
89 2 2 4 4 20 77 83
90 2 2 4 4 27 28 78
91 2 2 4 4 27 78 82
92 2 2 4 4 78 28 81
93 2 2 4 4 35 36 80
94 2 2 4 4 35 80 84
95 2 2 4 4 80 36 82
96 2 2 4 4 38 37 79
97 2 2 4 4 79 37 83
98 2 2 4 4 38 79 81
99 2 2 4 4 77 78 79
100 2 2 4 4 78 77 80
101 2 2 4 4 77 79 83
102 2 2 4 4 80 77 84
103 2 2 4 4 79 78 81
104 2 2 4 4 78 80 82
105 2 2 5 5 21 3 90
106 2 2 5 5 3 37 90
107 2 2 5 5 4 22 91
108 2 2 5 5 39 4 91
109 2 2 5 5 7 29 92
110 2 2 5 5 38 7 92
111 2 2 5 5 30 8 89
112 2 2 5 5 8 40 89
113 2 2 5 5 22 21 85
114 2 2 5 5 85 21 90
115 2 2 5 5 22 85 91
116 2 2 5 5 29 30 86
117 2 2 5 5 29 86 92
118 2 2 5 5 86 30 89
119 2 2 5 5 37 38 88
120 2 2 5 5 37 88 90
121 2 2 5 5 88 38 92
122 2 2 5 5 40 39 87
123 2 2 5 5 87 39 91
124 2 2 5 5 40 87 89
125 2 2 5 5 85 86 87
126 2 2 5 5 86 85 88
127 2 2 5 5 85 87 91
128 2 2 5 5 88 85 90
129 2 2 5 5 87 86 89
130 2 2 5 5 86 88 92
131 2 2 6 6 1 24 99
132 2 2 6 6 33 1 99
133 2 2 6 6 23 4 98
134 2 2 6 6 4 39 98
135 2 2 6 6 32 5 97
136 2 2 6 6 5 34 97
137 2 2 6 6 8 31 100
138 2 2 6 6 40 8 100
139 2 2 6 6 24 23 93
140 2 2 6 6 93 23 98
141 2 2 6 6 24 93 99
142 2 2 6 6 31 32 94
143 2 2 6 6 31 94 100
144 2 2 6 6 94 32 97
145 2 2 6 6 34 33 95
146 2 2 6 6 95 33 99
147 2 2 6 6 34 95 97
148 2 2 6 6 39 40 96
149 2 2 6 6 39 96 98
150 2 2 6 6 96 40 100
151 2 2 6 6 93 94 95
152 2 2 6 6 94 93 96
153 2 2 6 6 93 95 99
154 2 2 6 6 96 93 98
155 2 2 6 6 95 94 97
156 2 2 6 6 94 96 100
157 2 2 7 7 9 41 102
158 2 2 7 7 44 9 102
159 2 2 7 7 41 10 104
160 2 2 7 7 10 42 104
161 2 2 7 7 42 11 103
162 2 2 7 7 11 43 103
163 2 2 7 7 43 12 101
164 2 2 7 7 12 44 101
165 2 2 7 7 102 41 104
166 2 2 7 7 42 103 104
167 2 2 7 7 43 101 103
168 2 2 7 7 101 44 102
169 2 2 7 7 101 102 103
170 2 2 7 7 103 102 104
171 2 2 8 8 13 45 106
172 2 2 8 8 48 13 106
173 2 2 8 8 45 14 108
174 2 2 8 8 14 46 108
175 2 2 8 8 46 15 107
176 2 2 8 8 15 47 107
177 2 2 8 8 47 16 105
178 2 2 8 8 16 48 105
179 2 2 8 8 106 45 108
180 2 2 8 8 46 107 108
181 2 2 8 8 47 105 107
182 2 2 8 8 105 48 106
183 2 2 8 8 105 106 107
184 2 2 8 8 107 106 108
185 2 2 9 9 41 9 109
186 2 2 9 9 9 49 109
187 2 2 9 9 10 41 110
188 2 2 9 9 50 10 110
189 2 2 9 9 13 45 111
190 2 2 9 9 49 13 111
191 2 2 9 9 45 14 112
192 2 2 9 9 14 50 112
193 2 2 9 9 41 109 110
194 2 2 9 9 111 45 112
195 2 2 9 9 109 49 111
196 2 2 9 9 50 110 112
197 2 2 9 9 110 109 111
198 2 2 9 9 110 111 112
199 2 2 10 10 42 10 113
200 2 2 10 10 10 50 113
201 2 2 10 10 11 42 114
202 2 2 10 10 51 11 114
203 2 2 10 10 14 46 115
204 2 2 10 10 50 14 115
205 2 2 10 10 46 15 116
206 2 2 10 10 15 51 116
207 2 2 10 10 42 113 114
208 2 2 10 10 115 46 116
209 2 2 10 10 113 50 115
210 2 2 10 10 51 114 116
211 2 2 10 10 114 113 115
212 2 2 10 10 114 115 116
213 2 2 11 11 43 11 119
214 2 2 11 11 11 51 119
215 2 2 11 11 12 43 117
216 2 2 11 11 52 12 117
217 2 2 11 11 15 47 120
218 2 2 11 11 51 15 120
219 2 2 11 11 47 16 118
220 2 2 11 11 16 52 118
221 2 2 11 11 117 43 119
222 2 2 11 11 47 118 120
223 2 2 11 11 119 51 120
224 2 2 11 11 52 117 118
225 2 2 11 11 118 117 119
226 2 2 11 11 118 119 120
227 2 2 12 12 9 44 121
228 2 2 12 12 49 9 121
229 2 2 12 12 44 12 123
230 2 2 12 12 12 52 123
231 2 2 12 12 48 13 122
232 2 2 12 12 13 49 122
233 2 2 12 12 16 48 124
234 2 2 12 12 52 16 124
235 2 2 12 12 121 44 123
236 2 2 12 12 48 122 124
237 2 2 12 12 49 121 122
238 2 2 12 12 123 52 124
239 2 2 12 12 122 121 123
240 2 2 12 12 122 123 124
241 4 2 1 1 105 62 106 107
242 4 2 1 1 102 54 101 103
243 4 2 1 1 118 86 120 119
244 4 2 1 1 124 94 123 122
245 4 2 1 1 52 39 12 96
246 4 2 1 1 52 12 39 87
247 4 2 1 1 88 38 15 51
248 4 2 1 1 79 15 38 51
249 4 2 1 1 80 14 50 36
250 4 2 1 1 71 50 14 36
251 4 2 1 1 120 88 15 51
252 4 2 1 1 116 15 79 51
253 4 2 1 1 50 14 112 71
254 4 2 1 1 52 117 12 87
255 4 2 1 1 111 69 109 110
256 4 2 1 1 114 77 115 113
257 4 2 1 1 124 123 94 96
258 4 2 1 1 120 86 88 119
259 4 2 1 1 14 45 64 108
260 4 2 1 1 103 54 101 55
261 4 2 1 1 62 105 63 107
262 4 2 1 1 63 29 15 47
263 4 2 1 1 14 45 26 64
264 4 2 1 1 51 11 88 37
265 4 2 1 1 80 50 10 35
266 4 2 1 1 49 33 9 72
267 4 2 1 1 71 10 50 35
268 4 2 1 1 79 11 51 37
269 4 2 1 1 49 9 33 95
270 4 2 1 1 12 123 52 96
271 4 2 1 1 96 16 52 40
272 4 2 1 1 49 13 34 72
273 4 2 1 1 87 52 16 40
274 4 2 1 1 13 49 34 95
275 4 2 1 1 43 101 12 55
276 4 2 1 1 43 12 22 55
277 4 2 1 1 106 61 108 107
278 4 2 1 1 102 103 104 53
279 4 2 1 1 50 115 14 80
280 4 2 1 1 80 10 50 113
281 4 2 1 1 9 109 49 72
282 4 2 1 1 47 15 63 107
283 4 2 1 1 124 16 52 96
284 4 2 1 1 49 121 9 95
285 4 2 1 1 15 28 68 46
286 4 2 1 1 81 28 15 46
287 4 2 1 1 92 15 29 47
288 4 2 1 1 27 82 14 46
289 4 2 1 1 14 66 27 46
290 4 2 1 1 26 45 14 73
291 4 2 1 1 78 115 116 114
292 4 2 1 1 112 70 111 110
293 4 2 1 1 79 11 114 51
294 4 2 1 1 71 50 10 110
295 4 2 1 1 10 41 104 56
296 4 2 1 1 77 80 115 113
297 4 2 1 1 69 109 72 111
298 4 2 1 1 63 16 30 47
299 4 2 1 1 25 45 13 64
300 4 2 1 1 52 16 118 87
301 4 2 1 1 95 13 49 122
302 4 2 1 1 57 12 23 44
303 4 2 1 1 119 11 88 51
304 4 2 1 1 41 17 9 56
305 4 2 1 1 18 41 10 56
306 4 2 1 1 21 11 43 55
307 4 2 1 1 63 105 16 47
308 4 2 1 1 13 45 106 64
309 4 2 1 1 9 102 41 56
310 4 2 1 1 72 49 13 111
311 4 2 1 1 11 103 43 55
312 4 2 1 1 70 112 71 110
313 4 2 1 1 116 79 78 114
314 4 2 1 1 32 13 48 67
315 4 2 1 1 16 89 30 47
316 4 2 1 1 32 48 13 97
317 4 2 1 1 48 16 100 31
318 4 2 1 1 121 123 93 122
319 4 2 1 1 119 117 118 85
320 4 2 1 1 10 19 42 84
321 4 2 1 1 17 9 76 41
322 4 2 1 1 24 9 59 44
323 4 2 1 1 42 19 10 58
324 4 2 1 1 10 41 18 75
325 4 2 1 1 20 83 11 42
326 4 2 1 1 90 11 43 21
327 4 2 1 1 53 104 102 56
328 4 2 1 1 106 61 64 108
329 4 2 1 1 87 118 117 85
330 4 2 1 1 95 121 93 122
331 4 2 1 1 80 50 115 113
332 4 2 1 1 109 49 72 111
333 4 2 1 1 78 27 28 46
334 4 2 1 1 31 94 48 32
335 4 2 1 1 19 77 20 42
336 4 2 1 1 23 24 44 93
337 4 2 1 1 123 124 52 96
338 4 2 1 1 88 120 119 51
339 4 2 1 1 91 43 85 117
340 4 2 1 1 39 98 12 96
341 4 2 1 1 12 91 39 87
342 4 2 1 1 15 92 38 88
343 4 2 1 1 38 81 15 79
344 4 2 1 1 36 80 14 82
345 4 2 1 1 14 71 36 73
346 4 2 1 1 103 102 54 53
347 4 2 1 1 61 106 62 107
348 4 2 1 1 86 29 30 47
349 4 2 1 1 26 45 70 25
350 4 2 1 1 21 85 22 43
351 4 2 1 1 94 93 123 122
352 4 2 1 1 119 118 86 85
353 4 2 1 1 12 96 93 123
354 4 2 1 1 115 78 77 114
355 4 2 1 1 69 111 70 110
356 4 2 1 1 28 27 61 46
357 4 2 1 1 48 62 31 32
358 4 2 1 1 88 15 92 120
359 4 2 1 1 79 81 15 116
360 4 2 1 1 71 14 112 73
361 4 2 1 1 114 116 79 51
362 4 2 1 1 112 50 71 110
363 4 2 1 1 91 12 117 87
364 4 2 1 1 80 115 14 82
365 4 2 1 1 63 15 29 68
366 4 2 1 1 64 26 14 66
367 4 2 1 1 47 63 105 107
368 4 2 1 1 43 103 101 55
369 4 2 1 1 35 10 80 84
370 4 2 1 1 33 76 9 72
371 4 2 1 1 37 11 88 90
372 4 2 1 1 35 71 10 75
373 4 2 1 1 37 79 11 83
374 4 2 1 1 9 99 33 95
375 4 2 1 1 68 63 15 107
376 4 2 1 1 99 44 93 121
377 4 2 1 1 18 69 41 17
378 4 2 1 1 10 35 2 84
379 4 2 1 1 15 7 28 81
380 4 2 1 1 3 37 90 11
381 4 2 1 1 108 64 14 66
382 4 2 1 1 51 38 79 37
383 4 2 1 1 12 4 23 98
384 4 2 1 1 12 39 91 4
385 4 2 1 1 49 33 72 34
386 4 2 1 1 10 71 110 75
387 4 2 1 1 114 11 79 83
388 4 2 1 1 40 96 16 100
389 4 2 1 1 34 13 74 72
390 4 2 1 1 16 87 40 89
391 4 2 1 1 34 97 13 95
392 4 2 1 1 20 19 42 53
393 4 2 1 1 54 23 24 44
394 4 2 1 1 87 16 118 89
395 4 2 1 1 95 97 13 122
396 4 2 1 1 57 22 12 55
397 4 2 1 1 41 102 104 56
398 4 2 1 1 64 45 106 108
399 4 2 1 1 101 57 12 55
400 4 2 1 1 74 13 5 25
401 4 2 1 1 13 34 74 5
402 4 2 1 1 8 40 16 100
403 4 2 1 1 16 31 65 8
404 4 2 1 1 82 36 6 14
405 4 2 1 1 102 9 59 56
406 4 2 1 1 106 67 13 64
407 4 2 1 1 119 88 11 90
408 4 2 1 1 95 49 121 122
409 4 2 1 1 52 118 117 87
410 4 2 1 1 124 16 96 100
411 4 2 1 1 80 10 113 84
412 4 2 1 1 109 9 76 72
413 4 2 1 1 72 13 70 111
414 4 2 1 1 99 9 121 95
415 4 2 1 1 30 65 16 63
416 4 2 1 1 13 67 25 64
417 4 2 1 1 11 55 53 103
418 4 2 1 1 123 12 44 93
419 4 2 1 1 56 104 10 58
420 4 2 1 1 10 18 56 58
421 4 2 1 1 21 60 11 55
422 4 2 1 1 17 59 9 56
423 4 2 1 1 16 63 62 105
424 4 2 1 1 48 105 16 62
425 4 2 1 1 111 45 13 70
426 4 2 1 1 119 86 88 85
427 4 2 1 1 123 93 94 96
428 4 2 1 1 115 78 80 77
429 4 2 1 1 111 72 69 70
430 4 2 1 1 55 103 54 53
431 4 2 1 1 63 61 62 107
432 4 2 1 1 14 27 66 6
433 4 2 1 1 99 9 1 24
434 4 2 1 1 49 33 34 95
435 4 2 1 1 96 52 39 40
436 4 2 1 1 22 12 91 4
437 4 2 1 1 20 11 83 3
438 4 2 1 1 29 15 92 7
439 4 2 1 1 18 10 75 2
440 4 2 1 1 36 50 80 35
441 4 2 1 1 18 17 41 56
442 4 2 1 1 51 88 38 37
443 4 2 1 1 62 31 16 48
444 4 2 1 1 62 16 31 65
445 4 2 1 1 26 45 25 64
446 4 2 1 1 104 58 42 10
447 4 2 1 1 29 63 30 47
448 4 2 1 1 77 78 79 114
449 4 2 1 1 71 69 70 110
450 4 2 1 1 42 103 11 53
451 4 2 1 1 93 94 95 122
452 4 2 1 1 118 87 86 85
453 4 2 1 1 61 62 106 64
454 4 2 1 1 53 102 54 56
455 4 2 1 1 78 82 27 46
456 4 2 1 1 94 48 32 97
457 4 2 1 1 86 30 89 47
458 4 2 1 1 78 28 81 46
459 4 2 1 1 26 70 45 73
460 4 2 1 1 100 94 48 31
461 4 2 1 1 86 92 29 47
462 4 2 1 1 76 69 17 41
463 4 2 1 1 42 19 77 84
464 4 2 1 1 77 83 20 42
465 4 2 1 1 18 41 69 75
466 4 2 1 1 43 85 90 21
467 4 2 1 1 91 43 117 12
468 4 2 1 1 48 62 32 67
469 4 2 1 1 27 66 61 46
470 4 2 1 1 68 28 61 46
471 4 2 1 1 54 24 59 44
472 4 2 1 1 19 42 53 58
473 4 2 1 1 23 54 57 44
474 4 2 1 1 12 93 96 98
475 4 2 1 1 99 44 121 9
476 4 2 1 1 42 58 104 53
477 4 2 1 1 21 43 22 55
478 4 2 1 1 53 20 11 42
479 4 2 1 1 53 11 20 60
480 4 2 1 1 106 67 48 13
481 4 2 1 1 36 71 50 35
482 4 2 1 1 52 87 39 40
483 4 2 1 1 48 67 106 62
484 4 2 1 1 72 70 13 74
485 4 2 1 1 9 33 99 1
486 4 2 1 1 3 37 11 83
487 4 2 1 1 13 34 5 97
488 4 2 1 1 12 39 4 98
489 4 2 1 1 14 36 6 73
490 4 2 1 1 15 38 92 7
491 4 2 1 1 8 40 89 16
492 4 2 1 1 75 35 2 10
493 4 2 1 1 30 16 65 8
494 4 2 1 1 14 27 6 82
495 4 2 1 1 26 14 66 6
496 4 2 1 1 16 31 8 100
497 4 2 1 1 32 13 67 5
498 4 2 1 1 13 67 5 25
499 4 2 1 1 68 7 28 15
500 4 2 1 1 29 15 7 68
501 4 2 1 1 57 4 23 12
502 4 2 1 1 22 12 4 57
503 4 2 1 1 10 19 84 2
504 4 2 1 1 9 59 1 24
505 4 2 1 1 17 9 59 1
506 4 2 1 1 18 10 2 58
507 4 2 1 1 11 21 90 3
508 4 2 1 1 20 11 3 60
509 4 2 1 1 16 62 63 65
510 4 2 1 1 11 53 55 60
511 4 2 1 1 78 81 79 116
512 4 2 1 1 71 112 70 73
513 4 2 1 1 86 88 92 120
514 4 2 1 1 96 94 124 100
515 4 2 1 1 94 97 95 122
516 4 2 1 1 89 118 87 86
517 4 2 1 1 15 38 7 81
518 4 2 1 1 102 44 59 9
519 4 2 1 1 59 44 102 54
520 4 2 1 1 110 71 69 75
521 4 2 1 1 79 77 114 83
522 4 2 1 1 76 69 109 72
523 4 2 1 1 113 77 80 84
524 4 2 1 1 93 99 121 95
525 4 2 1 1 87 117 91 85
526 4 2 1 1 90 119 43 11
527 4 2 1 1 80 78 115 82
528 4 2 1 1 88 119 85 90
529 4 2 1 1 90 43 119 85
530 4 2 1 1 46 108 66 61
531 4 2 1 1 46 66 108 14
532 4 2 1 1 26 14 6 73
533 4 2 1 1 44 101 57 12
534 4 2 1 1 44 57 101 54
535 4 2 1 1 46 107 68 15
536 4 2 1 1 46 68 107 61
537 4 2 1 1 9 33 1 76
538 4 2 1 1 100 124 48 94
539 4 2 1 1 100 48 124 16
540 4 2 1 1 25 70 13 45
541 4 2 1 1 56 53 104 58
542 4 2 1 1 61 64 108 66
543 4 2 1 1 106 62 67 64
544 4 2 1 1 102 59 54 56
545 4 2 1 1 32 13 5 97
546 4 2 1 1 30 16 8 89
547 4 2 1 1 13 70 25 74
548 4 2 1 1 97 122 48 13
549 4 2 1 1 75 110 41 69
550 4 2 1 1 54 57 101 55
551 4 2 1 1 97 48 122 94
552 4 2 1 1 47 118 89 86
553 4 2 1 1 47 89 118 16
554 4 2 1 1 17 9 1 76
555 4 2 1 1 11 21 3 60
556 4 2 1 1 10 19 2 58
557 4 2 1 1 75 41 110 10
558 4 2 1 1 61 63 68 107
559 4 2 1 1 70 111 45 112
560 4 2 1 1 78 115 46 116
561 4 2 1 1 77 114 42 113
562 4 2 1 1 69 41 109 110
563 4 2 1 1 107 61 108 46
564 4 2 1 1 103 42 104 53
565 4 2 1 1 118 86 47 120
566 4 2 1 1 94 124 48 122
567 4 2 1 1 117 119 43 85
568 4 2 1 1 123 44 121 93
569 4 2 1 1 101 54 102 44
570 4 2 1 1 48 105 62 106
571 4 2 1 1 91 43 12 22
572 4 2 1 1 91 43 22 85
573 4 2 1 1 93 12 23 98
574 4 2 1 1 93 23 12 44
575 4 2 1 1 46 81 116 15
576 4 2 1 1 46 116 81 78
577 4 2 1 1 92 47 120 15
578 4 2 1 1 120 47 92 86
579 4 2 1 1 46 82 115 78
580 4 2 1 1 46 115 82 14
581 4 2 1 1 73 45 112 70
582 4 2 1 1 73 112 45 14
583 4 2 1 1 99 44 9 24
584 4 2 1 1 99 44 24 93
585 4 2 1 1 84 42 113 77
586 4 2 1 1 84 113 42 10
587 4 2 1 1 76 41 109 69
588 4 2 1 1 109 41 76 9
589 4 2 1 1 42 83 114 77
590 4 2 1 1 42 114 83 11
591 4 2 2 2 135 13 122 131
592 4 2 2 2 132 138 125 110
593 4 2 2 2 138 108 125 106
594 4 2 2 2 136 108 125 137
595 4 2 2 2 110 138 125 131
596 4 2 2 2 136 107 134 125
597 4 2 2 2 122 135 131 125
598 4 2 2 2 13 49 122 131
599 4 2 2 2 136 108 137 46
600 4 2 2 2 110 138 131 111
601 4 2 2 2 132 138 110 112
602 4 2 2 2 106 107 108 125
603 4 2 2 2 138 112 137 45
604 4 2 2 2 104 127 125 128
605 4 2 2 2 113 115 114 125
606 4 2 2 2 131 122 125 121
607 4 2 2 2 112 45 14 137
608 4 2 2 2 110 131 125 130
609 4 2 2 2 115 116 114 125
610 4 2 2 2 132 138 112 137
611 4 2 2 2 50 137 115 132
612 4 2 2 2 104 127 128 42
613 4 2 2 2 105 134 125 133
614 4 2 2 2 129 120 51 119
615 4 2 2 2 104 130 128 125
616 4 2 2 2 105 134 133 47
617 4 2 2 2 135 106 138 125
618 4 2 2 2 106 105 107 125
619 4 2 2 2 138 108 106 45
620 4 2 2 2 104 130 125 102
621 4 2 2 2 104 127 42 103
622 4 2 2 2 109 131 49 111
623 4 2 2 2 106 135 48 105
624 4 2 2 2 126 101 43 125
625 4 2 2 2 138 137 125 108
626 4 2 2 2 104 127 103 125
627 4 2 2 2 116 114 129 51
628 4 2 2 2 136 108 46 107
629 4 2 2 2 136 108 107 125
630 4 2 2 2 52 124 123 125
631 4 2 2 2 132 112 110 50
632 4 2 2 2 44 121 125 123
633 4 2 2 2 131 122 121 49
634 4 2 2 2 105 134 47 107
635 4 2 2 2 106 13 135 138
636 4 2 2 2 52 133 125 118
637 4 2 2 2 103 43 101 125
638 4 2 2 2 117 43 119 125
639 4 2 2 2 105 134 107 125
640 4 2 2 2 137 115 132 125
641 4 2 2 2 50 137 132 112
642 4 2 2 2 106 135 105 125
643 4 2 2 2 102 101 44 125
644 4 2 2 2 104 130 102 41
645 4 2 2 2 138 45 13 111
646 4 2 2 2 131 138 135 13
647 4 2 2 2 138 112 45 111
648 4 2 2 2 131 138 13 111
649 4 2 2 2 107 134 15 136
650 4 2 2 2 118 52 117 125
651 4 2 2 2 101 126 12 44
652 4 2 2 2 104 130 41 128
653 4 2 2 2 138 137 108 45
654 4 2 2 2 133 118 16 47
655 4 2 2 2 44 101 126 125
656 4 2 2 2 133 134 118 47
657 4 2 2 2 133 134 125 118
658 4 2 2 2 124 122 123 125
659 4 2 2 2 129 118 119 125
660 4 2 2 2 102 104 103 125
661 4 2 2 2 135 48 13 106
662 4 2 2 2 123 122 121 125
663 4 2 2 2 103 127 43 125
664 4 2 2 2 106 45 13 138
665 4 2 2 2 103 11 127 42
666 4 2 2 2 127 51 129 114
667 4 2 2 2 13 49 131 111
668 4 2 2 2 137 108 14 46
669 4 2 2 2 15 47 134 107
670 4 2 2 2 10 41 128 104
671 4 2 2 2 16 118 133 52
672 4 2 2 2 46 107 15 136
673 4 2 2 2 41 102 9 130
674 4 2 2 2 128 50 132 110
675 4 2 2 2 10 41 110 128
676 4 2 2 2 50 137 112 14
677 4 2 2 2 130 102 9 44
678 4 2 2 2 105 133 16 47
679 4 2 2 2 127 11 103 43
680 4 2 2 2 128 130 41 110
681 4 2 2 2 116 134 15 51
682 4 2 2 2 137 45 14 108
683 4 2 2 2 12 126 101 43
684 4 2 2 2 133 48 135 105
685 4 2 2 2 128 42 10 104
686 4 2 2 2 131 9 109 49
687 4 2 2 2 118 117 119 125
688 4 2 2 2 102 103 101 125
689 4 2 2 2 129 118 125 134
690 4 2 2 2 125 129 116 114
691 4 2 2 2 117 126 43 125
692 4 2 2 2 126 117 12 52
693 4 2 2 2 52 126 117 125
694 4 2 2 2 12 117 126 43
695 4 2 2 2 127 51 114 11
696 4 2 2 2 114 129 127 125
697 4 2 2 2 127 119 43 125
698 4 2 2 2 119 11 127 43
699 4 2 2 2 127 113 114 125
700 4 2 2 2 127 113 42 114
701 4 2 2 2 128 130 110 125
702 4 2 2 2 15 47 120 134
703 4 2 2 2 127 11 114 42
704 4 2 2 2 120 47 118 134
705 4 2 2 2 44 130 102 125
706 4 2 2 2 44 126 12 123
707 4 2 2 2 123 44 126 125
708 4 2 2 2 129 118 134 120
709 4 2 2 2 110 132 128 125
710 4 2 2 2 128 127 125 113
711 4 2 2 2 128 127 113 42
712 4 2 2 2 128 50 110 10
713 4 2 2 2 113 42 10 128
714 4 2 2 2 129 134 125 116
715 4 2 2 2 127 51 11 119
716 4 2 2 2 129 134 116 51
717 4 2 2 2 127 51 119 129
718 4 2 2 2 129 119 127 125
719 4 2 2 2 110 131 130 109
720 4 2 2 2 109 130 9 131
721 4 2 2 2 121 130 9 44
722 4 2 2 2 122 48 13 135
723 4 2 2 2 44 121 130 125
724 4 2 2 2 110 138 111 112
725 4 2 2 2 131 138 125 135
726 4 2 2 2 115 137 14 46
727 4 2 2 2 126 52 12 123
728 4 2 2 2 123 126 52 125
729 4 2 2 2 50 137 14 115
730 4 2 2 2 136 137 115 46
731 4 2 2 2 121 9 131 49
732 4 2 2 2 121 131 130 125
733 4 2 2 2 131 130 9 121
734 4 2 2 2 113 132 115 125
735 4 2 2 2 113 50 115 132
736 4 2 2 2 128 50 10 113
737 4 2 2 2 132 113 128 125
738 4 2 2 2 128 50 113 132
739 4 2 2 2 133 105 135 125
740 4 2 2 2 52 124 125 133
741 4 2 2 2 133 48 105 16
742 4 2 2 2 16 133 124 52
743 4 2 2 2 136 137 125 115
744 4 2 2 2 132 138 137 125
745 4 2 2 2 41 130 9 109
746 4 2 2 2 134 120 15 51
747 4 2 2 2 129 118 120 119
748 4 2 2 2 129 120 134 51
749 4 2 2 2 135 122 124 125
750 4 2 2 2 135 48 124 122
751 4 2 2 2 133 48 16 124
752 4 2 2 2 133 135 124 125
753 4 2 2 2 133 48 124 135
754 4 2 2 2 116 136 134 125
755 4 2 2 2 115 136 116 125
756 4 2 2 2 46 115 136 116
757 4 2 2 2 136 134 15 116
758 4 2 2 2 46 136 15 116
759 4 2 2 2 109 41 130 110
760 4 2 2 2 110 131 109 111
$EndElements
-77
View File
@@ -1,77 +0,0 @@
// Square-in-square 2D geometry for MFEM
// Creates concentric squares with different material attributes
// Define the square sizes
L_outer = 2.0;
L_inner = 0.5;
// Set mesh size and algorithm
mesh_size = 1.0;
Mesh.Algorithm = 6; // Frontal-Delaunay for 2D triangular mesh
Mesh.CharacteristicLengthFactor = 1.0;
Mesh.MshFileVersion = 2.2;
// Define center point for concentric squares
cx = 0.0;
cy = 0.0;
// Define the points (vertices of the outer square)
Point(1) = {cx-L_outer/2, cy-L_outer/2, 0, mesh_size}; // bottom-left outer
Point(2) = {cx+L_outer/2, cy-L_outer/2, 0, mesh_size}; // bottom-right outer
Point(3) = {cx+L_outer/2, cy+L_outer/2, 0, mesh_size}; // top-right outer
Point(4) = {cx-L_outer/2, cy+L_outer/2, 0, mesh_size}; // top-left outer
// Define the points (vertices of the inner square)
Point(5) = {cx-L_inner/2, cy-L_inner/2, 0, mesh_size}; // bottom-left inner
Point(6) = {cx+L_inner/2, cy-L_inner/2, 0, mesh_size}; // bottom-right inner
Point(7) = {cx+L_inner/2, cy+L_inner/2, 0, mesh_size}; // top-right inner
Point(8) = {cx-L_inner/2, cy+L_inner/2, 0, mesh_size}; // top-left inner
// Define the lines (edges of the outer square)
Line(1) = {1, 2}; // bottom edge
Line(2) = {2, 3}; // right edge
Line(3) = {3, 4}; // top edge
Line(4) = {4, 1}; // left edge
// Define the lines (edges of the inner square)
Line(5) = {5, 6}; // bottom edge
Line(6) = {6, 7}; // right edge
Line(7) = {7, 8}; // top edge
Line(8) = {8, 5}; // left edge
// Define the surfaces
// Outer square boundary
Line Loop(1) = {1, 2, 3, 4};
// Inner square boundary (hole in the outer region)
Line Loop(2) = {5, 6, 7, 8};
// Define the surface areas
// Outer region (annular region between squares)
Plane Surface(1) = {1, 2}; // Outer loop minus inner loop (creates hole)
// Inner region (solid inner square)
Plane Surface(2) = {2}; // Inner loop only
// Assign physical groups for materials
Physical Surface(1) = {1}; // Outer material (annular region)
Physical Surface(2) = {2}; // Inner material (solid square)
// Physical lines for boundary conditions
// Outer square boundary edges
Physical Line(1) = {1}; // outer bottom
Physical Line(2) = {2}; // outer right
Physical Line(3) = {3}; // outer top
Physical Line(4) = {4}; // outer left
// Inner square boundary edges
Physical Line(5) = {5}; // inner bottom
Physical Line(6) = {6}; // inner right
Physical Line(7) = {7}; // inner top
Physical Line(8) = {8}; // inner left
// Mesh control for quality
Mesh.OptimizeNetgen = 1;
Mesh.Optimize = 1;
Mesh.ElementOrder = 1;
Mesh.RecombineAll = 0; // Keep triangular elements (don't recombine to quads)
-50
View File
@@ -1,50 +0,0 @@
$MeshFormat
2.2 0 8
$EndMeshFormat
$Nodes
13
1 -1 -1 0
2 1 -1 0
3 1 1 0
4 -1 1 0
5 -0.25 -0.25 0
6 0.25 -0.25 0
7 0.25 0.25 0
8 -0.25 0.25 0
9 -2.752797989558076e-12 -1 0
10 1 -2.752797989558076e-12 0
11 2.752797989558076e-12 1 0
12 -1 2.752797989558076e-12 0
13 0 0 0
$EndNodes
$Elements
28
1 1 2 1 1 1 9
2 1 2 1 1 9 2
3 1 2 2 2 2 10
4 1 2 2 2 10 3
5 1 2 3 3 3 11
6 1 2 3 3 11 4
7 1 2 4 4 4 12
8 1 2 4 4 12 1
9 1 2 5 5 5 6
10 1 2 6 6 6 7
11 1 2 7 7 7 8
12 1 2 8 8 8 5
13 2 2 1 1 6 5 9
14 2 2 1 1 5 8 12
15 2 2 1 1 7 6 10
16 2 2 1 1 8 7 11
17 2 2 1 1 9 5 1
18 2 2 1 1 5 12 1
19 2 2 1 1 6 9 2
20 2 2 1 1 10 6 2
21 2 2 1 1 7 10 3
22 2 2 1 1 11 7 3
23 2 2 1 1 8 11 4
24 2 2 1 1 8 4 12
25 2 2 2 2 5 6 13
26 2 2 2 2 8 5 13
27 2 2 2 2 6 7 13
28 2 2 2 2 7 8 13
$EndElements
+5
View File
@@ -262,3 +262,8 @@ endif()
if(MFEM_USE_MOONOLITH)
add_subdirectory(moonolith)
endif()
# Include the examples/magma directory if MAGMA is enabled.
if (MFEM_USE_MAGMA)
add_subdirectory(magma)
endif()
+27
View File
@@ -0,0 +1,27 @@
# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
# LICENSE and NOTICE for details. LLNL-CODE-806117.
#
# This file is part of the MFEM library. For more information and source code
# availability visit https://mfem.org.
#
# MFEM is free software; you can redistribute it and/or modify it under the
# terms of the BSD-3 license. We welcome feedback and contributions, see file
# CONTRIBUTING.md for details.
set(MAGMA_EXAMPLES_SRCS)
list(APPEND MAGMA_EXAMPLES_SRCS
batched_mass_bench.cpp
)
# Include the source directory where mfem.hpp and mfem-performance.hpp are.
include_directories(BEFORE ${PROJECT_BINARY_DIR})
# Add "test_magma" target (build-only by default).
add_custom_target(test_magma
COMMENT "Building MAGMA examples ...")
# Add one executable per cpp file, adding "magma_" as prefix to avoid
# collisions in a unified build tree.
set(PFX magma_)
add_mfem_examples(MAGMA_EXAMPLES_SRCS ${PFX} "" test_magma)
+12
View File
@@ -0,0 +1,12 @@
This directory contains an MFEM example focused on benchmarking MAGMA batched
element-mass-matrix inverses/solves on GPU backends (CUDA/HIP).
- `batched_mass_bench.cpp`: Assemble packed element mass matrices and compare:
- MFEM tripack inverse apply (packed lower-triangular)
- MAGMA packed Cholesky factor + solve
- MAGMA full (dense) batched Cholesky factor + solve
Build and run from this directory (after building MFEM with MAGMA enabled):
- `make batched_mass_bench`
- `./batched_mass_bench -m ../../data/beam-hex.mesh -o 3 -r 2 -d hip -reps 100`
+798
View File
@@ -0,0 +1,798 @@
// MFEM Batched Mass Benchmark
//
// Compile with: make batched_mass_bench
//
// Sample runs: batched_mass_bench -m ../../data/beam-hex.mesh -o 3 -r 2 -d hip -reps 100
// batched_mass_bench -m ../../data/beam-hex.mesh -o 3 -r 3 -d cuda -reps 50
//
// Description: This example isolates the element mass-matrix
// inverse paths used by element matrix kernels. It
// assembles packed lower-triangular L2 element mass matrices,
// then times repeated mass inverse applications using MFEM's
// tripack path and MAGMA (packed and full-matrix variants).
#include "mfem.hpp"
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <memory>
using namespace mfem;
using namespace std;
namespace
{
#ifdef MFEM_USE_MAGMA
inline void MagmaSetPointer(float **output, float *input,
magma_int_t inc, magma_int_t row, magma_int_t col,
magma_int_t lda, magma_int_t batch_count,
magma_queue_t queue)
{
magma_sset_pointer(output, input, inc, row, col, lda, batch_count, queue);
}
inline void MagmaSetPointer(double **output, double *input,
magma_int_t inc, magma_int_t row, magma_int_t col,
magma_int_t lda, magma_int_t batch_count,
magma_queue_t queue)
{
magma_dset_pointer(output, input, inc, row, col, lda, batch_count, queue);
}
inline magma_int_t MagmaPotrfBatched(magma_uplo_t uplo, magma_int_t n,
float **dA, magma_int_t ldda,
magma_int_t *info, magma_int_t batch_count,
magma_queue_t queue)
{
return magma_spotrf_batched(uplo, n, dA, ldda, info, batch_count, queue);
}
inline magma_int_t MagmaPotrfBatched(magma_uplo_t uplo, magma_int_t n,
double **dA, magma_int_t ldda,
magma_int_t *info, magma_int_t batch_count,
magma_queue_t queue)
{
return magma_dpotrf_batched(uplo, n, dA, ldda, info, batch_count, queue);
}
inline magma_int_t MagmaPotrsBatched(magma_uplo_t uplo, magma_int_t n,
magma_int_t nrhs, float **dA,
magma_int_t ldda, float **dB,
magma_int_t lddb, magma_int_t batch_count,
magma_queue_t queue)
{
return magma_spotrs_batched(uplo, n, nrhs, dA, ldda, dB, lddb,
batch_count, queue);
}
inline magma_int_t MagmaPotrsBatched(magma_uplo_t uplo, magma_int_t n,
magma_int_t nrhs, double **dA,
magma_int_t ldda, double **dB,
magma_int_t lddb, magma_int_t batch_count,
magma_queue_t queue)
{
return magma_dpotrs_batched(uplo, n, nrhs, dA, ldda, dB, lddb,
batch_count, queue);
}
void PrintMagmaFasterCondition(const double eq_fixed_ms,
const double eq_apply_ms,
const double magma_fixed_ms,
const double magma_apply_ms)
{
cout << "MAGMA faster condition (assembly+setup+N applies): ";
if (magma_fixed_ms <= eq_fixed_ms && magma_apply_ms <= eq_apply_ms)
{
if (magma_fixed_ms == eq_fixed_ms && magma_apply_ms == eq_apply_ms)
{
cout << "modeled totals are equal for all positive integer N.\n";
}
else
{
cout << "faster for every positive integer N.\n";
}
return;
}
if (magma_fixed_ms >= eq_fixed_ms && magma_apply_ms >= eq_apply_ms)
{
cout << "not faster for any positive integer N.\n";
return;
}
if (magma_apply_ms > eq_apply_ms)
{
const double crossover =
(eq_fixed_ms - magma_fixed_ms)/(magma_apply_ms - eq_apply_ms);
if (crossover <= 1.0)
{
cout << "only for N < " << crossover
<< ", so not for any positive integer N.\n";
}
else
{
const double last_n =
floor(nextafter(crossover, -numeric_limits<double>::infinity()));
cout << "faster for N < " << crossover
<< " applies (positive integer N <= " << (long long)last_n
<< "); eq-iter is faster above that.\n";
}
}
else
{
const double crossover =
(magma_fixed_ms - eq_fixed_ms)/(eq_apply_ms - magma_apply_ms);
const double first_n = floor(crossover) + 1.0;
cout << "faster for N > " << crossover
<< " applies (positive integer N >= " << (long long)first_n
<< "); eq-iter is faster below that.\n";
}
}
real_t **SetMagmaPackedPointerArray(Array<real_t *> &ptrs, real_t *data,
const int stride,
const int batch_size,
const magma_queue_t queue)
{
if (ptrs.Size() != batch_size)
{
if (ptrs.Size() != 0) { magma_queue_sync(queue); }
ptrs.SetSize(batch_size, Device::GetDeviceMemoryType());
}
real_t **d_ptrs = ptrs.Write();
MagmaSetPointer(d_ptrs, data, 1, 0, 0, stride, batch_size, queue);
return d_ptrs;
}
void ComputeMagmaFullCholeskyLower(
const Vector &full,
const int n,
Vector &factor,
Array<real_t *> &factor_ptrs,
const magma_queue_t queue)
{
const int batch_size = full.Size()/(n*n);
MFEM_VERIFY(full.Size() == batch_size*n*n, "Invalid full matrix storage.");
factor = full;
if (batch_size == 0) { return; }
real_t *factor_data = factor.ReadWrite();
real_t **dA = SetMagmaPackedPointerArray(factor_ptrs, factor_data, n*n,
batch_size, queue);
Array<magma_int_t> info_array(batch_size, Device::GetDeviceMemoryType());
magma_int_t *d_info = info_array.Write();
magma_memset(d_info, 0, batch_size*sizeof(magma_int_t));
const magma_int_t status =
MagmaPotrfBatched(MagmaLower, n, dA, n, d_info, batch_size, queue);
MFEM_VERIFY(status == MAGMA_SUCCESS, "MAGMA full potrf batched failed.");
magma_queue_sync(queue);
const magma_int_t *info = info_array.HostRead();
for (int e = 0; e < batch_size; ++e)
{
MFEM_VERIFY(info[e] == 0, "MAGMA full potrf failed on matrix " << e << '.');
}
}
void SolveMagmaFullCholeskyLowerInPlace(
const Vector &full_factor,
const int n,
const Array<real_t *> &factor_ptrs,
Array<real_t *> &rhs_ptrs,
Vector &rhs_sol,
const magma_queue_t queue)
{
const int batch_size = rhs_sol.Size()/n;
MFEM_VERIFY(rhs_sol.Size() == batch_size*n, "Invalid RHS size.");
MFEM_VERIFY(factor_ptrs.Size() == batch_size,
"Factor pointer array has the wrong size.");
if (batch_size == 0) { return; }
real_t *factor_data = const_cast<real_t *>(full_factor.Read());
real_t **dA = const_cast<real_t **>(factor_ptrs.Read());
(void)factor_data; // Factor data is referenced by dA (for clarity).
real_t *rhs_data = rhs_sol.ReadWrite();
real_t **dB = SetMagmaPackedPointerArray(rhs_ptrs, rhs_data, n, batch_size,
queue);
const magma_int_t status =
MagmaPotrsBatched(MagmaLower, n, 1, dA, n, dB, n, batch_size, queue);
MFEM_VERIFY(status == MAGMA_SUCCESS, "MAGMA full potrs batched failed.");
}
#endif
void FillRHS(Vector &rhs)
{
rhs.UseDevice(true);
real_t *x = rhs.HostWrite();
for (int i = 0; i < rhs.Size(); ++i)
{
x[i] = 1.0 + real_t((13*i + 7) % 29)/real_t(29);
}
}
void ComputeLowerPackedResidual(
const TriPackLowerMatrix &lower,
const Vector &x,
const Vector &rhs,
double &l2_residual,
double &relative_l2_residual,
real_t &max_residual,
real_t &relative_max_residual)
{
const int n = lower.GetNumRows();
const int batch_size = lower.GetNumMatrices();
const int packed_size = lower.GetPackedSize();
MFEM_VERIFY(x.Size() == batch_size*n, "Solution vector has the wrong size.");
MFEM_VERIFY(rhs.Size() == batch_size*n, "Right-hand side has the wrong size.");
const real_t *A = lower.Data().HostRead();
const real_t *X = x.HostRead();
const real_t *B = rhs.HostRead();
long double l2_sum = 0.0;
long double rhs_l2_sum = 0.0;
real_t max_abs = 0.0;
real_t rhs_max_abs = 0.0;
for (int e = 0; e < batch_size; ++e)
{
const real_t *Ae = A + e*packed_size;
const real_t *Xe = X + e*n;
const real_t *Be = B + e*n;
for (int i = 0; i < n; ++i)
{
long double ax = 0.0;
for (int j = 0; j < n; ++j)
{
const real_t aij =
(i >= j) ?
Ae[TriPackLowerMatrix::LowerIndex(i, j, n)] :
Ae[TriPackLowerMatrix::LowerIndex(j, i, n)];
ax += (long double)aij * (long double)Xe[j];
}
const long double residual = ax - (long double)Be[i];
l2_sum += residual*residual;
rhs_l2_sum += (long double)Be[i]*(long double)Be[i];
max_abs = max(max_abs, (real_t)fabs((double)residual));
rhs_max_abs = max(rhs_max_abs, fabs(Be[i]));
}
}
l2_residual = sqrt((double)l2_sum);
const double rhs_l2_norm = sqrt((double)rhs_l2_sum);
relative_l2_residual =
(rhs_l2_norm > 0.0) ? l2_residual/rhs_l2_norm : l2_residual;
max_residual = max_abs;
relative_max_residual =
(rhs_max_abs > 0.0) ? max_residual/rhs_max_abs : max_residual;
}
void ApplyLowerInverseInPlace(
const TriPackLowerMatrix &lower_inverse,
Vector &x,
Vector &work)
{
const int n = lower_inverse.GetNumRows();
const int batch_size = lower_inverse.GetNumMatrices();
const int packed_size = lower_inverse.GetPackedSize();
MFEM_VERIFY(x.Size() == batch_size*n, "Input vector has the wrong size.");
work.SetSize(batch_size*n);
work.UseDevice(true);
const real_t *L = lower_inverse.Data().Read();
const real_t *X = x.Read();
real_t *T = work.Write();
mfem::forall(batch_size*n, [=] MFEM_HOST_DEVICE (int idx)
{
const int i = idx % n;
const int e = idx / n;
const real_t *Le = L + e*packed_size;
const real_t *Xe = X + e*n;
real_t sum = 0.0;
for (int j = 0; j <= i; ++j)
{
sum += Le[TriPackLowerMatrix::LowerIndex(i, j, n)]*Xe[j];
}
T[idx] = sum;
});
const real_t *T_in = work.Read();
real_t *Y = x.Write();
mfem::forall(batch_size*n, [=] MFEM_HOST_DEVICE (int idx)
{
const int i = idx % n;
const int e = idx / n;
const real_t *Le = L + e*packed_size;
const real_t *Te = T_in + e*n;
real_t sum = 0.0;
for (int j = i; j < n; ++j)
{
sum += Le[TriPackLowerMatrix::LowerIndex(j, i, n)]*Te[j];
}
Y[idx] = sum;
});
}
double TimeLowerInverseApply(
const TriPackLowerMatrix &inverse,
const Vector &rhs,
const int reps,
Vector &x,
Vector &work)
{
StopWatch sw;
// Dry run to remove first-use kernel and workspace allocation costs.
x = rhs;
ApplyLowerInverseInPlace(inverse, x, work);
MFEM_DEVICE_SYNC;
sw.Start();
for (int r = 0; r < reps; ++r)
{
x = rhs;
ApplyLowerInverseInPlace(inverse, x, work);
}
MFEM_DEVICE_SYNC;
sw.Stop();
return 1000.0*sw.RealTime()/reps;
}
#ifdef MFEM_USE_MAGMA
double TimeMagmaSolve(
const TriPackLowerMatrix &lower_factor,
const Vector &rhs,
const int reps,
Vector &x,
MagmaPackedLowerCholesky &ws)
{
StopWatch sw;
// Dry run to remove first-use MAGMA and RHS pointer-array setup costs.
x = rhs;
ws.SolveInPlace(lower_factor, x);
MFEM_DEVICE_SYNC;
sw.Start();
for (int r = 0; r < reps; ++r)
{
x = rhs;
ws.SolveInPlace(lower_factor, x);
}
MFEM_DEVICE_SYNC;
sw.Stop();
return 1000.0*sw.RealTime()/reps;
}
double TimeMagmaInverseApply(
const TriPackLowerMatrix &lower_inverse,
const Vector &rhs,
const int reps,
Vector &x,
MagmaPackedLowerInverse &ws)
{
StopWatch sw;
// Dry run to remove first-use MAGMA and RHS pointer-array setup costs.
x = rhs;
ws.ApplyInPlace(lower_inverse, x);
MFEM_DEVICE_SYNC;
sw.Start();
for (int r = 0; r < reps; ++r)
{
x = rhs;
ws.ApplyInPlace(lower_inverse, x);
}
MFEM_DEVICE_SYNC;
sw.Stop();
return 1000.0*sw.RealTime()/reps;
}
double TimeMagmaFullSolve(
const Vector &full_factor,
const int n,
const Array<real_t *> &factor_ptrs,
const Vector &rhs,
const int reps,
Vector &x,
const magma_queue_t queue)
{
StopWatch sw;
Array<real_t *> rhs_ptrs;
x = rhs;
SolveMagmaFullCholeskyLowerInPlace(full_factor, n, factor_ptrs, rhs_ptrs, x,
queue);
MFEM_DEVICE_SYNC;
sw.Start();
for (int r = 0; r < reps; ++r)
{
x = rhs;
SolveMagmaFullCholeskyLowerInPlace(full_factor, n, factor_ptrs, rhs_ptrs,
x, queue);
}
MFEM_DEVICE_SYNC;
sw.Stop();
return 1000.0*sw.RealTime()/reps;
}
#endif
} // namespace
int main(int argc, char *argv[])
{
const char *mesh_file = "../../data/beam-hex.mesh";
int order = 3;
int ref_levels = 1;
int reps = 100;
int setup_reps = 10;
const char *device_config = "cpu";
bool use_magma = true;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
args.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree).");
args.AddOption(&ref_levels, "-r", "--refine",
"Number of uniform refinements.");
args.AddOption(&reps, "-reps", "--apply-repetitions",
"Number of mass inverse applications to time.");
args.AddOption(&setup_reps, "-sreps", "--setup-repetitions",
"Number of setup repetitions to time.");
args.AddOption(&device_config, "-d", "--device",
"Device configuration string, see Device::Configure().");
args.AddOption(&use_magma, "-magma", "--magma-solve",
"-no-magma", "--no-magma-solve",
"Time the MAGMA packed Cholesky solve when available.");
args.ParseCheck();
MFEM_VERIFY(reps > 0, "The number of apply repetitions must be positive.");
MFEM_VERIFY(setup_reps > 0,
"The number of setup repetitions must be positive.");
Device device(device_config);
device.Print();
#ifndef MFEM_USE_MAGMA
MFEM_VERIFY(!use_magma, "MFEM was built without MAGMA support.");
#endif
Mesh mesh(mesh_file, 1, 1);
for (int l = 0; l < ref_levels; ++l)
{
mesh.UniformRefinement();
}
const int dim = mesh.Dimension();
L2_FECollection fec(order, dim, BasisType::Positive);
FiniteElementSpace fespace(&mesh, &fec);
MFEM_VERIFY(UsesTensorBasis(fespace),
"This example requires a tensor-product finite element space.");
const int ne = mesh.GetNE();
const int elem_dofs = fespace.GetTypicalFE()->GetDof();
MassIntegrator mass;
StopWatch sw;
// Dry run assembly before timing steady-state work.
TriPackLowerMatrix packed_ea;
mass.AssembleEATriangular(fespace, packed_ea, false);
MFEM_DEVICE_SYNC;
sw.Clear();
sw.Start();
mass.AssembleEATriangular(fespace, packed_ea, false);
MFEM_DEVICE_SYNC;
sw.Stop();
const double assemble_ms = 1000.0*sw.RealTime();
TriPackLowerMatrix tripack_inverse;
// Dry run setup before timing steady-state setup work.
tripack::ComputeCholeskyLowerInverse(packed_ea, tripack_inverse);
MFEM_DEVICE_SYNC;
sw.Clear();
sw.Start();
for (int r = 0; r < setup_reps; ++r)
{
tripack::ComputeCholeskyLowerInverse(packed_ea, tripack_inverse);
}
MFEM_DEVICE_SYNC;
sw.Stop();
const double tripack_inverse_setup_ms = 1000.0*sw.RealTime()/setup_reps;
#ifdef MFEM_USE_MAGMA
magma_queue_t magma_queue = nullptr;
if (use_magma) { magma_queue = Magma::Queue(); }
TriPackLowerMatrix magma_factor;
std::unique_ptr<MagmaPackedLowerCholesky> magma_chol_ws;
double magma_factor_ms = 0.0;
TriPackLowerMatrix magma_inverse;
std::unique_ptr<MagmaPackedLowerInverse> magma_inv_ws;
double magma_inverse_ms = 0.0;
bool magma_ppinv_enabled = false;
Vector magma_full_factor;
Array<real_t *> magma_full_factor_ptrs;
double magma_full_factor_ms = 0.0;
if (use_magma)
{
magma_chol_ws.reset(new MagmaPackedLowerCholesky());
magma_chol_ws->SetQueue(magma_queue);
// Dry run setup before timing steady-state setup work.
magma_chol_ws->Factor(packed_ea, magma_factor);
MFEM_DEVICE_SYNC;
sw.Clear();
sw.Start();
for (int r = 0; r < setup_reps; ++r)
{
magma_chol_ws->Factor(packed_ea, magma_factor);
}
MFEM_DEVICE_SYNC;
sw.Stop();
magma_factor_ms = 1000.0*sw.RealTime()/setup_reps;
// Benchmark packed inverse (ppinv) only for sizes supported by MAGMA's
// current packed-inverse apply kernel.
if (elem_dofs <= 64)
{
magma_ppinv_enabled = true;
magma_inv_ws.reset(new MagmaPackedLowerInverse());
magma_inv_ws->SetQueue(magma_queue);
// Dry run setup before timing steady-state setup work.
magma_inv_ws->Compute(packed_ea, magma_inverse);
MFEM_DEVICE_SYNC;
sw.Clear();
sw.Start();
for (int r = 0; r < setup_reps; ++r)
{
magma_inv_ws->Compute(packed_ea, magma_inverse);
}
MFEM_DEVICE_SYNC;
sw.Stop();
magma_inverse_ms = 1000.0*sw.RealTime()/setup_reps;
}
// Full (dense) batched Cholesky factorization for comparison.
Vector full_ea(ne*elem_dofs*elem_dofs);
full_ea.UseDevice(true);
mass.AssembleEA(fespace, full_ea, false);
magma_full_factor.SetSize(full_ea.Size(), Device::GetDeviceMemoryType());
magma_full_factor.UseDevice(true);
ComputeMagmaFullCholeskyLower(full_ea, elem_dofs, magma_full_factor,
magma_full_factor_ptrs, magma_queue);
MFEM_DEVICE_SYNC;
sw.Clear();
sw.Start();
for (int r = 0; r < setup_reps; ++r)
{
ComputeMagmaFullCholeskyLower(full_ea, elem_dofs, magma_full_factor,
magma_full_factor_ptrs, magma_queue);
}
MFEM_DEVICE_SYNC;
sw.Stop();
magma_full_factor_ms = 1000.0*sw.RealTime()/setup_reps;
}
#endif
Vector rhs(ne*elem_dofs);
FillRHS(rhs);
Vector tripack_x(rhs.Size()), work;
tripack_x.UseDevice(true);
const double tripack_apply_ms =
TimeLowerInverseApply(tripack_inverse, rhs, reps, tripack_x, work);
double tripack_res_l2 = 0.0, tripack_rel_res_l2 = 0.0;
real_t tripack_res_max = 0.0, tripack_rel_res_max = 0.0;
ComputeLowerPackedResidual(packed_ea, tripack_x, rhs,
tripack_res_l2, tripack_rel_res_l2,
tripack_res_max, tripack_rel_res_max);
#ifdef MFEM_USE_MAGMA
double magma_solve_ms = 0.0;
double magma_full_solve_ms = 0.0;
double magma_res_l2 = 0.0, magma_rel_res_l2 = 0.0;
real_t magma_res_max = 0.0, magma_rel_res_max = 0.0;
double magma_full_res_l2 = 0.0, magma_full_rel_res_l2 = 0.0;
real_t magma_full_res_max = 0.0, magma_full_rel_res_max = 0.0;
Vector magma_x;
Vector magma_full_x;
if (use_magma)
{
magma_x.SetSize(rhs.Size());
magma_x.UseDevice(true);
magma_solve_ms =
TimeMagmaSolve(magma_factor, rhs, reps, magma_x, *magma_chol_ws);
ComputeLowerPackedResidual(packed_ea, magma_x, rhs,
magma_res_l2, magma_rel_res_l2,
magma_res_max, magma_rel_res_max);
magma_full_x.SetSize(rhs.Size());
magma_full_x.UseDevice(true);
magma_full_solve_ms =
TimeMagmaFullSolve(magma_full_factor, elem_dofs,
magma_full_factor_ptrs, rhs, reps, magma_full_x,
magma_queue);
ComputeLowerPackedResidual(packed_ea, magma_full_x, rhs,
magma_full_res_l2, magma_full_rel_res_l2,
magma_full_res_max, magma_full_rel_res_max);
}
double magma_ppinv_apply_ms = 0.0;
double magma_ppinv_res_l2 = 0.0, magma_ppinv_rel_res_l2 = 0.0;
real_t magma_ppinv_res_max = 0.0, magma_ppinv_rel_res_max = 0.0;
Vector magma_ppinv_x;
if (use_magma && magma_ppinv_enabled)
{
magma_ppinv_x.SetSize(rhs.Size());
magma_ppinv_x.UseDevice(true);
magma_ppinv_apply_ms =
TimeMagmaInverseApply(magma_inverse, rhs, reps, magma_ppinv_x,
*magma_inv_ws);
ComputeLowerPackedResidual(packed_ea, magma_ppinv_x, rhs,
magma_ppinv_res_l2, magma_ppinv_rel_res_l2,
magma_ppinv_res_max,
magma_ppinv_rel_res_max);
}
#endif
cout << fixed << setprecision(6);
cout << "Mass matrix inverse microbenchmark" << '\n';
cout << "Mesh: " << mesh_file << '\n';
cout << "Dimension: " << dim << '\n';
cout << "Elements: " << ne << '\n';
cout << "Element dofs: " << elem_dofs << '\n';
cout << "Scalar element unknowns: " << ne*elem_dofs << '\n';
cout << "Apply repetitions: " << reps << '\n';
cout << "Setup repetitions: " << setup_reps << '\n';
cout << '\n';
cout << "Assembly packed EA (ms): " << assemble_ms << '\n';
cout << "Setup MFEM tripack inverse (ms): " << tripack_inverse_setup_ms
<< '\n';
#ifdef MFEM_USE_MAGMA
if (use_magma)
{
cout << "Setup MAGMA packed Cholesky factor (ms): " << magma_factor_ms
<< '\n';
cout << "Setup MAGMA full Cholesky factor (ms): " << magma_full_factor_ms
<< '\n';
if (magma_ppinv_enabled)
{
cout << "Setup MAGMA packed inverse (ppinv) (ms): " << magma_inverse_ms
<< '\n';
}
else
{
cout << "Setup MAGMA packed inverse (ppinv) (ms): skipped "
<< "(requires element dofs <= 64)\n";
}
}
#endif
cout << '\n';
cout << "Apply MFEM tripack inverse (ms/apply): " << tripack_apply_ms
<< '\n';
#ifdef MFEM_USE_MAGMA
if (use_magma)
{
cout << "Apply MAGMA packed Cholesky solve (ms/apply): "
<< magma_solve_ms << '\n';
cout << "Apply MAGMA full Cholesky solve (ms/apply): "
<< magma_full_solve_ms << '\n';
cout << "MAGMA solve / MFEM tripack inverse apply: "
<< magma_solve_ms/tripack_apply_ms << '\n';
cout << "MAGMA full solve / MFEM tripack inverse apply: "
<< magma_full_solve_ms/tripack_apply_ms << '\n';
if (magma_ppinv_enabled)
{
cout << "Apply MAGMA packed inverse (ppinv) (ms/apply): "
<< magma_ppinv_apply_ms << '\n';
cout << "MAGMA ppinv apply / MFEM tripack inverse apply: "
<< magma_ppinv_apply_ms/tripack_apply_ms << '\n';
}
else
{
cout << "Apply MAGMA packed inverse (ppinv) (ms/apply): skipped "
<< "(requires element dofs <= 64)\n";
}
const double tripack_fixed_ms = assemble_ms + tripack_inverse_setup_ms;
const double magma_fixed_ms = assemble_ms + magma_factor_ms;
const double magma_full_fixed_ms = assemble_ms + magma_full_factor_ms;
const double tripack_total_ms =
tripack_fixed_ms + reps*tripack_apply_ms;
const double magma_total_ms = magma_fixed_ms + reps*magma_solve_ms;
const double magma_full_total_ms =
magma_full_fixed_ms + reps*magma_full_solve_ms;
cout << "Total MFEM tripack inverse for current repetitions "
<< "(assembly+setup+applies, ms): " << tripack_total_ms << '\n';
cout << "Total MAGMA packed Cholesky solve for current repetitions "
<< "(assembly+setup+applies, ms): " << magma_total_ms << '\n';
cout << "Total MAGMA full Cholesky solve for current repetitions "
<< "(assembly+setup+applies, ms): " << magma_full_total_ms << '\n';
cout << "Faster approach for current repetitions: "
<< ((magma_total_ms < tripack_total_ms &&
magma_total_ms <= magma_full_total_ms) ? "MAGMA packed" :
((magma_full_total_ms < tripack_total_ms &&
magma_full_total_ms < magma_total_ms) ? "MAGMA full" :
((tripack_total_ms < magma_total_ms &&
tripack_total_ms <= magma_full_total_ms) ? "tripack" :
"tie")))
<< '\n';
PrintMagmaFasterCondition(tripack_fixed_ms, tripack_apply_ms,
magma_fixed_ms, magma_solve_ms);
PrintMagmaFasterCondition(tripack_fixed_ms, tripack_apply_ms,
magma_full_fixed_ms, magma_full_solve_ms);
if (magma_ppinv_enabled)
{
const double magma_ppinv_fixed_ms =
assemble_ms + magma_inverse_ms;
const double magma_ppinv_total_ms =
magma_ppinv_fixed_ms + reps*magma_ppinv_apply_ms;
cout << "Total MAGMA packed inverse (ppinv) for current "
<< "repetitions (assembly+setup+applies, ms): "
<< magma_ppinv_total_ms << '\n';
}
}
#endif
cout << '\n';
cout << scientific << setprecision(12);
cout << "Residual, MFEM tripack inverse, max: "
<< tripack_res_max << " (relative "
<< tripack_rel_res_max << "), L2: "
<< tripack_res_l2 << " (relative "
<< tripack_rel_res_l2 << ")\n";
#ifdef MFEM_USE_MAGMA
if (use_magma)
{
cout << "Residual, MAGMA packed Cholesky solve, max: "
<< magma_res_max << " (relative "
<< magma_rel_res_max << "), L2: "
<< magma_res_l2 << " (relative "
<< magma_rel_res_l2 << ")\n";
cout << "Residual, MAGMA full Cholesky solve, max: "
<< magma_full_res_max << " (relative "
<< magma_full_rel_res_max << "), L2: "
<< magma_full_res_l2 << " (relative "
<< magma_full_rel_res_l2 << ")\n";
if (magma_ppinv_enabled)
{
cout << "Residual, MAGMA packed inverse (ppinv), max: "
<< magma_ppinv_res_max << " (relative "
<< magma_ppinv_rel_res_max << "), L2: "
<< magma_ppinv_res_l2 << " (relative "
<< magma_ppinv_rel_res_l2 << ")\n";
}
}
#endif
return 0;
}
+55
View File
@@ -0,0 +1,55 @@
# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
# LICENSE and NOTICE for details. LLNL-CODE-806117.
#
# This file is part of the MFEM library. For more information and source code
# availability visit https://mfem.org.
#
# MFEM is free software; you can redistribute it and/or modify it under the
# terms of the BSD-3 license. We welcome feedback and contributions, see file
# CONTRIBUTING.md for details.
# Use the MFEM build directory
MFEM_DIR ?= ../..
MFEM_BUILD_DIR ?= ../..
MFEM_INSTALL_DIR ?= ../../mfem
SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/examples/magma/,)
CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
$(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
MFEM_LIB_FILE = mfem_is_not_built
-include $(CONFIG_MK)
SEQ_EXAMPLES = batched_mass_bench
EXAMPLES = $(SEQ_EXAMPLES)
.SUFFIXES:
.SUFFIXES: .o .cpp .mk
.PHONY: all clean clean-build
# Remove built-in rule
%: %.cpp
# Replace the default implicit rule for *.cpp files
%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
all: $(EXAMPLES)
ifeq ($(MFEM_USE_MAGMA),NO)
$(EXAMPLES):
$(error MFEM is not configured with MAGMA)
endif
MFEM_TESTS = EXAMPLES
include $(MFEM_TEST_MK)
# Generate an error message if the MFEM library is not built and exit
$(MFEM_LIB_FILE):
$(error The MFEM library is not built)
clean: clean-build
clean-build:
rm -f *.o *~ $(SEQ_EXAMPLES)
rm -rf *.dSYM *.TVD.*breakpoints
+5 -1
View File
@@ -68,6 +68,9 @@ endif
ifeq ($(MFEM_USE_CALIPER),YES)
SUBDIRS += caliper
endif
ifeq ($(MFEM_USE_MAGMA),YES)
SUBDIRS += magma
endif
SUBDIRS_ALL = $(addsuffix /all,$(SUBDIRS))
SUBDIRS_TEST = $(addsuffix /test,$(SUBDIRS))
@@ -84,7 +87,8 @@ SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))
# Replace the default implicit rule for *.cpp files
%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@.o
$(MFEM_CXX) $(MFEM_LINK_FLAGS) $@.o -o $@ $(MFEM_LIBS)
all: $(EXAMPLES) $(SUBDIRS_ALL)
+7
View File
@@ -17,6 +17,7 @@
#include "fespace.hpp"
#include "ceed/interface/util.hpp"
#include "qfunction.hpp"
#include "../linalg/tripack.hpp"
#include <memory>
#include "kernel_dispatch.hpp"
@@ -2463,6 +2464,8 @@ protected:
int dim, ne, nq, dofs1D, quad1D;
void AssembleEA_(Vector &ea, const bool add);
void AssembleEATriangular_(TriPackLowerMatrix &ea,
const bool add);
public:
@@ -2515,6 +2518,10 @@ public:
void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
const bool add) override;
void AssembleEATriangular(const FiniteElementSpace &fes,
TriPackLowerMatrix &emat,
const bool add = false);
void AssembleEABoundary(const FiniteElementSpace &fes, Vector &emat,
const bool add) override;
-206
View File
@@ -22,8 +22,6 @@
#include <algorithm>
#include <cmath>
#include <cstdarg>
#include <unordered_map>
#include <unordered_set>
using namespace std;
@@ -4529,210 +4527,6 @@ void FiniteElementSpace
}
}
void FiniteElementSpace::GetBoundaryLoopEdgeDofs(
const Array<int> &boundary_element_indices,
Array<int> &boundary_edge_dofs,
Array<int> *dof_edges,
Array<int> *dof_boundary_elements) const
{
MFEM_VERIFY(mesh->Dimension() >= 2,
"GetBoundaryLoopEdgeDofs requires 2D or 3D meshes to find edge objects");
boundary_edge_dofs.SetSize(0);
if (dof_edges) { dof_edges->SetSize(0); }
if (dof_boundary_elements) { dof_boundary_elements->SetSize(0); }
// A DOF that appears in exactly one selected boundary element lies on the
// bounding loop; one appearing in two or more is interior to the boundary
// region and is dropped. Count occurrences of each DOF (using scratch maps,
// exposed only as parallel-indexed Array<int> below) and record, on first
// sight, the local edge and boundary element carrying it.
//
// The count is over GetEdgeDofs, which returns endpoint vertex DOFs as well
// as edge-interior DOFs (relevant for collections such as ND_R2D that carry
// vertex DOFs). Edge-interior DOFs occur once per edge, so the count mainly
// resolves vertex DOFs: a vertex shared by several elements is interior and
// dropped, while a genuine loop-corner (open-curve endpoint) vertex is kept.
// This is why we count GetEdgeDofs rather than collecting GetEdgeInteriorDofs,
// which would omit the endpoint vertex DOFs the method is documented to keep.
// The 3D removal criterion (any edge in two or more faces) matches the
// parallel version rather than a parity toggle.
std::unordered_map<int, int> dof_count, dof_edge, dof_belem;
Array<int> edge_dofs, edges, edge_orientations;
const int dim = mesh->Dimension();
for (int i = 0; i < boundary_element_indices.Size(); ++i)
{
const int boundary_element_idx = boundary_element_indices[i];
std::unordered_set<int> boundary_element_dofs;
if (dim == 3)
{
// Boundary elements are 2D faces; extract their 1D edges.
int face_index, face_orientation;
mesh->GetBdrElementFace(boundary_element_idx, &face_index,
&face_orientation);
mesh->GetFaceEdges(face_index, edges, edge_orientations);
}
else
{
// Boundary elements are 1D segments, each being a single edge.
mesh->GetBdrElementEdges(boundary_element_idx, edges, edge_orientations);
MFEM_VERIFY(edges.Size() == 1,
"2D boundary element should have exactly one edge");
}
for (int j = 0; j < edges.Size(); ++j)
{
GetEdgeDofs(edges[j], edge_dofs);
for (int k = 0; k < edge_dofs.Size(); ++k)
{
const int dof = edge_dofs[k];
// Count each DOF once per boundary element and record metadata the
// first time it is seen, so H1 DOFs shared by multiple edges of the
// same element are not double counted.
if (boundary_element_dofs.insert(dof).second &&
dof_count[dof]++ == 0)
{
dof_edge[dof] = edges[j];
dof_belem[dof] = boundary_element_idx;
}
}
}
}
// Emit the DOFs seen in exactly one selected boundary element, in a
// deterministic (increasing DOF index) order shared by all output arrays.
std::vector<int> kept;
kept.reserve(dof_count.size());
for (const auto &[dof, count] : dof_count)
{
if (count == 1) { kept.push_back(dof); }
}
std::sort(kept.begin(), kept.end());
boundary_edge_dofs.Reserve(static_cast<int>(kept.size()));
if (dof_edges) { dof_edges->Reserve(static_cast<int>(kept.size())); }
if (dof_boundary_elements)
{
dof_boundary_elements->Reserve(static_cast<int>(kept.size()));
}
for (int dof : kept)
{
boundary_edge_dofs.Append(dof);
if (dof_edges) { dof_edges->Append(dof_edge[dof]); }
if (dof_boundary_elements) { dof_boundary_elements->Append(dof_belem[dof]); }
}
}
void FiniteElementSpace::GetBoundaryElementsByAttribute(
const Array<int> &bdr_attrs,
std::vector<Array<int>> &attr_to_elements)
{
// One (initially empty) list of boundary elements per requested attribute,
// indexed to match bdr_attrs.
attr_to_elements.assign(bdr_attrs.Size(), Array<int>());
// Map attribute value -> position in bdr_attrs for quick lookup.
std::unordered_map<int, int> attr_to_index;
for (int i = 0; i < bdr_attrs.Size(); ++i)
{
attr_to_index[bdr_attrs[i]] = i;
}
// Bucket boundary elements by their attribute.
for (int i = 0; i < mesh->GetNBE(); ++i)
{
int attr = mesh->GetBdrElement(i)->GetAttribute();
auto it = attr_to_index.find(attr);
if (it != attr_to_index.end())
{
attr_to_elements[it->second].Append(i);
}
}
}
void FiniteElementSpace::GetBoundaryElementsByAttribute(int bdr_attr,
Array<int> &boundary_elements)
{
boundary_elements.SetSize(0);
for (int i = 0; i < mesh->GetNBE(); ++i)
{
if (mesh->GetBdrElement(i)->GetAttribute() == bdr_attr)
{
boundary_elements.Append(i);
}
}
}
void FiniteElementSpace::ComputeLoopEdgeOrientations(
const Array<int> &dof_edges,
const Array<int> &dof_boundary_elements,
const Vector &loop_normal,
Array<int> &dof_orientations) const
{
MFEM_VERIFY(dof_edges.Size() == dof_boundary_elements.Size(),
"dof_edges and dof_boundary_elements must be parallel-indexed");
const int ndof = dof_edges.Size();
dof_orientations.SetSize(ndof);
Array<int> edge_verts, bdr_elem_verts;
Vector edge_vec(3), to_edge_vec(3), cross_product(3);
for (int i = 0; i < ndof; i++)
{
const int edge_id = dof_edges[i];
const int bdr_elem_idx = dof_boundary_elements[i];
// Get edge vertices
mesh->GetEdgeVertices(edge_id, edge_verts);
const real_t *v0 = mesh->GetVertex(edge_verts[0]);
const real_t *v1 = mesh->GetVertex(edge_verts[1]);
// Get boundary element vertices
mesh->GetBdrElement(bdr_elem_idx)->GetVertices(bdr_elem_verts);
// Find the third vertex (not part of the edge)
int third_vertex = -1;
for (int j = 0; j < bdr_elem_verts.Size(); j++)
{
int v = bdr_elem_verts[j];
if (v != edge_verts[0] && v != edge_verts[1])
{
third_vertex = v;
break;
}
}
if (third_vertex == -1)
{
MFEM_ABORT("Boundary element " << bdr_elem_idx << " has only 2 vertices, "
"but 3D boundary elements must have at least 3 vertices");
}
const real_t *v2 = mesh->GetVertex(third_vertex);
// Edge vector
for (int j = 0; j < 3; j++) { edge_vec[j] = v1[j] - v0[j]; }
// Vector from third vertex to edge (use edge midpoint)
for (int j = 0; j < 3; j++)
{
real_t edge_midpoint = (v0[j] + v1[j]) * 0.5;
to_edge_vec[j] = edge_midpoint - v2[j];
}
// Cross product: to_edge × edge
to_edge_vec.cross3D(edge_vec, cross_product);
// Check alignment with loop normal
real_t dot_product = cross_product * loop_normal;
dof_orientations[i] = (dot_product > 0) ? 1 : -1;
}
}
FiniteElementCollection *FiniteElementSpace::Load(Mesh *m, std::istream &input)
{
string buff;
-75
View File
@@ -22,7 +22,6 @@
#include "restriction.hpp"
#include <iostream>
#include <unordered_map>
#include <vector>
namespace mfem
{
@@ -1390,80 +1389,6 @@ public:
virtual void GetExteriorTrueDofs(Array<int> &exterior_dofs,
int component = -1) const;
/** @brief Extract the edge degrees of freedom of a boundary "loop".
Here a "loop" is the set of boundary edges bounding the region covered by
@a boundary_element_indices: in 3D the outer edges of a patch of boundary
faces, in 2D the boundary segments themselves. An edge that is shared by
two (or more) of the selected boundary elements is interior to that region
rather than on its bounding loop, so its DOFs are excluded from the result.
This exclusion of interior DOFs is the defining feature of the method.
The three output arrays share a single indexing: for each valid index @a i,
@a dof_edges[i] and @a dof_boundary_elements[i] describe the DOF
@a boundary_edge_dofs[i].
@param[in] boundary_element_indices Boundary element indices spanning a
boundary surface (3D) or curve (2D).
@param[out] boundary_edge_dofs Local DOF indices on the boundary loop.
@param[out] dof_edges Optional; local edge index carrying each DOF.
@param[out] dof_boundary_elements Optional; a boundary element containing
each DOF.
@note In 3D the edge DOFs are extracted from the 1D edges of the 2D
boundary faces; in 2D they come directly from the 1D boundary segments, so
@a dof_edges then holds the boundary element (segment) edge indices.
@note This method uses GetEdgeDofs internally, which returns both vertex and
edge DOFs. Standard Nédélec elements (ND_FECollection) have no vertex DOFs,
so only genuine edge DOFs appear. Collections that carry vertex DOFs (e.g.
ND_R2D_FECollection) additionally contribute the vertex DOFs at loop
endpoints.
@note This is the serial version. For parallel meshes, use the parallel
version in ParFiniteElementSpace which handles processor boundaries
correctly.
@note Requires a 2D or 3D mesh to identify edge objects. The method will
assert if called on 1D meshes.
@note Only supports conforming meshes; non-conforming meshes are not
supported. */
void GetBoundaryLoopEdgeDofs(const Array<int> &boundary_element_indices,
Array<int> &boundary_edge_dofs,
Array<int> *dof_edges = nullptr,
Array<int> *dof_boundary_elements = nullptr) const;
/** @brief Get boundary elements grouped by attribute.
For each attribute in @a bdr_attrs, collect the indices of all boundary
elements carrying that attribute. The result is indexed to match
@a bdr_attrs: @a attr_to_elements[i] holds the boundary elements with
attribute @a bdr_attrs[i]. */
void GetBoundaryElementsByAttribute(
const Array<int> &bdr_attrs,
std::vector<Array<int>> &attr_to_elements);
/** @brief Get all boundary elements with a specific attribute. */
void GetBoundaryElementsByAttribute(int bdr_attr,
Array<int> &boundary_elements);
/** @brief Compute edge orientations relative to a boundary loop direction.
For each boundary-loop DOF described by @a dof_edges and
@a dof_boundary_elements (see GetBoundaryLoopEdgeDofs), determine whether
the carrying edge is
traversed in the direction consistent with @a loop_normal, following the
right-hand rule. Intended for 3D meshes.
@param[in] dof_edges Local edge index of each DOF (parallel-indexed with
the boundary_edge_dofs output of GetBoundaryLoopEdgeDofs).
@param[in] dof_boundary_elements A boundary element containing each DOF,
using the same indexing as @a dof_edges.
@param[in] loop_normal Normal vector defining the loop orientation.
@param[out] dof_orientations Orientation (+1 or -1) for each DOF, using the
same indexing as @a dof_edges. */
void ComputeLoopEdgeOrientations(const Array<int> &dof_edges,
const Array<int> &dof_boundary_elements,
const Vector &loop_normal,
Array<int> &dof_orientations) const;
/// Convert a Boolean marker array to a list containing all marked indices.
static void MarkerToList(const Array<int> &marker, Array<int> &list);
+104
View File
@@ -16,6 +16,71 @@
namespace mfem
{
namespace
{
void AssembleEATriangularImpl(const int dim, const int ne,
const int dofs1D, const int quad1D,
const Array<real_t> &B,
const Vector &pa_data,
Vector &data,
const bool add)
{
using internal::EAMassAssembleTriangular1DLower;
using internal::EAMassAssembleTriangular2DLower;
using internal::EAMassAssembleTriangular3DLower;
if (dim == 1)
{
auto kernel = EAMassAssembleTriangular1DLower<0,0>;
switch ((dofs1D << 4 ) | quad1D)
{
case 0x22: kernel = EAMassAssembleTriangular1DLower<2,2>; break;
case 0x33: kernel = EAMassAssembleTriangular1DLower<3,3>; break;
case 0x44: kernel = EAMassAssembleTriangular1DLower<4,4>; break;
case 0x55: kernel = EAMassAssembleTriangular1DLower<5,5>; break;
case 0x66: kernel = EAMassAssembleTriangular1DLower<6,6>; break;
case 0x77: kernel = EAMassAssembleTriangular1DLower<7,7>; break;
case 0x88: kernel = EAMassAssembleTriangular1DLower<8,8>; break;
case 0x99: kernel = EAMassAssembleTriangular1DLower<9,9>; break;
}
return kernel(ne, B, pa_data, data, add, dofs1D, quad1D);
}
else if (dim == 2)
{
auto kernel = EAMassAssembleTriangular2DLower<0,0>;
switch ((dofs1D << 4 ) | quad1D)
{
case 0x22: kernel = EAMassAssembleTriangular2DLower<2,2>; break;
case 0x33: kernel = EAMassAssembleTriangular2DLower<3,3>; break;
case 0x44: kernel = EAMassAssembleTriangular2DLower<4,4>; break;
case 0x55: kernel = EAMassAssembleTriangular2DLower<5,5>; break;
case 0x66: kernel = EAMassAssembleTriangular2DLower<6,6>; break;
case 0x77: kernel = EAMassAssembleTriangular2DLower<7,7>; break;
case 0x88: kernel = EAMassAssembleTriangular2DLower<8,8>; break;
case 0x99: kernel = EAMassAssembleTriangular2DLower<9,9>; break;
}
return kernel(ne, B, pa_data, data, add, dofs1D, quad1D);
}
else if (dim == 3)
{
auto kernel = EAMassAssembleTriangular3DLower<0,0>;
switch ((dofs1D << 4 ) | quad1D)
{
case 0x23: kernel = EAMassAssembleTriangular3DLower<2,3>; break;
case 0x34: kernel = EAMassAssembleTriangular3DLower<3,4>; break;
case 0x45: kernel = EAMassAssembleTriangular3DLower<4,5>; break;
case 0x56: kernel = EAMassAssembleTriangular3DLower<5,6>; break;
case 0x67: kernel = EAMassAssembleTriangular3DLower<6,7>; break;
case 0x78: kernel = EAMassAssembleTriangular3DLower<7,8>; break;
case 0x89: kernel = EAMassAssembleTriangular3DLower<8,9>; break;
}
return kernel(ne, B, pa_data, data, add, dofs1D, quad1D);
}
MFEM_ABORT("Unknown kernel.");
}
}
void MassIntegrator::AssembleEA_(Vector &ea_data,
const bool add)
@@ -75,6 +140,16 @@ void MassIntegrator::AssembleEA_(Vector &ea_data,
MFEM_ABORT("Unknown kernel.");
}
void MassIntegrator::AssembleEATriangular_(
TriPackLowerMatrix &ea_data,
const bool add)
{
Vector &data = ea_data.Data();
const Array<real_t> &B = maps->B;
return AssembleEATriangularImpl(dim, ne, dofs1D, quad1D, B, pa_data, data,
add);
}
void MassIntegrator::AssembleEA(const FiniteElementSpace &fes,
Vector &ea_data,
const bool add)
@@ -83,6 +158,35 @@ void MassIntegrator::AssembleEA(const FiniteElementSpace &fes,
if (ne > 0) { AssembleEA_(ea_data, add); }
}
void MassIntegrator::AssembleEATriangular(const FiniteElementSpace &fes,
TriPackLowerMatrix &ea_data,
const bool add)
{
AssemblePA(fes);
if (ne == 0) { return; }
int elem_dofs = 1;
for (int d = 0; d < dim; ++d)
{
elem_dofs *= dofs1D;
}
if (add)
{
MFEM_VERIFY(ea_data.GetNumRows() == elem_dofs,
"Invalid triangular EA element size.");
MFEM_VERIFY(ea_data.GetNumMatrices() == ne,
"Invalid triangular EA element count.");
}
else
{
ea_data.SetSize(elem_dofs, ne);
ea_data.UseDevice(true);
}
AssembleEATriangular_(ea_data, add);
}
void MassIntegrator::AssembleEABoundary(const FiniteElementSpace &fes,
Vector &ea_data,
const bool add)
+474
View File
@@ -1402,6 +1402,480 @@ inline void EAMassAssemble3D(const int NE,
});
}
template <int T_D1D = 0, int T_Q1D = 0>
inline void EAMassAssembleTriangular1DLower(const int NE,
const Array<real_t> &basis,
const Vector &padata,
Vector &eadata,
const bool add,
const int d1d = 0,
const int q1d = 0)
{
const int D1D = T_D1D ? T_D1D : d1d;
const int Q1D = T_Q1D ? T_Q1D : q1d;
MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
auto B = Reshape(basis.Read(), Q1D, D1D);
auto D = Reshape(padata.Read(), Q1D, NE);
auto M = Reshape(add ? eadata.ReadWrite() : eadata.Write(),
TriPackLowerMatrix::PackedSize(D1D), NE);
mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
{
const int D1D = T_D1D ? T_D1D : d1d;
const int Q1D = T_Q1D ? T_Q1D : q1d;
for (int i1 = 0; i1 < D1D; ++i1)
{
for (int j1 = 0; j1 <= i1; ++j1)
{
real_t val = 0.0;
for (int k1 = 0; k1 < Q1D; ++k1)
{
val += B(k1, i1) * B(k1, j1) * D(k1, e);
}
const int idx = TriPackLowerMatrix::LowerIndex(i1, j1, D1D);
if (add)
{
M(idx, e) += val;
}
else
{
M(idx, e) = val;
}
}
}
});
}
template <int T_D1D = 0, int T_Q1D = 0>
inline void EAMassAssembleTriangular2DLower(const int NE,
const Array<real_t> &basis,
const Vector &padata,
Vector &eadata,
const bool add,
const int d1d = 0,
const int q1d = 0)
{
const int D1D = T_D1D ? T_D1D : d1d;
const int Q1D = T_Q1D ? T_Q1D : q1d;
MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
const int ndofs = D1D*D1D;
auto B = Reshape(basis.Read(), Q1D, D1D);
auto D = Reshape(padata.Read(), Q1D, Q1D, NE);
auto M = Reshape(add ? eadata.ReadWrite() : eadata.Write(),
TriPackLowerMatrix::PackedSize(ndofs), NE);
mfem::forall_2D(NE, D1D, D1D, [=] MFEM_HOST_DEVICE (int e)
{
const int D1D = T_D1D ? T_D1D : d1d;
const int Q1D = T_Q1D ? T_Q1D : q1d;
constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
const int ndofs = D1D*D1D;
real_t r_B[MQ1][MD1];
for (int d = 0; d < D1D; d++)
{
for (int q = 0; q < Q1D; q++)
{
r_B[q][d] = B(q,d);
}
}
MFEM_SHARED real_t s_D[MQ1][MQ1];
MFEM_FOREACH_THREAD(k1,x,Q1D)
{
MFEM_FOREACH_THREAD(k2,y,Q1D)
{
s_D[k1][k2] = D(k1,k2,e);
}
}
MFEM_SYNC_THREAD;
MFEM_FOREACH_THREAD(i1,x,D1D)
{
MFEM_FOREACH_THREAD(i2,y,D1D)
{
const int row = i1 + D1D*i2;
for (int j2 = 0; j2 < D1D; ++j2)
{
for (int j1 = 0; j1 < D1D; ++j1)
{
const int col = j1 + D1D*j2;
if (row < col)
{
continue;
}
real_t val = 0.0;
for (int k1 = 0; k1 < Q1D; ++k1)
{
for (int k2 = 0; k2 < Q1D; ++k2)
{
val += r_B[k1][i1] * r_B[k1][j1]
* r_B[k2][i2] * r_B[k2][j2]
* s_D[k1][k2];
}
}
const int idx = TriPackLowerMatrix::LowerIndex(row, col, ndofs);
if (add)
{
M(idx, e) += val;
}
else
{
M(idx, e) = val;
}
}
}
}
}
});
}
template <int T_D1D, int T_Q1D, int T_COLB, int T_NT = 32>
inline void EAMassAssembleTriangular3D_LowerBlockCols_Impl(
const int NE,
const Array<real_t> &basis,
const Vector &padata,
Vector &eadata,
const bool add,
const int,
const int)
{
static_assert(T_D1D > 0 && T_Q1D > 0, "");
// Specialized packed lower-triangular hex mass assembly using block-column
// sum-factorization.
constexpr int D1D = T_D1D;
constexpr int Q1D = T_Q1D;
constexpr int COLB = T_COLB;
constexpr int NT = T_NT;
constexpr int ND = D1D*D1D*D1D;
constexpr int NQ = Q1D*Q1D*Q1D;
auto B = Reshape(basis.Read(), Q1D, D1D);
auto D = Reshape(padata.Read(), Q1D, Q1D, Q1D, NE);
auto M = Reshape(add ? eadata.ReadWrite() : eadata.Write(),
TriPackLowerMatrix::PackedSize(ND), NE);
mfem::forall_3D_grid(NE, NT, 1, 1, 0, [=] MFEM_HOST_DEVICE (int e)
{
const int tid = MFEM_THREAD_ID(x);
MFEM_SHARED real_t s_B[Q1D][D1D];
MFEM_SHARED real_t uW[NQ*COLB];
MFEM_SHARED real_t t1[D1D*Q1D*Q1D*COLB];
MFEM_SHARED real_t t2[D1D*D1D*Q1D*COLB];
for (int qb = tid; qb < Q1D*D1D; qb += NT)
{
const int q = qb % Q1D;
const int d = qb / Q1D;
s_B[q][d] = B(q, d);
}
MFEM_SYNC_THREAD;
for (int j0 = 0; j0 < ND; j0 += COLB)
{
const int b = (j0 + COLB <= ND) ? COLB : (ND - j0);
int j1[COLB], j2[COLB], j3[COLB];
for (int c = 0; c < COLB; ++c)
{
if (c < b)
{
const int jj = j0 + c;
j1[c] = jj % D1D;
const int tmp = jj / D1D;
j2[c] = tmp % D1D;
j3[c] = tmp / D1D;
}
}
for (int q = tid; q < NQ; q += NT)
{
const int q1 = q % Q1D;
const int tmp = q / Q1D;
const int q2 = tmp % Q1D;
const int q3 = tmp / Q1D;
const real_t Dq = D(q1, q2, q3, e);
for (int c = 0; c < b; ++c)
{
uW[q + NQ*c] = s_B[q1][j1[c]] * s_B[q2][j2[c]]
* s_B[q3][j3[c]] * Dq;
}
}
MFEM_SYNC_THREAD;
constexpr int T1S = D1D*Q1D*Q1D;
for (int a = tid; a < T1S; a += NT)
{
const int i1 = a % D1D;
const int tmp = a / D1D;
const int q2 = tmp % Q1D;
const int q3 = tmp / Q1D;
for (int c = 0; c < b; ++c)
{
real_t sum = 0.0;
for (int q1 = 0; q1 < Q1D; ++q1)
{
const int q = q1 + Q1D*(q2 + Q1D*q3);
sum += s_B[q1][i1] * uW[q + NQ*c];
}
t1[a + T1S*c] = sum;
}
}
MFEM_SYNC_THREAD;
constexpr int T2S = D1D*D1D*Q1D;
for (int a = tid; a < T2S; a += NT)
{
const int i1 = a % D1D;
const int tmp = a / D1D;
const int i2 = tmp % D1D;
const int q3 = tmp / D1D;
for (int c = 0; c < b; ++c)
{
real_t sum = 0.0;
for (int q2 = 0; q2 < Q1D; ++q2)
{
const int a1 = i1 + D1D*(q2 + Q1D*q3);
sum += s_B[q2][i2] * t1[a1 + T1S*c];
}
t2[a + T2S*c] = sum;
}
}
MFEM_SYNC_THREAD;
for (int c = 0; c < b; ++c)
{
const int col = j0 + c;
const int jj1 = j1[c];
const int jj2 = j2[c];
const int jj3 = j3[c];
for (int i3 = jj3 + 1; i3 < D1D; ++i3)
{
for (int a = tid; a < D1D*D1D; a += NT)
{
const int i1 = a % D1D;
const int i2 = a / D1D;
real_t sum = 0.0;
for (int q3 = 0; q3 < Q1D; ++q3)
{
const int a2 = i1 + D1D*(i2 + D1D*q3);
sum += s_B[q3][i3] * t2[a2 + T2S*c];
}
const int row = i1 + D1D*(i2 + D1D*i3);
const int idx =
TriPackLowerMatrix::LowerIndex(row, col, ND);
if (add) { M(idx, e) += sum; }
else { M(idx, e) = sum; }
}
MFEM_SYNC_THREAD;
}
for (int i2 = jj2 + 1; i2 < D1D; ++i2)
{
const int i3 = jj3;
for (int i1 = tid; i1 < D1D; i1 += NT)
{
real_t sum = 0.0;
for (int q3 = 0; q3 < Q1D; ++q3)
{
const int a2 = i1 + D1D*(i2 + D1D*q3);
sum += s_B[q3][i3] * t2[a2 + T2S*c];
}
const int row = i1 + D1D*(i2 + D1D*i3);
const int idx =
TriPackLowerMatrix::LowerIndex(row, col, ND);
if (add) { M(idx, e) += sum; }
else { M(idx, e) = sum; }
}
MFEM_SYNC_THREAD;
}
{
const int i3 = jj3;
const int i2 = jj2;
for (int i1 = tid + jj1; i1 < D1D; i1 += NT)
{
real_t sum = 0.0;
for (int q3 = 0; q3 < Q1D; ++q3)
{
const int a2 = i1 + D1D*(i2 + D1D*q3);
sum += s_B[q3][i3] * t2[a2 + T2S*c];
}
const int row = i1 + D1D*(i2 + D1D*i3);
const int idx =
TriPackLowerMatrix::LowerIndex(row, col, ND);
if (add) { M(idx, e) += sum; }
else { M(idx, e) = sum; }
}
MFEM_SYNC_THREAD;
}
}
}
});
}
template <int T_D1D, int T_Q1D, int T_NT = 32>
inline void EAMassAssembleTriangular3D_LowerBlockCols(
const int NE,
const Array<real_t> &basis,
const Vector &padata,
Vector &eadata,
const bool add,
const int d1d = 0,
const int q1d = 0)
{
static_assert(T_D1D > 0 && T_Q1D > 0, "");
constexpr int D1D = T_D1D;
constexpr int Q1D = T_Q1D;
constexpr int NQ = Q1D*Q1D*Q1D;
constexpr int SharedBytesPerCol =
sizeof(real_t)*(NQ + D1D*Q1D*Q1D + D1D*D1D*Q1D);
constexpr int SharedBytesBase = sizeof(real_t)*(Q1D*D1D);
constexpr int MaxSharedBytes = 48*1024;
constexpr int COLB =
(SharedBytesBase + 4*SharedBytesPerCol <= MaxSharedBytes) ? 4 :
(SharedBytesBase + 2*SharedBytesPerCol <= MaxSharedBytes) ? 2 : 1;
return EAMassAssembleTriangular3D_LowerBlockCols_Impl<T_D1D, T_Q1D, COLB, T_NT>(
NE, basis, padata, eadata, add, d1d, q1d);
}
template <int T_D1D = 0, int T_Q1D = 0>
inline void EAMassAssembleTriangular3DLower(const int NE,
const Array<real_t> &basis,
const Vector &padata,
Vector &eadata,
const bool add,
const int d1d = 0,
const int q1d = 0)
{
const int D1D = T_D1D ? T_D1D : d1d;
const int Q1D = T_Q1D ? T_Q1D : q1d;
MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
const int ndofs = D1D*D1D*D1D;
auto B = Reshape(basis.Read(), Q1D, D1D);
auto D = Reshape(padata.Read(), Q1D, Q1D, Q1D, NE);
auto M = Reshape(add ? eadata.ReadWrite() : eadata.Write(),
TriPackLowerMatrix::PackedSize(ndofs), NE);
if constexpr (T_D1D > 0 && T_Q1D > 0)
{
// Use the sum-factorized packed paths when the tensor dimensions are
// known at compile time. The generic path below handles dynamic sizes.
return EAMassAssembleTriangular3D_LowerBlockCols<T_D1D, T_Q1D>(
NE, basis, padata, eadata, add, d1d, q1d);
}
mfem::forall_3D(NE, D1D, D1D, D1D, [=] MFEM_HOST_DEVICE (int e)
{
const int D1D = T_D1D ? T_D1D : d1d;
const int Q1D = T_Q1D ? T_Q1D : q1d;
constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
constexpr int DQ = T_D1D * T_Q1D;
const int ndofs = D1D*D1D*D1D;
constexpr bool USE_REG = DQ != 0 && DQ <= 12;
constexpr int MD1r = USE_REG ? MD1 : 1;
constexpr int MQ1r = USE_REG ? MQ1 : 1;
constexpr int MD1s = USE_REG ? 1 : MD1;
constexpr int MQ1s = USE_REG ? 1 : MQ1;
MFEM_SHARED real_t s_B[MQ1s][MD1s];
real_t r_B[MQ1r][MD1r];
real_t (*l_B)[MD1] = nullptr;
if (USE_REG)
{
for (int d = 0; d < D1D; d++)
{
for (int q = 0; q < Q1D; q++)
{
r_B[q][d] = B(q,d);
}
}
l_B = (real_t (*)[MD1])r_B;
}
else
{
if (MFEM_THREAD_ID(z) == 0)
{
MFEM_FOREACH_THREAD(d,x,D1D)
{
MFEM_FOREACH_THREAD(q,y,Q1D)
{
s_B[q][d] = B(q,d);
}
}
}
l_B = (real_t (*)[MD1])s_B;
}
MFEM_SHARED real_t s_D[MQ1][MQ1][MQ1];
MFEM_FOREACH_THREAD(k1,x,Q1D)
{
MFEM_FOREACH_THREAD(k2,y,Q1D)
{
MFEM_FOREACH_THREAD(k3,z,Q1D)
{
s_D[k1][k2][k3] = D(k1,k2,k3,e);
}
}
}
MFEM_SYNC_THREAD;
MFEM_FOREACH_THREAD(i1,x,D1D)
{
MFEM_FOREACH_THREAD(i2,y,D1D)
{
MFEM_FOREACH_THREAD(i3,z,D1D)
{
const int row = i1 + D1D*(i2 + D1D*i3);
for (int j3 = 0; j3 < D1D; ++j3)
{
for (int j2 = 0; j2 < D1D; ++j2)
{
for (int j1 = 0; j1 < D1D; ++j1)
{
const int col = j1 + D1D*(j2 + D1D*j3);
if (row < col)
{
continue;
}
real_t val = 0.0;
for (int k1 = 0; k1 < Q1D; ++k1)
{
for (int k2 = 0; k2 < Q1D; ++k2)
{
for (int k3 = 0; k3 < Q1D; ++k3)
{
val += l_B[k1][i1] * l_B[k1][j1]
* l_B[k2][i2] * l_B[k2][j2]
* l_B[k3][i3] * l_B[k3][j3]
* s_D[k1][k2][k3];
}
}
}
const int idx = TriPackLowerMatrix::LowerIndex(row, col, ndofs);
if (add)
{
M(idx, e) += val;
}
else
{
M(idx, e) = val;
}
}
}
}
}
}
}
});
}
} // namespace internal
namespace
-338
View File
@@ -26,8 +26,6 @@
#include <limits>
#include <list>
#include <unordered_map>
#include <unordered_set>
namespace mfem
{
@@ -1287,342 +1285,6 @@ void ParFiniteElementSpace::GetExteriorVDofs(Array<int> &ext_dofs,
Synchronize(ext_dofs);
}
void ParFiniteElementSpace::GetBoundaryLoopEdgeDofs(
const Array<int> &boundary_element_indices,
Array<int> &ess_tdof_list,
Array<int> &boundary_edge_dofs_out,
Array<int> *ldof_marker,
Array<int> *dof_edges,
Array<int> *dof_boundary_elements,
Array<int> *ess_edge_list)
{
MFEM_VERIFY(!pmesh->Nonconforming(),
"GetBoundaryLoopEdgeDofs does not support nonconforming meshes");
MFEM_VERIFY(pmesh->Dimension() >= 2,
"GetBoundaryLoopEdgeDofs requires 2D or 3D meshes to find 1D edge objects");
// Call the serial version, then rebuild scratch maps/set from the returned
// arrays for the O(1) lookups the parallel reconciliation below needs.
Array<int> loc_dofs, loc_edges, loc_belems;
FiniteElementSpace::GetBoundaryLoopEdgeDofs(boundary_element_indices, loc_dofs,
&loc_edges, &loc_belems);
std::unordered_set<int> boundary_edge_dofs;
std::unordered_map<int, int> dof_to_edge_map;
std::unordered_map<int, int> dof_to_boundary_element;
boundary_edge_dofs.reserve(loc_dofs.Size());
dof_to_edge_map.reserve(loc_dofs.Size());
dof_to_boundary_element.reserve(loc_dofs.Size());
for (int i = 0; i < loc_dofs.Size(); i++)
{
const int dof = loc_dofs[i];
boundary_edge_dofs.insert(dof);
dof_to_edge_map[dof] = loc_edges[i];
dof_to_boundary_element[dof] = loc_belems[i];
}
// Parallel processing: Build edge sharing lookup table
std::unordered_map<int, int> edge_to_group_size;
int num_groups = pmesh->GetNGroups();
int total_shared_edges = 0;
for (int group = 1; group < num_groups; group++)
{
total_shared_edges += pmesh->GroupNEdges(group);
}
edge_to_group_size.reserve(total_shared_edges);
for (int group = 1; group < num_groups; group++)
{
int group_size = pmesh->gtopo.GetGroupSize(group);
int num_edges_in_group = pmesh->GroupNEdges(group);
for (int i = 0; i < num_edges_in_group; i++)
{
edge_to_group_size.emplace(pmesh->GroupEdge(group, i), group_size);
}
}
// Get global indices
Array<HYPRE_BigInt> global_edge_indices;
pmesh->GetGlobalEdgeIndices(global_edge_indices);
// Handle dimension-specific boundary element relationships
Array<HYPRE_BigInt> global_face_indices;
std::unordered_map<int, int> boundary_element_to_companion;
std::unordered_set<int> dofs_to_remove;
const int dim = pmesh->Dimension();
if (dim == 3)
{
// In 3D: boundary elements are faces, we track which face each boundary element is
pmesh->GetGlobalFaceIndices(global_face_indices);
for (int boundary_element_idx : boundary_element_indices)
{
int face_index, face_orientation;
pmesh->GetBdrElementFace(boundary_element_idx, &face_index, &face_orientation);
boundary_element_to_companion[boundary_element_idx] = face_index;
}
std::vector<HYPRE_BigInt> local_data;
local_data.reserve(boundary_edge_dofs.size() * 2);
std::unordered_set<int> processed_edges;
processed_edges.reserve(boundary_edge_dofs.size());
for (const auto& [dof, local_edge] : dof_to_edge_map)
{
// Skip if already processed this edge
if (!processed_edges.insert(local_edge).second) { continue; }
// Check if edge is shared (fast lookup)
auto it = edge_to_group_size.find(local_edge);
if (it != edge_to_group_size.end() && it->second > 1)
{
// Get boundary element and companion index directly from pre-computed map
int boundary_element_idx = dof_to_boundary_element[dof];
int companion_index = boundary_element_to_companion[boundary_element_idx];
// Store edge-face pair for 3D artificial boundary detection
local_data.push_back(global_edge_indices[local_edge]);
local_data.push_back(global_face_indices[companion_index]);
}
}
// MPI communication for 3D artificial boundary detection
int num_procs = pmesh->GetNRanks();
int local_size = local_data.size();
std::vector<int> mpi_arrays(num_procs * 4);
int* all_sizes = mpi_arrays.data();
int* displs = all_sizes + num_procs;
int* byte_sizes = displs + num_procs;
int* byte_displs = byte_sizes + num_procs;
MPI_Allgather(&local_size, 1, MPI_INT, all_sizes, 1, MPI_INT, pmesh->GetComm());
int total_size = 0;
constexpr int hypre_size = sizeof(HYPRE_BigInt);
for (int i = 0; i < num_procs; i++)
{
displs[i] = total_size;
byte_displs[i] = total_size * hypre_size;
total_size += all_sizes[i];
byte_sizes[i] = all_sizes[i] * hypre_size;
}
if (total_size > 0)
{
std::vector<HYPRE_BigInt> all_data(total_size);
MPI_Allgatherv(local_data.data(), local_size * hypre_size, MPI_BYTE,
all_data.data(), byte_sizes, byte_displs, MPI_BYTE, pmesh->GetComm());
// Build global-to-local edge mapping
std::unordered_map<HYPRE_BigInt, int> global_to_local_edge;
global_to_local_edge.reserve(global_edge_indices.Size());
for (int i = 0; i < global_edge_indices.Size(); ++i)
{
global_to_local_edge[global_edge_indices[i]] = i;
}
// Process collected data to find edges in multiple faces (artificial boundaries)
std::unordered_map<HYPRE_BigInt, std::unordered_set<HYPRE_BigInt>>edge_to_faces;
edge_to_faces.reserve(total_size / 2);
for (size_t i = 0; i < all_data.size(); i += 2)
{
edge_to_faces[all_data[i]].insert(all_data[i + 1]);
}
// Mark DOFs from artificial edges for removal
dofs_to_remove.reserve(local_data.size() / 4);
for (size_t i = 0; i < local_data.size(); i += 2)
{
HYPRE_BigInt global_edge_id = local_data[i];
// If this edge appears in 2+ distinct faces, it's artificial
if (edge_to_faces[global_edge_id].size() >= 2)
{
int local_edge = global_to_local_edge[global_edge_id];
Array<int> local_edge_dofs;
GetEdgeDofs(local_edge, local_edge_dofs);
// Mark boundary DOFs of this edge for removal
for (int k = 0; k < local_edge_dofs.Size(); ++k)
{
int dof = local_edge_dofs[k];
if (boundary_edge_dofs.count(dof))
{
dofs_to_remove.insert(dof);
}
}
}
}
}
}
else if (dim == 2)
{
// In 2D the boundary elements are themselves the edges, so there are no
// artificial boundary edges to detect. However, for collections with
// vertex DOFs (e.g. ND_R2D), a vertex shared by two boundary segments is
// interior to the boundary curve and must be dropped. The serial code
// does this by erasing a DOF on its second occurrence, which only sees
// the occurrences local to this rank. When the two segments meeting at a
// vertex live on different ranks, each rank sees a single occurrence and
// wrongly keeps the DOF. Reconcile the occurrence parity across each
// sharing group: membership in boundary_edge_dofs is the local parity,
// and the parities sum (mod 2) to the global occurrence parity.
Array<int> boundary_dof_count(GetVSize());
boundary_dof_count = 0;
for (const int dof : boundary_edge_dofs)
{
boundary_dof_count[dof] = 1;
}
// implement allreduce(+) as reduce(+) + broadcast
gcomm->Reduce<int>(boundary_dof_count, GroupCommunicator::Sum);
gcomm->Bcast(boundary_dof_count);
for (const int dof : boundary_edge_dofs)
{
if (boundary_dof_count[dof] % 2 == 0)
{
dofs_to_remove.insert(dof);
}
}
}
// Remove artificial DOFs
for (int dof : dofs_to_remove)
{
boundary_edge_dofs.erase(dof);
dof_to_edge_map.erase(dof);
dof_to_boundary_element.erase(dof);
}
// Convert to true DOFs and output
ess_tdof_list.SetSize(0);
ess_tdof_list.Reserve(boundary_edge_dofs.size());
if (ess_edge_list)
{
// Reset as well, so that it stays in correspondence with ess_tdof_list
// when the same output array is reused across calls.
ess_edge_list->SetSize(0);
ess_edge_list->Reserve(boundary_edge_dofs.size());
}
// Marker of the boundary edge DOFs. Always computed locally because the
// parallel reconciliation below needs it; only copied to the caller's output
// if requested (see the ldof_marker parameter).
Array<int> local_ldof_marker(GetVSize());
local_ldof_marker = 0;
for (int dof : boundary_edge_dofs)
{
local_ldof_marker[dof] = 1; // Mark all boundary edge dofs
}
// Make sure that a selected shared DOF is marked on every rank of its
// sharing group, including ranks holding none of the selected boundary
// elements. Only the group master owns the corresponding true DOF, so
// without this the true DOF would be emitted by no rank at all: the
// non-master ranks get -1 from GetLocalTDofNumber(), while the master may
// not have selected the DOF locally.
Synchronize(local_ldof_marker);
// A DOF marked only through the synchronization above has no local
// dof_to_edge_map entry, but the shared edge carrying it is still present in
// the local mesh. Build the missing DOF -> edge entries from the shared
// edges of the groups, so that ess_edge_list stays in correspondence with
// ess_tdof_list. Note that a vertex DOF is not associated with a unique
// edge, so it is only resolved when it is an interior DOF of an edge.
std::unordered_map<int, int> shared_dof_to_edge;
Array<int> shared_edge_dofs;
for (int group = 1; group < num_groups; group++)
{
const int num_edges_in_group = pmesh->GroupNEdges(group);
for (int i = 0; i < num_edges_in_group; i++)
{
const int edge = pmesh->GroupEdge(group, i);
GetEdgeInteriorDofs(edge, shared_edge_dofs);
for (int k = 0; k < shared_edge_dofs.Size(); k++)
{
shared_dof_to_edge.emplace(shared_edge_dofs[k], edge);
}
}
}
// Build parallel arrays for DOFs and corresponding edges
std::vector<std::pair<int, int>> tdof_edge_pairs;
tdof_edge_pairs.reserve(boundary_edge_dofs.size());
for (int dof = 0; dof < local_ldof_marker.Size(); dof++)
{
if (!local_ldof_marker[dof]) { continue; }
const int tdof = GetLocalTDofNumber(dof);
if (tdof < 0) { continue; } // tdof == -1 means not owned by this rank
int edge = -1;
auto it = dof_to_edge_map.find(dof);
if (it != dof_to_edge_map.end())
{
edge = it->second;
}
else
{
auto shared_it = shared_dof_to_edge.find(dof);
if (shared_it != shared_dof_to_edge.end())
{
edge = shared_it->second;
}
}
tdof_edge_pairs.push_back({tdof, edge});
}
// Sort by true DOF index to maintain consistent ordering
std::sort(tdof_edge_pairs.begin(), tdof_edge_pairs.end());
// Extract sorted true DOFs and edges
for (const auto& pair : tdof_edge_pairs)
{
ess_tdof_list.Append(pair.first);
if (ess_edge_list)
{
ess_edge_list->Append(pair.second);
}
}
// Emit the local boundary-loop DOFs in a deterministic (increasing DOF
// index) order shared by all output arrays.
std::vector<int> kept(boundary_edge_dofs.begin(), boundary_edge_dofs.end());
std::sort(kept.begin(), kept.end());
boundary_edge_dofs_out.SetSize(0);
boundary_edge_dofs_out.Reserve(static_cast<int>(kept.size()));
if (dof_edges)
{
dof_edges->SetSize(0);
dof_edges->Reserve(static_cast<int>(kept.size()));
}
if (dof_boundary_elements)
{
dof_boundary_elements->SetSize(0);
dof_boundary_elements->Reserve(static_cast<int>(kept.size()));
}
for (int dof : kept)
{
boundary_edge_dofs_out.Append(dof);
if (dof_edges) { dof_edges->Append(dof_to_edge_map[dof]); }
if (dof_boundary_elements)
{
dof_boundary_elements->Append(dof_to_boundary_element[dof]);
}
}
if (ldof_marker) { ldof_marker->Swap(local_ldof_marker); }
}
void ParFiniteElementSpace::GetExteriorTrueDofs(Array<int> &ext_tdof_list,
int component) const
{
-35
View File
@@ -460,41 +460,6 @@ public:
void GetExteriorTrueDofs(Array<int> &ext_tdof_list,
int component = -1) const override;
/** @brief Extract the edge degrees of freedom of a boundary "loop" on a
parallel mesh (see the serial FiniteElementSpace::GetBoundaryLoopEdgeDofs
for the definition of a loop). This version removes the artificial
boundary edges that appear at processor boundaries, so the selected DOFs
are independent of the mesh partitioning.
As in the serial version, the @a boundary_edge_dofs_out, @a dof_edges and
@a dof_boundary_elements outputs share a single indexing describing the
same local DOF at each position.
Requirements:
- Mesh must be conforming (no hanging nodes)
- Mesh dimension must be >= 2
@param[in] boundary_element_indices Array of boundary element indices.
@param[out] ess_tdof_list Essential true DOF indices, sorted ascending.
@param[out] boundary_edge_dofs_out Local boundary-loop DOF indices.
@param[out] ldof_marker Optional; marker of the boundary edge DOFs,
derivable from @a boundary_edge_dofs_out via ListToMarker().
@param[out] dof_edges Optional; local edge index of each DOF.
@param[out] dof_boundary_elements Optional; a boundary element containing
each DOF.
@param[out] ess_edge_list Optional array of edge indices, in one-to-one
correspondence with @a ess_tdof_list. An entry
is -1 when the true DOF is owned by this rank
but no local edge can be associated with it,
which can happen for a shared vertex DOF whose
boundary elements are all on other ranks. */
void GetBoundaryLoopEdgeDofs(const Array<int> &boundary_element_indices,
Array<int> &ess_tdof_list,
Array<int> &boundary_edge_dofs_out,
Array<int> *ldof_marker = nullptr,
Array<int> *dof_edges = nullptr,
Array<int> *dof_boundary_elements = nullptr,
Array<int> *ess_edge_list = nullptr);
/** If the given ldof is owned by the current processor, return its local
tdof number, otherwise return -1 */
int GetLocalTDofNumber(int ldof) const;
+173 -304
View File
@@ -231,11 +231,9 @@ const Operator &InterpolationGridTransfer::BackwardOperator()
L2ProjectionGridTransfer::L2Projection::L2Projection(
const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_,
CoefficientWithOrder coeff_ho_, CoefficientWithOrder coeff_lor_,
MemoryType d_mt_)
: Operator(fes_lor_.GetVSize(), fes_ho_.GetVSize()),
fes_ho(fes_ho_), fes_lor(fes_lor_), coeff_ho(coeff_ho_),
coeff_lor(coeff_lor_), d_mt(d_mt_)
fes_ho(fes_ho_), fes_lor(fes_lor_), d_mt(d_mt_)
{ }
void L2ProjectionGridTransfer::L2Projection::BuildHo2Lor(
@@ -265,13 +263,12 @@ void L2ProjectionGridTransfer::L2Projection::ElemMixedMass(
IntegrationPointTransformation& ip_tr,
DenseMatrix& M_mixed_el) const
{
int order = fe_lor.GetOrder() + fe_ho.GetOrder() + tr_lor->OrderW() +
coeff_ho.order;
const IntegrationRule &ir = IntRules.Get(geom, order);
int order = fe_lor.GetOrder() + fe_ho.GetOrder() + tr_lor->OrderW();
const IntegrationRule* ir = &IntRules.Get(geom, order);
M_mixed_el = 0.0;
for (int i = 0; i < ir.GetNPoints(); i++)
for (int i = 0; i < ir->GetNPoints(); i++)
{
const IntegrationPoint& ip_lor = ir.IntPoint(i);
const IntegrationPoint& ip_lor = ir->IntPoint(i);
IntegrationPoint ip_ho;
ip_tr.Transform(ip_lor, ip_ho);
Vector shape_lor(fe_lor.GetDof());
@@ -287,23 +284,23 @@ void L2ProjectionGridTransfer::L2Projection::ElemMixedMass(
{
w *= tr_lor->Weight();
}
if (coeff_ho)
{
w *= coeff_ho.coeff->Eval(*tr_ho, ip_ho);
}
shape_lor *= w;
AddMultVWt(shape_lor, shape_ho, M_mixed_el);
}
}
void L2ProjectionGridTransfer::L2Projection::ElemMixedEvaluation(
Geometry::Type geom, const FiniteElement& fe_ho, const FiniteElement& fe_lor,
IntegrationPointTransformation& ip_tr, const IntegrationRule& ir,
void L2ProjectionGridTransfer::L2Projection::ElemMixedMass(
Geometry::Type geom, const FiniteElement& fe_ho,
const FiniteElement& fe_lor, ElementTransformation* el_tr,
IntegrationPointTransformation& ip_tr,
DenseMatrix& B_L, DenseMatrix& B_H) const
{
for (int i = 0; i < ir.GetNPoints(); i++)
int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr->OrderW();
const IntegrationRule* ir = &IntRules.Get(geom, order);
for (int i = 0; i < ir->GetNPoints(); i++)
{
const IntegrationPoint& ip_lor = ir.IntPoint(i);
const IntegrationPoint& ip_lor = ir->IntPoint(i);
IntegrationPoint ip_ho;
// maps integration point ip_lor -> ip_ho
@@ -323,6 +320,7 @@ void L2ProjectionGridTransfer::L2Projection::ElemMixedEvaluation(
B_H(i, j) = shape_ho(j);
}
}
}
void L2ProjectionGridTransfer::L2Projection::MixedMassEA(
@@ -330,11 +328,10 @@ void L2ProjectionGridTransfer::L2Projection::MixedMassEA(
const FiniteElementSpace& fes_lor_ea,
Vector &M_LH, MemoryType d_mt_)
{
Mesh &mesh_ho = *fes_ho_ea.GetMesh();
Mesh &mesh_lor = *fes_lor_ea.GetMesh();
const int nel_ho = mesh_ho.GetNE();
const int nel_lor = mesh_lor.GetNE();
Mesh* mesh_ho = fes_ho_ea.GetMesh();
Mesh* mesh_lor = fes_lor_ea.GetMesh();
int nel_ho = mesh_ho->GetNE();
int nel_lor = mesh_lor->GetNE();
if (nel_ho == 0)
{
@@ -342,11 +339,11 @@ void L2ProjectionGridTransfer::L2Projection::MixedMassEA(
return;
}
const CoarseFineTransformations& cf_tr = mesh_lor.GetRefinementTransforms();
const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
int nref_max = 0;
Array<Geometry::Type> geoms;
mesh_ho.GetGeometries(mesh_ho.Dimension(), geoms);
mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
for (int ig = 0; ig < geoms.Size(); ++ig)
{
Geometry::Type geom = geoms[ig];
@@ -363,226 +360,130 @@ void L2ProjectionGridTransfer::L2Projection::MixedMassEA(
{
// Assume all HO elements are LOR in the same way
const int iho = 0;
Array<int> lor_els;
ho2lor.GetRow(iho, lor_els);
const int nref = ho2lor.RowSize(iho);
MFEM_VERIFY(nel_ho*nref == nel_lor, "we expect nel_ho*nref == nel_lor");
Geometry::Type geom = mesh_ho.GetElementBaseGeometry(iho);
emb_tr.SetIdentityTransformation(geom);
const DenseTensor &pmats = cf_tr.point_matrices[geom];
const FiniteElement &fe_ho = *fes_ho_ea.GetFE(iho);
const FiniteElement &fe_lor = *fes_lor_ea.GetFE(lor_els[0]);
// Allocate space for DenseTensors
ElementTransformation &el_tr = *mesh_lor.GetTypicalElementTransformation();
const int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr.OrderW()
+ coeff_ho.order;
const IntegrationRule &ir_ea = IntRules.Get(geom, order);
const int qPts = ir_ea.GetNPoints();
// Containers for the basis functions sampled at quadrature points
B_L.SetSize(qPts, fe_lor.GetDof(), nref, d_mt);
B_H.SetSize(qPts, fe_ho.GetDof(), nref, d_mt);
D.SetSize(qPts, nref, nel_ho, d_mt);
const GeometricFactors *geo_facts =
mesh_lor.GetGeometricFactors(ir_ea, GeometricFactors::DETERMINANTS);
Vector coeff_vec(qPts*nel_lor);
coeff_vec.UseDevice(true);
const int dim = mesh_ho.Dimension();
const int nq1d = (int)floor(pow(ir_ea.Size(), 1.0/dim) + 0.5);
const int nref_1d = (int)floor(pow(nref, 1.0/dim) + 0.5);
if (!coeff_ho)
{
coeff_vec = 1.0;
}
else if (UsesTensorBasis(fes_ho) &&
nq1d*nref_1d <= DeviceDofQuadLimits::Get().MAX_Q1D)
{
// Fast coefficient evaluation for tensor-product case. We create a
// "composite" quadrature rule in the high-order element that is the
// union of the quadrature rules within each of the low-order-refined
// subelements.
//
// NOTE: if the integration rule order is high and there are many LOR
// subelements, this can create a very big quadrature rule. That is
// why we need to check that we do not exceed MAX_Q1D. If we do, then
// we fall back on the slower "legacy" evaluation.
Array<int> lor_els;
ho2lor.GetRow(iho, lor_els);
int nref = ho2lor.RowSize(iho);
// Construct the composite rule as a tensor-product of the 1D LOR rule.
IntegrationRule ir_ho = [&]()
Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
const FiniteElement &fe_ho = *fes_ho_ea.GetFE(iho);
const FiniteElement &fe_lor = *fes_lor_ea.GetFE(lor_els[0]);
// Allocate space for DenseTensors
ElementTransformation *el_tr = fes_lor_ea.GetElementTransformation(0);
int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr->OrderW();
const IntegrationRule* ir_ea = &IntRules.Get(geom, order);
int qPts = ir_ea->GetNPoints();
// Containers for the basis functions sampled at quadrature points
B_L.SetSize(qPts, fe_lor.GetDof(), nref, d_mt);
B_H.SetSize(qPts, fe_ho.GetDof(), nref, d_mt);
D.SetSize(qPts, nref, nel_ho, d_mt);
const GeometricFactors *geo_facts =
mesh_lor->GetGeometricFactors(*ir_ea, GeometricFactors::DETERMINANTS);
MFEM_ASSERT(nel_ho*nref == nel_lor, "we expect nel_ho*nref == nel_lor");
// Setup data at quadrature points
// TODO add support for user coefficient
const auto W = Reshape(ir_ea->GetWeights().Read(), qPts);
const auto J = Reshape(geo_facts->detJ.Read(), qPts, nel_lor);
const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);
mfem::forall(qPts * nref * nel_ho, [=] MFEM_HOST_DEVICE (int tid)
{
IntegrationRule ir_ho_1d(nq1d * nref_1d);
for (int iref = 0; iref < nref_1d; ++iref)
{
const real_t a = pmats(cf_tr.embeddings[iref].matrix)(0,0);
const real_t b = pmats(cf_tr.embeddings[iref].matrix)(0,1);
for (int iq = 0; iq < nq1d; ++iq)
{
ir_ho_1d[iq + iref*nq1d].x = a + ir_ea[iq].x*(b - a);
}
}
if (dim == 1) { return ir_ho_1d; }
else if (dim == 2) { return IntegrationRule(ir_ho_1d, ir_ho_1d); }
else { return IntegrationRule(ir_ho_1d, ir_ho_1d, ir_ho_1d); }
}();
const int q = tid % qPts;
const int iref = (tid / qPts) % nref;
const int iho = (tid / (qPts * nref)) % nel_ho;
// Project the high-order coefficient on the high-order composite rule.
QuadratureSpace qs(mesh_ho, ir_ho);
CoefficientVector coeff_vec_ho(*coeff_ho.coeff, qs);
const int lo_el_id = iref + nref*iho;
const real_t detJ = J(q, lo_el_id);
// Permute the coefficient values to the expected LOR ordering.
const int nq_ho = ir_ho.Size();
const auto d_Q_ho = Reshape(coeff_vec_ho.Read(), nq_ho, nel_ho);
const auto d_Q = Reshape(coeff_vec.Write(), qPts, nel_lor);
d_D(q, iref, iho) = W(q) * detJ;
mfem::forall(nq_ho * nel_ho, [=] MFEM_HOST_DEVICE (int ii)
{
const int e_ho = ii / nq_ho;
const int iq_ho = ii % nq_ho;
int iq_tensor = iq_ho;
int iq_lor = 0;
int iref = 0;
int iq_stride = 1;
int iref_stride = 1;
const int nq_ho_1d = nq1d*nref_1d;
for (int d = 0; d < dim; ++d)
{
const int iq_ho_1d = iq_tensor % nq_ho_1d;
iq_tensor /= nq_ho_1d;
iq_lor += (iq_ho_1d % nq1d)*iq_stride;
iref += (iq_ho_1d / nq1d)*iref_stride;
iq_stride *= nq1d;
iref_stride *= nref_1d;
}
const int e_lor = iref + e_ho*nref;
d_Q(iq_lor, e_lor) = d_Q_ho(iq_ho, e_ho);
});
}
else
{
// Legacy/fallback coefficient evaluation for non-tensor-product cases
// or when the number of quadrature points is too large for the device
// kernels.
IntegrationPoint ip_ho;
for (int e_ho = 0; e_ho < nel_ho; ++e_ho)
emb_tr.SetIdentityTransformation(geom);
const DenseTensor &pmats = cf_tr.point_matrices[geom];
// Collect the basis functions
for (int iref = 0; iref < nref; ++iref)
{
ElementTransformation &ho_tr = *mesh_ho.GetElementTransformation(e_ho);
for (int iref = 0; iref < nref; ++iref)
{
const int e_lor = iref + e_ho*nref;
emb_tr.SetPointMat(pmats(cf_tr.embeddings[e_lor].matrix));
int ilor = lor_els[iref];
// Now assemble the block-row of the mixed mass matrix associated
// with integrating HO functions against LOR functions on the LOR
// sub-element.
for (int iq = 0; iq < qPts; ++iq)
{
const IntegrationPoint &ip_lor = ir_ea[iq];
ip_tr.Transform(ip_lor, ip_ho);
ho_tr.SetIntPoint(&ip_ho);
coeff_vec[iq + e_lor*qPts] = coeff_ho.coeff->Eval(ho_tr, ip_ho);
}
}
}
// Create the transformation that embeds the fine low-order element
// within the coarse high-order element in reference space
emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
DenseMatrix &b_lo = B_L(ilor);
DenseMatrix &b_ho = B_H(ilor);
ElemMixedMass(geom, fe_ho, fe_lor, el_tr, ip_tr, b_lo, b_ho);
} // loop over subcells of ho element
// end of quadrature point setup
}
// Setup data at quadrature points
const auto W = Reshape(ir_ea.GetWeights().Read(), qPts);
const auto J = Reshape(geo_facts->detJ.Read(), qPts, nel_lor);
const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);
const auto d_Q = Reshape(coeff_vec.Read(), qPts, nel_lor);
mfem::forall(qPts * nref * nel_ho, [=] MFEM_HOST_DEVICE (int tid)
{
const int q = tid % qPts;
const int iref = (tid / qPts) % nref;
const int iho = (tid / (qPts * nref)) % nel_ho;
const int lo_el_id = iref + nref*iho;
const real_t detJ = J(q, lo_el_id);
d_D(q, iref, iho) = W(q) * d_Q(q, lo_el_id) * detJ;
});
// Collect the basis functions
for (int iref = 0; iref < nref; ++iref)
{
int ilor = lor_els[iref];
// Now assemble the block-row of the mixed mass matrix associated
// with integrating HO functions against LOR functions on the LOR
// sub-element.
// Create the transformation that embeds the fine low-order element
// within the coarse high-order element in reference space
emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
DenseMatrix &b_lo = B_L(ilor);
DenseMatrix &b_ho = B_H(ilor);
ElemMixedEvaluation(geom, fe_ho, fe_lor, ip_tr, ir_ea, b_lo, b_ho);
} // loop over subcells of ho element
// end of quadrature point setup
} // completed setup of basis function and quadrature point
// Assemble mixed mass matrix
int iho = 0;
Array<int> lor_els;
ho2lor.GetRow(iho, lor_els);
int nref = ho2lor.RowSize(iho);
const FiniteElement &fe_ho = *fes_ho_ea.GetFE(iho);
const FiniteElement &fe_lor = *fes_lor_ea.GetFE(lor_els[0]);
const int ndof_ho = fe_ho.GetDof();
const int ndof_lor = fe_lor.GetDof();
const int qPts = D.SizeI();
M_LH.SetSize(ndof_lor*ndof_ho*nref*nel_ho, d_mt);
// Rows x columns
// Recall MFEM is column major
// rows x columns is inverted - matrix is ndof_lor x ndof_ho
auto v_M_LH = Reshape(M_LH.Write(), ndof_lor, ndof_ho, nref,
nel_ho);
const int fe_ho_ndof = fe_ho.GetDof();
const int fe_lor_ndof = fe_lor.GetDof();
auto d_B_L = Reshape(B_L.Read(), qPts, fe_lor_ndof, nref);
auto d_B_H = Reshape(B_H.Read(), qPts, fe_ho_ndof, nref);
auto d_D = Reshape(D.Read(), qPts, nref, nel_ho);
mfem::forall(fe_ho_ndof*nref*nel_ho, [=] MFEM_HOST_DEVICE (int idx)
{
const int bh = idx % fe_ho_ndof;
const int iref = (idx / fe_ho_ndof) % nref;
const int iho = idx / fe_ho_ndof / nref;
// (B_lo_dofs x Q) x (Q x B_ho_dofs)
for (int bl = 0; bl < fe_lor_ndof; ++bl)
int iho = 0;
Array<int> lor_els;
ho2lor.GetRow(iho, lor_els);
int nref = ho2lor.RowSize(iho);
const FiniteElement &fe_ho = *fes_ho_ea.GetFE(iho);
const FiniteElement &fe_lor = *fes_lor_ea.GetFE(lor_els[0]);
const int ndof_ho = fe_ho.GetDof();
const int ndof_lor = fe_lor.GetDof();
const int qPts = D.SizeI();
M_LH.SetSize(ndof_lor*ndof_ho*nref*nel_ho, d_mt);
// Rows x columns
// Recall MFEM is column major
// rows x columns is inverted - matrix is ndof_lor x ndof_ho
auto v_M_LH = Reshape(M_LH.Write(), ndof_lor, ndof_ho, nref,
nel_ho);
const int fe_ho_ndof = fe_ho.GetDof();
const int fe_lor_ndof = fe_lor.GetDof();
auto d_B_L = Reshape(B_L.Read(), qPts, fe_lor_ndof, nref);
auto d_B_H = Reshape(B_H.Read(), qPts, fe_ho_ndof, nref);
auto d_D = Reshape(D.Read(), qPts, nref, nel_ho);
mfem::forall(fe_ho_ndof*nref*nel_ho, [=] MFEM_HOST_DEVICE (int idx)
{
real_t dot = 0.0;
for (int qi=0; qi<qPts; ++qi)
const int bh = idx % fe_ho_ndof;
const int iref = (idx / fe_ho_ndof) % nref;
const int iho = idx / fe_ho_ndof / nref;
// (B_lo_dofs x Q) x (Q x B_ho_dofs)
for (int bl = 0; bl < fe_lor_ndof; ++bl)
{
dot += d_B_L(qi, bl, iref) * d_D(qi, iref, iho) * d_B_H(qi, bh, iref);
real_t dot = 0.0;
for (int qi=0; qi<qPts; ++qi)
{
dot += d_B_L(qi, bl, iref) * d_D(qi, iref, iho) * d_B_H(qi, bh, iref);
}
// column major storage
v_M_LH(bl, bh, iref, iho) = dot;
}
// column major storage
v_M_LH(bl, bh, iref, iho) = dot;
}
});
});
} // end of mixed assembly mass matrix
}
L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space
(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_,
CoefficientWithOrder coeff_ho_, CoefficientWithOrder coeff_lor_,
const bool use_ea_, MemoryType d_mt_)
: L2Projection(fes_ho_, fes_lor_, coeff_ho_, coeff_lor_, d_mt_), use_ea(use_ea_)
: L2Projection(fes_ho_, fes_lor_, d_mt_),
use_ea(use_ea_)
{
if (use_ea)
{
@@ -658,11 +559,7 @@ L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space
DenseMatrix Minv_lor(ndof_lor*nref, ndof_lor*nref);
DenseMatrix M_mixed(ndof_lor*nref, ndof_ho);
MassIntegrator mi = [&]()
{
return coeff_lor ? MassIntegrator(*coeff_lor.coeff) : MassIntegrator();
}();
MassIntegrator mi;
DenseMatrix M_lor_el(ndof_lor, ndof_lor);
DenseMatrixInverse Minv_lor_el(&M_lor_el);
DenseMatrix M_lor(ndof_lor*nref, ndof_lor*nref);
@@ -680,10 +577,6 @@ L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space
// Assemble the low-order refined mass matrix and invert locally
int ilor = lor_els[iref];
ElementTransformation *tr_lor = fes_lor.GetElementTransformation(ilor);
const int order = 2*fe_lor.GetOrder() + tr_lor->OrderW() + coeff_lor.order;
mi.SetIntegrationRule(IntRules.Get(geom, order));
mi.AssembleElementMatrix(fe_lor, *tr_lor, M_lor_el);
M_lor.CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
Minv_lor_el.Factor();
@@ -775,22 +668,25 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::EAL2ProjectionL2Space()
// Need to compute M_L
// Note: Using user-inputted M_LH IntegrationRule ir
// (higher order than needed) in order to re-use coeff
MassIntegrator mi = [&]()
{
return coeff_lor ? MassIntegrator(*coeff_lor.coeff) : MassIntegrator();
}();
const int order = 2*fes_lor.GetMaxElementOrder()
+ mesh_lor->GetTypicalElementTransformation()->OrderW()
+ coeff_lor.order;
mi.SetIntegrationRule(
IntRules.Get(mesh_lor->GetTypicalElementGeometry(), order));
MassIntegrator mi;
Vector M_ea_lor;
const int ndof_lor = fes_lor.GetTypicalFE()->GetDof();
const int ndof_ho = fes_ho.GetTypicalFE()->GetDof();
const int nref = ho2lor.RowSize(0);
M_ea_lor.SetSize(ndof_lor*ndof_lor*nel_lor, d_mt);
int ndof_lor;
int ndof_ho;
int nref;
{
int iho = 0;
Array<int> lor_els;
ho2lor.GetRow(iho, lor_els);
nref = ho2lor.RowSize(iho);
const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
const FiniteElement &fe_lor = *fes_lor.GetFE(lor_els[0]);
ndof_ho = fe_ho.GetDof();
ndof_lor = fe_lor.GetDof();
M_ea_lor.SetSize(ndof_lor*ndof_lor*nel_lor, d_mt);
}
const bool add = false;
mi.AssembleEA(fes_lor, M_ea_lor, add);
@@ -1136,9 +1032,8 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::EAProlongateTranspose(
L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_,
CoefficientWithOrder coeff_ho_, CoefficientWithOrder coeff_lor_,
const bool use_ea_, MemoryType d_mt_)
: L2Projection(fes_ho_, fes_lor_, coeff_ho_, coeff_lor_, d_mt_),
: L2Projection(fes_ho_, fes_lor_, d_mt_),
use_ea(use_ea_)
{
@@ -1197,9 +1092,8 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor,
CoefficientWithOrder coeff_ho_, CoefficientWithOrder coeff_lor_,
const bool use_ea_, MemoryType d_mt_)
: L2Projection(pfes_ho, pfes_lor, coeff_ho_, coeff_lor_, d_mt_),
: L2Projection(pfes_ho, pfes_lor, d_mt_),
use_ea(use_ea_), pcg(pfes_ho.GetComm())
{
@@ -1271,12 +1165,12 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetupPCG()
void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
{
Mesh &mesh_ho = *fes_ho.GetMesh();
Mesh &mesh_lor = *fes_lor.GetMesh();
const int nel_ho = mesh_ho.GetNE();
const int nel_lor = mesh_lor.GetNE();
const int ndof_ho = fes_ho.GetNDofs();
const int ndof_lor = fes_lor.GetNDofs();
Mesh* mesh_ho = fes_ho.GetMesh();
Mesh* mesh_lor = fes_lor.GetMesh();
int nel_ho = mesh_ho->GetNE();
int nel_lor = mesh_lor->GetNE();
int ndof_ho = fes_ho.GetNDofs();
int ndof_lor = fes_lor.GetNDofs();
// If the local mesh is empty, skip all computations
if (nel_ho == 0)
@@ -1284,11 +1178,11 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
return;
}
const CoarseFineTransformations& cf_tr = mesh_lor.GetRefinementTransforms();
const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
int nref_max = 0;
Array<Geometry::Type> geoms;
mesh_ho.GetGeometries(mesh_ho.Dimension(), geoms);
mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
for (int ig = 0; ig < geoms.Size(); ++ig)
{
Geometry::Type geom = geoms[ig];
@@ -1311,8 +1205,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
BilinearForm Mho(fes_ho_scalar.get());
Mho.SetAssemblyLevel(AssemblyLevel::PARTIAL);
Mho.AddDomainIntegrator(coeff_ho ? new MassIntegrator(*coeff_ho.coeff)
: new MassIntegrator);
Mho.AddDomainIntegrator(new MassIntegrator);
Mho.Assemble();
// Processor local lumped Mass
@@ -1322,16 +1215,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
BilinearForm Mlor(fes_lor_scalar.get());
Mlor.SetAssemblyLevel(AssemblyLevel::PARTIAL);
{
MassIntegrator *mi = coeff_lor ? new MassIntegrator(*coeff_lor.coeff)
: new MassIntegrator;
const int order = 2*fes_lor.GetMaxElementOrder()
+ mesh_lor.GetTypicalElementTransformation()->OrderW()
+ coeff_lor.order;
mi->SetIntegrationRule(
IntRules.Get(mesh_lor.GetTypicalElementGeometry(), order));
Mlor.AddDomainIntegrator(mi);
}
Mlor.AddDomainIntegrator(new MassIntegrator);
Mlor.Assemble();
Vector ones_lor(Mlor.Width()); ones_lor = 1.0;
@@ -1344,14 +1228,15 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
MixedMassEA(fes_ho, fes_lor, M_LH_ea, d_mt);
// Set ownership
M_LH.reset(new H1SpaceMixedMassOperator(fes_ho_scalar.get(),
fes_lor_scalar.get(),
&ho2lor,
&M_LH_ea));
M_LH_local_op = new H1SpaceMixedMassOperator(fes_ho_scalar.get(),
fes_lor_scalar.get(),
&ho2lor,
&M_LH_ea);
ML_inv_vea.reset(new H1SpaceLumpedMassOperator(fes_ho_scalar.get(),
fes_lor_scalar.get(),
ML_inv_ea));
M_LH.reset(M_LH_local_op);
R.reset(new ProductOperator(ML_inv_vea.get(), M_LH.get(), false,
false));
@@ -1368,18 +1253,18 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
(const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor)
{
Mesh &mesh_ho = *pfes_ho.GetParMesh();
Mesh &mesh_lor = *pfes_lor.GetParMesh();
int nel_ho = mesh_ho.GetNE();
int nel_lor = mesh_lor.GetNE();
Mesh* mesh_ho = pfes_ho.GetParMesh();
Mesh* mesh_lor = pfes_lor.GetParMesh();
int nel_ho = mesh_ho->GetNE();
int nel_lor = mesh_lor->GetNE();
int ndof_ho = pfes_ho.GetNDofs();
int ndof_lor = pfes_lor.GetNDofs();
const CoarseFineTransformations& cf_tr = mesh_lor.GetRefinementTransforms();
const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
int nref_max = 0;
Array<Geometry::Type> geoms;
mesh_ho.GetGeometries(mesh_ho.Dimension(), geoms);
mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
for (int ig = 0; ig < geoms.Size(); ++ig)
{
Geometry::Type geom = geoms[ig];
@@ -1402,8 +1287,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
ParBilinearForm pMho(pfes_ho_scalar.get());
pMho.SetAssemblyLevel(AssemblyLevel::PARTIAL);
pMho.AddDomainIntegrator(coeff_ho ? new MassIntegrator(*coeff_ho.coeff)
: new MassIntegrator);
pMho.AddDomainIntegrator(new MassIntegrator);
pMho.Assemble();
// Processor local lumped Mass
@@ -1413,16 +1297,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
ParBilinearForm pMlor(pfes_lor_scalar.get());
pMlor.SetAssemblyLevel(AssemblyLevel::PARTIAL);
{
MassIntegrator *mi = coeff_lor ? new MassIntegrator(*coeff_lor.coeff)
: new MassIntegrator;
const int order = 2*fes_lor.GetMaxElementOrder()
+ mesh_lor.GetTypicalElementTransformation()->OrderW()
+ coeff_lor.order;
mi->SetIntegrationRule(
IntRules.Get(mesh_lor.GetTypicalElementGeometry(), order));
pMlor.AddDomainIntegrator(mi);
}
pMlor.AddDomainIntegrator(new MassIntegrator);
pMlor.Assemble();
Vector ones_lor(pMlor.Width()); ones_lor = 1.0;
@@ -1695,7 +1570,7 @@ std::unique_ptr<SparseMatrix>>
int ilor = lor_els[iref];
ElementTransformation* el_tr = fes_lor.GetElementTransformation(ilor);
int order = 2 * fe_lor.GetOrder() + el_tr->OrderW() + coeff_lor.order;
int order = 2 * fe_lor.GetOrder() + el_tr->OrderW();
const IntegrationRule* ir = &IntRules.Get(geom, order);
ML_el = 0.0;
for (int i = 0; i < ir->GetNPoints(); ++i)
@@ -1703,13 +1578,7 @@ std::unique_ptr<SparseMatrix>>
const IntegrationPoint& ip_lor = ir->IntPoint(i);
fe_lor.CalcShape(ip_lor, shape_lor);
el_tr->SetIntPoint(&ip_lor);
real_t w = ip_lor.weight;
if (coeff_lor)
{
w *= coeff_lor.coeff->Eval(*el_tr, ip_lor);
}
shape_lor *= el_tr->Weight() * w;
ML_el += shape_lor;
ML_el += (shape_lor *= (el_tr->Weight() * ip_lor.weight));
}
fes_lor.GetElementDofs(ilor, dofs_lor);
ML_inv.AddElementVector(dofs_lor, ML_el);
@@ -2155,8 +2024,8 @@ void L2ProjectionGridTransfer::BuildF()
{
if (!Parallel())
{
F = new L2ProjectionH1Space(
dom_fes, ran_fes, coeff_ho, coeff_lor, use_ea, d_mt);
F = new L2ProjectionH1Space(dom_fes, ran_fes,
use_ea, d_mt);
}
else
{
@@ -2165,15 +2034,15 @@ void L2ProjectionGridTransfer::BuildF()
static_cast<mfem::ParFiniteElementSpace&>(dom_fes);
const mfem::ParFiniteElementSpace& ran_pfes =
static_cast<mfem::ParFiniteElementSpace&>(ran_fes);
F = new L2ProjectionH1Space(
dom_pfes, ran_pfes, coeff_ho, coeff_lor, use_ea, d_mt);
F = new L2ProjectionH1Space(dom_pfes, ran_pfes,
use_ea, d_mt);
#endif
}
}
else
{
F = new L2ProjectionL2Space(
dom_fes, ran_fes, coeff_ho, coeff_lor, use_ea, d_mt);
F = new L2ProjectionL2Space(dom_fes, ran_fes,
use_ea, d_mt);
}
}
+7 -76
View File
@@ -19,8 +19,6 @@
#include "pfespace.hpp"
#endif
#include <cstddef>
namespace mfem
{
@@ -164,18 +162,6 @@ public:
};
struct CoefficientWithOrder
{
Coefficient *coeff;
int order;
CoefficientWithOrder() : coeff(nullptr), order(0) { }
CoefficientWithOrder(std::nullptr_t) : coeff(nullptr), order(0) { }
CoefficientWithOrder(Coefficient &coeff_) : coeff(&coeff_), order(1) { }
CoefficientWithOrder(Coefficient &coeff_, int order_)
: coeff(&coeff_), order(order_) { }
operator bool() const { return coeff != nullptr; }
};
/** @brief Transfer data in L2 and H1 finite element spaces between a coarse
mesh and an embedded refined mesh using L2 projection. */
/** The forward, coarse-to-fine, transfer uses L2 projection. The backward,
@@ -221,8 +207,6 @@ public:
protected:
const FiniteElementSpace& fes_ho;
const FiniteElementSpace& fes_lor;
CoefficientWithOrder coeff_ho;
CoefficientWithOrder coeff_lor;
MemoryType d_mt;
Array<int> offsets;
@@ -230,15 +214,8 @@ public:
L2Projection(const FiniteElementSpace& fes_ho_,
const FiniteElementSpace& fes_lor_,
CoefficientWithOrder coeff_ho_,
CoefficientWithOrder coeff_lor_,
MemoryType d_mt_ = Device::GetHostMemoryType());
L2Projection(const FiniteElementSpace& fes_ho_,
const FiniteElementSpace& fes_lor_,
MemoryType d_mt_ = Device::GetHostMemoryType())
: L2Projection(fes_ho_, fes_lor_, nullptr, nullptr, d_mt_) { }
void BuildHo2Lor(int nel_ho, int nel_lor,
const CoarseFineTransformations& cf_tr);
@@ -248,11 +225,11 @@ public:
IntegrationPointTransformation& ip_tr,
DenseMatrix& M_mixed_el) const;
void ElemMixedEvaluation(Geometry::Type geom, const FiniteElement& fe_ho,
const FiniteElement& fe_lor,
IntegrationPointTransformation& ip_tr,
const IntegrationRule& ir,
DenseMatrix& B_L, DenseMatrix& B_H) const;
void ElemMixedMass(Geometry::Type geom, const FiniteElement& fe_ho,
const FiniteElement& fe_lor,
ElementTransformation* el_tr,
IntegrationPointTransformation& ip_tr,
DenseMatrix& B_L, DenseMatrix& B_H) const;
public:
/* Returns the Mixed Mass M_LH via device element assembly by building the
basis functions and data at the quadrature points. */
@@ -310,17 +287,9 @@ public:
public:
L2ProjectionL2Space(const FiniteElementSpace& fes_ho_,
const FiniteElementSpace& fes_lor_,
CoefficientWithOrder coeff_ho_,
CoefficientWithOrder coeff_lor_,
const bool use_ea_,
MemoryType d_mt_ = Device::GetHostMemoryType());
L2ProjectionL2Space(const FiniteElementSpace& fes_ho_,
const FiniteElementSpace& fes_lor_,
const bool use_ea_,
MemoryType d_mt_ = Device::GetHostMemoryType())
: L2ProjectionL2Space(fes_ho_, fes_lor_, nullptr, nullptr, use_ea_, d_mt_) { }
/*Same as above but assembles and stores R_ea, P_ea */
void EAL2ProjectionL2Space();
@@ -387,30 +356,13 @@ public:
public:
L2ProjectionH1Space(const FiniteElementSpace &fes_ho_,
const FiniteElementSpace &fes_lor_,
CoefficientWithOrder coeff_ho_,
CoefficientWithOrder coeff_lor_,
const bool use_ea_,
MemoryType d_mt_ = Device::GetHostMemoryType());
L2ProjectionH1Space(const FiniteElementSpace& fes_ho_,
const FiniteElementSpace& fes_lor_,
const bool use_ea_,
MemoryType d_mt_ = Device::GetHostMemoryType())
: L2ProjectionH1Space(fes_ho_, fes_lor_, nullptr, nullptr, use_ea_, d_mt_) { }
#ifdef MFEM_USE_MPI
L2ProjectionH1Space(const ParFiniteElementSpace &pfes_ho_,
const ParFiniteElementSpace &pfes_lor_,
CoefficientWithOrder coeff_ho_,
CoefficientWithOrder coeff_lor_,
const bool use_ea_,
MemoryType d_mt_ = Device::GetHostMemoryType());
L2ProjectionH1Space(const ParFiniteElementSpace& fes_ho_,
const ParFiniteElementSpace& fes_lor_,
const bool use_ea_,
MemoryType d_mt_ = Device::GetHostMemoryType())
: L2ProjectionH1Space(fes_ho_, fes_lor_, nullptr, nullptr, use_ea_, d_mt_) { }
#endif
/// Same as above but assembles action of R through 4 parts:
/// ( ) inv( lumped(M_L) ), which is a diagonal matrix (essentially a vector)
@@ -556,38 +508,18 @@ public:
virtual ~L2Prolongation() { }
};
/// Coefficient for the mixed L2 inner product.
CoefficientWithOrder coeff_ho;
/// Coefficient for the low-order L2 inner product.
CoefficientWithOrder coeff_lor;
L2Projection *F; ///< Forward, coarse-to-fine, operator
L2Prolongation *B; ///< Backward, fine-to-coarse, operator
bool force_l2_space;
public:
/// Construct the unweighted L2 projection grid transfer.
L2ProjectionGridTransfer(FiniteElementSpace &coarse_fes_,
FiniteElementSpace &fine_fes_,
bool force_l2_space_ = false,
MemoryType d_mt_ = Device::GetHostMemoryType()) // move to method
: GridTransfer(coarse_fes_, fine_fes_),
coeff_ho(nullptr), coeff_lor(nullptr), F(nullptr), B(nullptr),
force_l2_space(force_l2_space_) { }
/// @brief Construct the weighted L2 projection grid transfer.
///
/// The low-order inner product is weighted by @a coeff_lor, and the mixed
/// inner product is weighted by @a coeff_ho.
L2ProjectionGridTransfer(FiniteElementSpace &coarse_fes_,
FiniteElementSpace &fine_fes_,
CoefficientWithOrder coeff_ho_,
CoefficientWithOrder coeff_lor_,
bool force_l2_space_ = false,
MemoryType d_mt_ = Device::GetHostMemoryType()) // move to method
: GridTransfer(coarse_fes_, fine_fes_),
coeff_ho(coeff_ho_), coeff_lor(coeff_lor_), F(nullptr), B(nullptr),
force_l2_space(force_l2_space_) { }
F(NULL), B(NULL), force_l2_space(force_l2_space_)
{ }
virtual ~L2ProjectionGridTransfer();
const Operator &ForwardOperator() override;
@@ -595,7 +527,6 @@ public:
const Operator &BackwardOperator() override;
bool SupportsBackwardsOperator() const override;
private:
void BuildF();
};
-158
View File
@@ -1108,126 +1108,6 @@ void GroupCommunicator::ReduceEnd(T *ldata, int layout,
num_requests = 0;
}
template <class T>
void GroupCommunicator::ReduceMarked(T *ldata, const Array<int> &marker,
int layout,
void (*Op)(OpData<T>)) const
{
if (comm_lock == 0) { return; }
// The above also handles the case (group_buf_size == 0).
MFEM_VERIFY(comm_lock == 2, "object is NOT locked for Reduce");
switch (mode)
{
case byGroup: // ***** Communication by groups *****
{
OpData<T> opd;
opd.ldata = ldata;
Array<int> group_num_req(group_ldof.Size());
for (int gr = 1; gr < group_ldof.Size(); gr++)
{
group_num_req[gr] =
gtopo.IAmMaster(gr) ? gtopo.GetGroupSize(gr)-1 : 0;
}
int idx;
while (MPI_Waitany(num_requests, requests, &idx, MPI_STATUS_IGNORE),
idx != MPI_UNDEFINED)
{
int gr = request_marker[idx];
if (gr == -1) { continue; } // skip send requests
// Delay the processing of a group until all receive requests, for
// that group, are done:
if ((--group_num_req[gr]) != 0) { continue; }
opd.nldofs = group_ldof.RowSize(gr);
// groups without dofs are skipped, so here nldofs > 0.
opd.buf = (T *)group_buf.GetData() + buf_offsets[gr];
opd.ldofs = (layout == 0) ?
group_ldof.GetRow(gr) : group_ltdof.GetRow(gr);
opd.nb = gtopo.GetGroupSize(gr)-1;
// Apply operation only to marked DOFs. The receive buffer is
// neighbor-major with stride opd.nldofs, i.e. the contributions to
// DOF i are buf[j*opd.nldofs + i] for j = 0 ... opd.nb-1. Setting
// nldofs = 1 for a single DOF changes that stride to 1, so the
// strided values must first be gathered into a contiguous buffer.
Array<T> single_buf(opd.nb);
for (int i = 0; i < opd.nldofs; i++)
{
if (marker[opd.ldofs[i]])
{
for (int j = 0; j < opd.nb; j++)
{
single_buf[j] = opd.buf[j*opd.nldofs + i];
}
// Create a temporary OpData with just this one DOF
OpData<T> single_opd;
single_opd.ldata = ldata;
single_opd.buf = single_buf.GetData();
single_opd.ldofs = opd.ldofs + i;
single_opd.nldofs = 1;
single_opd.nb = opd.nb;
// Apply the operation
Op(single_opd);
}
}
}
break;
}
case byNeighbor: // ***** Communication by neighbors *****
{
MPI_Waitall(num_requests, requests, MPI_STATUSES_IGNORE);
for (int nbr = 1; nbr < nbr_send_groups.Size(); nbr++)
{
// In Reduce operation: send_groups <--> recv_groups
const int num_recv_groups = nbr_send_groups.RowSize(nbr);
if (num_recv_groups > 0)
{
const int *grp_list = nbr_send_groups.GetRow(nbr);
const T *buf = (T*)group_buf.GetData() + buf_offsets[nbr];
for (int i = 0; i < num_recv_groups; i++)
{
// Custom version of ReduceGroupFromBuffer that checks marker
int gr = grp_list[i];
const int *ldofs = (layout == 0) ?
group_ldof.GetRow(gr) : group_ltdof.GetRow(gr);
const int nldofs = group_ldof.RowSize(gr);
for (int j = 0; j < nldofs; j++)
{
if (marker[ldofs[j]])
{
// Create a temporary OpData with just this one DOF
OpData<T> opd;
opd.ldata = ldata;
opd.buf = const_cast<T*>(buf) + j;
opd.ldofs = ldofs + j;
opd.nldofs = 1;
opd.nb = 1;
// Apply the operation
Op(opd);
}
}
buf += nldofs;
}
}
}
break;
}
}
comm_lock = 0; // 0 - no lock
num_requests = 0;
}
template <class T>
void GroupCommunicator::Sum(OpData<T> opd)
{
@@ -1291,8 +1171,6 @@ void GroupCommunicator::Max(OpData<T> opd)
template <class T>
void GroupCommunicator::BitOR(OpData<T> opd)
{
static_assert(std::is_integral<T>::value,
"BitOR reduction requires an integral type.");
for (int i = 0; i < opd.nldofs; i++)
{
T data = opd.ldata[opd.ldofs[i]];
@@ -1304,33 +1182,6 @@ void GroupCommunicator::BitOR(OpData<T> opd)
}
}
template <class T>
void GroupCommunicator::MaxAbs(OpData<T> opd)
{
for (int i = 0; i < opd.nldofs; i++)
{
T data = opd.ldata[opd.ldofs[i]];
T abs_data = std::abs(data);
for (int j = 0; j < opd.nb; j++)
{
T b = opd.buf[j*opd.nldofs+i];
T abs_b = std::abs(b);
// On an equal-magnitude tie keep the more positive value, so
// opposite-sign ties resolve deterministically to the positive one.
if (abs_data < abs_b || (abs_data == abs_b && data < b))
{
data = b;
abs_data = abs_b;
}
}
opd.ldata[opd.ldofs[i]] = data;
}
}
void GroupCommunicator::PrintInfo(std::ostream &os) const
{
char c = '\0';
@@ -1467,24 +1318,18 @@ template void GroupCommunicator::BcastEnd<int>(int *, int) const;
template void GroupCommunicator::ReduceBegin<int>(const int *) const;
template void GroupCommunicator::ReduceEnd<int>(
int *, int, void (*)(OpData<int>)) const;
template void GroupCommunicator::ReduceMarked<int>(
int*, const Array<int>&, int, void (*)(OpData<int>)) const;
template void GroupCommunicator::BcastBegin<double>(double *, int) const;
template void GroupCommunicator::BcastEnd<double>(double *, int) const;
template void GroupCommunicator::ReduceBegin<double>(const double *) const;
template void GroupCommunicator::ReduceEnd<double>(
double *, int, void (*)(OpData<double>)) const;
template void GroupCommunicator::ReduceMarked<double>(
double*, const Array<int>&, int, void (*)(OpData<double>)) const;
template void GroupCommunicator::BcastBegin<float>(float *, int) const;
template void GroupCommunicator::BcastEnd<float>(float *, int) const;
template void GroupCommunicator::ReduceBegin<float>(const float *) const;
template void GroupCommunicator::ReduceEnd<float>(
float *, int, void (*)(OpData<float>)) const;
template void GroupCommunicator::ReduceMarked<float>(
float*, const Array<int>&, int, void (*)(OpData<float>)) const;
// @endcond
@@ -1493,17 +1338,14 @@ template void GroupCommunicator::Sum<int>(OpData<int>);
template void GroupCommunicator::Min<int>(OpData<int>);
template void GroupCommunicator::Max<int>(OpData<int>);
template void GroupCommunicator::BitOR<int>(OpData<int>);
template void GroupCommunicator::MaxAbs<int>(OpData<int>);
template void GroupCommunicator::Sum<double>(OpData<double>);
template void GroupCommunicator::Min<double>(OpData<double>);
template void GroupCommunicator::Max<double>(OpData<double>);
template void GroupCommunicator::MaxAbs<double>(OpData<double>);
template void GroupCommunicator::Sum<float>(OpData<float>);
template void GroupCommunicator::Min<float>(OpData<float>);
template void GroupCommunicator::Max<float>(OpData<float>);
template void GroupCommunicator::MaxAbs<float>(OpData<float>);
#ifdef __bgq__
+3 -28
View File
@@ -22,7 +22,6 @@
#include "globals.hpp"
#include <mpi.h>
#include <cstdint>
#include <type_traits>
// can't directly use MPI_CXX_BOOL because Microsoft's MPI implementation
// doesn't include MPI_CXX_BOOL. Fallback to MPI_C_BOOL if unavailable.
@@ -409,38 +408,14 @@ public:
template <class T> void Reduce(Array<T> &ldata, void (*Op)(OpData<T>)) const
{ Reduce<T>((T *)ldata, Op); }
/// Reduce operation Sum, instantiated for int, double and float
/// Reduce operation Sum, instantiated for int and double
template <class T> static void Sum(OpData<T>);
/// Reduce operation Min, instantiated for int, double and float
/// Reduce operation Min, instantiated for int and double
template <class T> static void Min(OpData<T>);
/// Reduce operation Max, instantiated for int, double and float
/// Reduce operation Max, instantiated for int and double
template <class T> static void Max(OpData<T>);
/// Reduce operation bitwise OR, instantiated for int only
template <class T> static void BitOR(OpData<T>);
/// Reduce operation selecting the signed value with the largest absolute
/// value, instantiated for int, double and float. The result keeps its sign;
/// it is not the non-negative absolute value. Equal-magnitude ties are
/// broken deterministically toward the more positive value, so opposite-sign
/// ties resolve to the positive one regardless of accumulation order.
template <class T> static void MaxAbs(OpData<T>);
/** @brief Finalize reduction operation started with ReduceBegin(), but only apply
the reduction to DOFs marked in the marker array.
@note The reduction is carried out in the signed type @a T, so the result
is signed even for bitwise operations.
*/
template <class T>
void ReduceMarked(T *ldata, const Array<int> &marker, int layout,
void (*Op)(OpData<T>)) const;
/** @brief Reduce within each group where the master is the root, but only for marked DOFs. */
template <class T>
void Reduce(T *ldata, const Array<int> &marker, void (*Op)(OpData<T>)) const
{
ReduceBegin(ldata);
ReduceMarked(ldata, marker, 0, Op);
}
/// Print information about the GroupCommunicator from all MPI ranks.
void PrintInfo(std::ostream &out = mfem::out) const;
+4
View File
@@ -36,6 +36,8 @@ list(APPEND SRCS
sparsemat.cpp
sparsesmoothers.cpp
symmat.cpp
tripack.cpp
tripack_magma.cpp
vector.cpp
)
@@ -73,6 +75,8 @@ list(APPEND HDRS
sparsemat.hpp
sparsesmoothers.hpp
symmat.hpp
tripack.hpp
tripack_magma.hpp
tlayout.hpp
tmatrix.hpp
ttensor.hpp
+4
View File
@@ -27,6 +27,10 @@
#include "sparsesmoothers.hpp"
#include "densemat.hpp"
#include "symmat.hpp"
#include "tripack.hpp"
#ifdef MFEM_USE_MAGMA
#include "tripack_magma.hpp"
#endif
#include "ode.hpp"
#include "solvers.hpp"
#include "handle.hpp"
+471
View File
@@ -0,0 +1,471 @@
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#include "tripack.hpp"
#include "../general/forall.hpp"
#include <cmath>
#include <limits>
namespace mfem
{
namespace
{
MFEM_HOST_DEVICE inline int SymmetricIndex(const int i,
const int j,
const int n)
{
return TriPackLowerMatrix::Index(i, j, n);
}
MFEM_HOST_DEVICE inline bool TriPackIsFinite(const real_t val)
{
#ifdef isfinite
return isfinite(val);
#else
return std::isfinite(val);
#endif
}
void ComputeCholeskyFactorsLowerDevice(
const TriPackLowerMatrix &packed_lower,
Vector &factor)
{
const int n = packed_lower.GetNumRows();
const int batch_size = packed_lower.GetNumMatrices();
const int packed_size = packed_lower.GetPackedSize();
const real_t nan = std::numeric_limits<real_t>::quiet_NaN();
const real_t *A = packed_lower.Data().Read();
factor.SetSize(batch_size*packed_size);
factor.UseDevice(true);
real_t *L = factor.Write();
mfem::forall(batch_size, [=] MFEM_HOST_DEVICE (int e)
{
const int eoff = e*packed_size;
const real_t eps = std::numeric_limits<real_t>::epsilon();
bool bad = false;
// Copy packed-lower input into factor storage (no scaling).
for (int j = 0; j < n; ++j)
{
for (int i = j; i < n; ++i)
{
const int t = TriPackLowerMatrix::LowerIndex(i, j, n);
const real_t Aij = A[eoff + t];
if (!TriPackIsFinite(Aij)) { bad = true; }
L[eoff + t] = Aij;
}
}
if (bad)
{
for (int t = 0; t < packed_size; ++t) { L[eoff + t] = nan; }
return;
}
for (int k = 0; k < n; ++k)
{
const int kk = eoff + TriPackLowerMatrix::LowerIndex(k, k, n);
const real_t Lkk0 = L[kk];
real_t Lkk = Lkk0;
for (int s = 0; s < k; ++s)
{
const real_t Lks = L[eoff + TriPackLowerMatrix::LowerIndex(k, s, n)];
Lkk -= Lks*Lks;
}
const real_t tol = 64.0*eps*fabs(Lkk0);
if (!TriPackIsFinite(Lkk) || Lkk < -tol)
{
bad = true;
break;
}
if (Lkk < 0.0) { Lkk = 0.0; }
L[kk] = sqrt(Lkk);
const real_t Ldiag = L[kk];
for (int i = k + 1; i < n; ++i)
{
const int ik = eoff + TriPackLowerMatrix::LowerIndex(i, k, n);
real_t Aik = L[ik];
for (int s = 0; s < k; ++s)
{
Aik -= L[eoff + TriPackLowerMatrix::LowerIndex(i, s, n)] *
L[eoff + TriPackLowerMatrix::LowerIndex(k, s, n)];
}
L[ik] = Aik/Ldiag;
}
}
if (bad)
{
for (int t = 0; t < packed_size; ++t) { L[eoff + t] = nan; }
}
});
}
void ComputeCholeskyFactorsLower(
const TriPackLowerMatrix &packed_lower,
Vector &factor)
{
const int n = packed_lower.GetNumRows();
const int batch_size = packed_lower.GetNumMatrices();
const int packed_size = packed_lower.GetPackedSize();
factor.SetSize(batch_size*packed_size);
factor.UseDevice(true);
const real_t *A = packed_lower.Data().HostRead();
real_t *L = factor.HostWrite();
for (int e = 0; e < batch_size; ++e)
{
const int eoff = e*packed_size;
for (int j = 0; j < n; ++j)
{
for (int i = j; i < n; ++i)
{
L[eoff + TriPackLowerMatrix::LowerIndex(i, j, n)] =
A[eoff + TriPackLowerMatrix::LowerIndex(i, j, n)];
}
}
for (int k = 0; k < n; ++k)
{
const int kk = eoff + TriPackLowerMatrix::LowerIndex(k, k, n);
real_t Akk = L[kk];
for (int s = 0; s < k; ++s)
{
const real_t Lks = L[eoff + TriPackLowerMatrix::LowerIndex(k, s, n)];
Akk -= Lks*Lks;
}
MFEM_VERIFY(Akk > 0.0, "Matrix is not SPD.");
L[kk] = std::sqrt(Akk);
const real_t Lkk = L[kk];
for (int i = k + 1; i < n; ++i)
{
real_t Aik = L[eoff + TriPackLowerMatrix::LowerIndex(i, k, n)];
for (int s = 0; s < k; ++s)
{
Aik -= L[eoff + TriPackLowerMatrix::LowerIndex(i, s, n)] *
L[eoff + TriPackLowerMatrix::LowerIndex(k, s, n)];
}
L[eoff + TriPackLowerMatrix::LowerIndex(i, k, n)] = Aik / Lkk;
}
}
}
}
}
namespace tripack
{
bool CompareWithFull(const TriPackLowerMatrix &packed, const Vector &full,
real_t tol)
{
const int n = packed.GetNumRows();
const int batch_size = packed.GetNumMatrices();
const int packed_size = packed.GetPackedSize();
MFEM_VERIFY(full.Size() == batch_size*n*n,
"Full matrix data has the wrong size.");
const real_t *packed_data = packed.Data().HostRead();
const real_t *full_data = full.HostRead();
if (tol == 0.0)
{
tol = 256.0*std::numeric_limits<real_t>::epsilon();
}
for (int e = 0; e < batch_size; ++e)
{
const int packed_offset = e*packed_size;
const int full_offset = e*n*n;
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < n; ++j)
{
const real_t packed_val =
packed_data[packed_offset + SymmetricIndex(i, j, n)];
const real_t full_val = full_data[full_offset + i + n*j];
if (std::fabs(full_val - packed_val) > tol)
{
return false;
}
}
}
for (int i = 0; i < n; ++i)
{
for (int j = i + 1; j < n; ++j)
{
const real_t a = full_data[full_offset + i + n*j];
const real_t b = full_data[full_offset + j + n*i];
if (std::fabs(a - b) > tol)
{
return false;
}
}
}
}
return true;
}
void Mult(const TriPackLowerMatrix &packed, const Vector &x, Vector &y)
{
const int n = packed.GetNumRows();
const int batch_size = packed.GetNumMatrices();
MFEM_VERIFY(x.Size() == batch_size*n, "Input vector has the wrong size.");
y.SetSize(batch_size*n);
y.UseDevice(true);
const real_t *A = packed.Data().Read();
const real_t *X = x.Read();
real_t *Y = y.Write();
mfem::forall(batch_size*n, [=] MFEM_HOST_DEVICE (int idx)
{
const int i = idx % n;
const int e = idx / n;
const real_t *Ae = A + e*TriPackLowerMatrix::PackedSize(n);
const real_t *Xe = X + e*n;
real_t sum = 0.0;
for (int j = 0; j < n; ++j)
{
sum += Ae[SymmetricIndex(i, j, n)] * Xe[j];
}
Y[idx] = sum;
});
}
void Lump(const TriPackLowerMatrix &packed, Vector &lump)
{
const int n = packed.GetNumRows();
const int batch_size = packed.GetNumMatrices();
lump.SetSize(batch_size*n);
lump.UseDevice(true);
const real_t *A = packed.Data().Read();
real_t *L = lump.Write();
mfem::forall(batch_size*n, [=] MFEM_HOST_DEVICE (int idx)
{
const int i = idx % n;
const int e = idx / n;
const real_t *Ae = A + e*TriPackLowerMatrix::PackedSize(n);
real_t sum = 0.0;
for (int j = 0; j < n; ++j)
{
sum += Ae[SymmetricIndex(i, j, n)];
}
L[idx] = sum;
});
}
void ComputeCholeskyLower(const TriPackLowerMatrix &packed_lower,
TriPackLowerMatrix &lower_factor)
{
const int n = packed_lower.GetNumRows();
const int batch_size = packed_lower.GetNumMatrices();
MFEM_VERIFY(&packed_lower != &lower_factor,
"Input and output TriPackLowerMatrix objects must be distinct.");
if (batch_size == 0)
{
lower_factor.SetSize(n, batch_size);
return;
}
Vector factored;
ComputeCholeskyFactorsLower(packed_lower, factored);
lower_factor.SetSize(n, batch_size);
lower_factor.UseDevice(true);
lower_factor.Data() = factored;
}
void ComputeCholeskyLowerInverse(const TriPackLowerMatrix &packed_lower,
TriPackLowerMatrix &lower_inverse)
{
const int n = packed_lower.GetNumRows();
const int batch_size = packed_lower.GetNumMatrices();
const int packed_size = packed_lower.GetPackedSize();
const real_t nan = std::numeric_limits<real_t>::quiet_NaN();
MFEM_VERIFY(&packed_lower != &lower_inverse,
"Input and output TriPackLowerMatrix objects must be distinct.");
if (batch_size == 0)
{
lower_inverse.SetSize(n, batch_size);
return;
}
lower_inverse.SetSize(n, batch_size);
lower_inverse.UseDevice(true);
Vector factored;
Vector work(batch_size*packed_size);
work.UseDevice(true);
ComputeCholeskyFactorsLowerDevice(packed_lower, factored);
const real_t *L = factored.Read();
real_t *X = work.Write();
real_t *Linv = lower_inverse.Data().Write();
mfem::forall(batch_size, [=] MFEM_HOST_DEVICE (int e)
{
const int eoff = e*packed_size;
bool bad = false;
for (int t = 0; t < packed_size; ++t)
{
if (!TriPackIsFinite(L[eoff + t]))
{
bad = true;
break;
}
}
if (bad)
{
for (int t = 0; t < packed_size; ++t) { Linv[eoff + t] = nan; }
return;
}
for (int t = 0; t < packed_size; ++t) { X[eoff + t] = 0.0; }
// Compute X = L^{-1} (packed lower).
for (int j = 0; j < n; ++j)
{
const int jj = eoff + TriPackLowerMatrix::LowerIndex(j, j, n);
X[jj] = 1.0/L[jj];
const real_t invLjj = X[jj];
for (int i = j + 1; i < n; ++i)
{
real_t sum = 0.0;
for (int k = j; k < i; ++k)
{
sum += L[eoff + TriPackLowerMatrix::LowerIndex(i, k, n)] *
X[eoff + TriPackLowerMatrix::LowerIndex(k, j, n)];
}
X[eoff + TriPackLowerMatrix::LowerIndex(i, j, n)] = -invLjj*sum;
}
}
for (int t = 0; t < packed_size; ++t)
{
Linv[eoff + t] = X[eoff + t];
}
});
}
void SolveLower(const TriPackLowerMatrix &lower_factor,
const Vector &rhs,
Vector &sol)
{
const int n = lower_factor.GetNumRows();
const int batch_size = lower_factor.GetNumMatrices();
const int packed_size = lower_factor.GetPackedSize();
MFEM_VERIFY(rhs.Size() == batch_size*n, "Right-hand side has the wrong size.");
Vector out(batch_size*n);
out.UseDevice(true);
const real_t *L = lower_factor.Data().Read();
const real_t *B = rhs.Read();
real_t *X = out.Write();
mfem::forall(batch_size, [=] MFEM_HOST_DEVICE (int e)
{
const real_t *Le = L + e*packed_size;
const real_t *Be = B + e*n;
real_t *Xe = X + e*n;
for (int i = 0; i < n; ++i)
{
real_t sum = Be[i];
for (int j = 0; j < i; ++j)
{
sum -= Le[TriPackLowerMatrix::LowerIndex(i, j, n)] * Xe[j];
}
Xe[i] = sum / Le[TriPackLowerMatrix::LowerIndex(i, i, n)];
}
});
sol.SetSize(batch_size*n);
sol = out;
}
void SolveLowerTranspose(const TriPackLowerMatrix &lower_factor,
const Vector &rhs,
Vector &sol)
{
const int n = lower_factor.GetNumRows();
const int batch_size = lower_factor.GetNumMatrices();
const int packed_size = lower_factor.GetPackedSize();
MFEM_VERIFY(rhs.Size() == batch_size*n, "Right-hand side has the wrong size.");
Vector out(batch_size*n);
out.UseDevice(true);
const real_t *L = lower_factor.Data().Read();
const real_t *B = rhs.Read();
real_t *X = out.Write();
mfem::forall(batch_size, [=] MFEM_HOST_DEVICE (int e)
{
const real_t *Le = L + e*packed_size;
const real_t *Be = B + e*n;
real_t *Xe = X + e*n;
for (int i = n - 1; i >= 0; --i)
{
real_t sum = Be[i];
for (int j = i + 1; j < n; ++j)
{
sum -= Le[TriPackLowerMatrix::LowerIndex(j, i, n)] * Xe[j];
}
Xe[i] = sum / Le[TriPackLowerMatrix::LowerIndex(i, i, n)];
}
});
sol.SetSize(batch_size*n);
sol = out;
}
void SolveCholeskyLower(const TriPackLowerMatrix &lower_factor,
const Vector &rhs,
Vector &sol)
{
Vector tmp;
SolveLower(lower_factor, rhs, tmp);
SolveLowerTranspose(lower_factor, tmp, sol);
}
} // namespace tripack
} // namespace mfem
+120
View File
@@ -0,0 +1,120 @@
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#ifndef MFEM_TRIPACK
#define MFEM_TRIPACK
#include "../config/config.hpp"
#include "vector.hpp"
namespace mfem
{
/// Packed storage for a batch of symmetric matrices of fixed size.
/// Storage is packed lower-triangular in LAPACK/MAGMA column-major convention.
class TriPackLowerMatrix
{
private:
Vector data;
int nrows = 0;
int nmats = 0;
public:
TriPackLowerMatrix() = default;
TriPackLowerMatrix(int n, int batch_size)
{
SetSize(n, batch_size);
}
MFEM_HOST_DEVICE static int PackedSize(const int n)
{
return n*(n + 1)/2;
}
/// Packed index for (i,j) in the lower triangle (requires i >= j).
MFEM_HOST_DEVICE static int LowerIndex(const int i, const int j, const int n)
{
return j*(2*n + 1 - j)/2 + (i - j);
}
/// Packed index for (i,j) in symmetric storage (maps to lower triangle).
MFEM_HOST_DEVICE static int Index(const int i, const int j, const int n)
{
return (i >= j) ? LowerIndex(i, j, n) : LowerIndex(j, i, n);
}
void SetSize(const int n, const int batch_size)
{
nrows = n;
nmats = batch_size;
data.SetSize(batch_size*PackedSize(n));
}
int GetNumRows() const { return nrows; }
int GetNumMatrices() const { return nmats; }
int GetPackedSize() const { return PackedSize(nrows); }
int Size() const { return data.Size(); }
void UseDevice(bool use_dev) { data.UseDevice(use_dev); }
TriPackLowerMatrix &operator=(real_t value)
{
data = value;
return *this;
}
Vector &Data() { return data; }
const Vector &Data() const { return data; }
};
namespace tripack
{
bool CompareWithFull(const TriPackLowerMatrix &packed, const Vector &full,
real_t tol = 0.0);
void Mult(const TriPackLowerMatrix &packed, const Vector &x, Vector &y);
void Lump(const TriPackLowerMatrix &packed, Vector &lump);
void ComputeCholeskyLower(const TriPackLowerMatrix &packed_lower,
TriPackLowerMatrix &lower_factor);
void SolveLower(const TriPackLowerMatrix &lower_factor,
const Vector &rhs,
Vector &sol);
void SolveLowerTranspose(const TriPackLowerMatrix &lower_factor,
const Vector &rhs,
Vector &sol);
void SolveCholeskyLower(const TriPackLowerMatrix &lower_factor,
const Vector &rhs,
Vector &sol);
/// Compute the inverse of the Cholesky lower factor for a batch of SPD matrices.
///
/// Given packed lower-triangular matrices A (SPD), this routine computes the
/// packed lower-triangular matrices L^{-1}, where A = L L^T.
///
/// This is intended for fast inverse applications using (L^{-1})^T (L^{-1}).
void ComputeCholeskyLowerInverse(const TriPackLowerMatrix &packed_lower,
TriPackLowerMatrix &lower_inverse);
} // namespace tripack
} // namespace mfem
#endif
+267
View File
@@ -0,0 +1,267 @@
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#include "tripack_magma.hpp"
#ifdef MFEM_USE_MAGMA
#include "../general/forall.hpp"
namespace mfem
{
namespace
{
#ifdef MFEM_USE_SINGLE
#define MFEM_TRIPACK_MAGMA_PREFIX(stub) magma_s##stub
#define MFEM_TRIPACK_MAGMA_SET_POINTER magma_sset_pointer
#elif defined(MFEM_USE_DOUBLE)
#define MFEM_TRIPACK_MAGMA_PREFIX(stub) magma_d##stub
#define MFEM_TRIPACK_MAGMA_SET_POINTER magma_dset_pointer
#else
#error "Unsupported MFEM precision for MAGMA packed routines."
#endif
real_t **SetPackedPointerArray(Array<real_t *> &ptrs,
real_t *data,
const int stride,
const int batch_size,
const magma_queue_t queue)
{
if (ptrs.Size() != batch_size)
{
if (ptrs.Size() != 0) { magma_queue_sync(queue); }
ptrs.SetSize(batch_size, Device::GetDeviceMemoryType());
}
real_t **d_ptrs = ptrs.Write();
MFEM_TRIPACK_MAGMA_SET_POINTER(d_ptrs, data, 1, 0, 0, stride,
batch_size, queue);
return d_ptrs;
}
} // namespace
MagmaPackedLowerCholesky::MagmaPackedLowerCholesky()
{
queue = Magma::Queue();
}
void MagmaPackedLowerCholesky::Factor(
const TriPackLowerMatrix &A,
TriPackLowerMatrix &L)
{
MFEM_VERIFY(queue != nullptr, "MAGMA queue is not set.");
n = A.GetNumRows();
batch_size = A.GetNumMatrices();
packed_size = A.GetPackedSize();
L.SetSize(n, batch_size);
L.UseDevice(true);
if (batch_size == 0) { return; }
L.Data() = A.Data();
real_t *factor_data = L.Data().ReadWrite();
real_t **d_factor_ptrs =
SetPackedPointerArray(factor_ptrs, factor_data, packed_size,
batch_size, queue);
info.SetSize(batch_size, Device::GetDeviceMemoryType());
magma_int_t *d_info = info.Write();
magma_memset(d_info, 0, batch_size*sizeof(magma_int_t));
const magma_int_t status =
(n <= 8) ?
MFEM_TRIPACK_MAGMA_PREFIX(pptrf_batched_small)(
MagmaLower, n, d_factor_ptrs, d_info, batch_size, queue) :
MFEM_TRIPACK_MAGMA_PREFIX(pptf2_batched_small)(
MagmaLower, n, d_factor_ptrs, d_info, batch_size, queue);
MFEM_VERIFY(status == MAGMA_SUCCESS,
"MAGMA packed Cholesky factorization failed.");
magma_queue_sync(queue);
const magma_int_t *h_info = info.HostRead();
for (int e = 0; e < batch_size; ++e)
{
MFEM_VERIFY(h_info[e] == 0,
"MAGMA packed Cholesky factorization failed on matrix "
<< e << '.');
}
}
void MagmaPackedLowerCholesky::SolveInPlace(
const TriPackLowerMatrix &L,
Vector &rhs_sol) const
{
MFEM_VERIFY(queue != nullptr, "MAGMA queue is not set.");
MFEM_VERIFY(L.GetNumRows() > 0 || L.GetNumMatrices() == 0,
"Invalid factor dimensions.");
const int solve_n = L.GetNumRows();
const int solve_batch = L.GetNumMatrices();
const int solve_packed = L.GetPackedSize();
MFEM_VERIFY(rhs_sol.Size() == solve_batch*solve_n,
"Right-hand side has the wrong size.");
if (solve_batch == 0) { return; }
real_t *factor_data = const_cast<real_t *>(L.Data().Read());
real_t **d_factor_ptrs =
SetPackedPointerArray(factor_ptrs, factor_data, solve_packed,
solve_batch, queue);
real_t *rhs_data = rhs_sol.ReadWrite();
real_t **d_rhs_ptrs =
SetPackedPointerArray(rhs_ptrs, rhs_data, solve_n, solve_batch, queue);
const magma_int_t status =
MFEM_TRIPACK_MAGMA_PREFIX(pptrs_batched)(
MagmaLower, solve_n, 1, d_factor_ptrs, d_rhs_ptrs, solve_n,
solve_batch, queue);
MFEM_VERIFY(status == MAGMA_SUCCESS,
"MAGMA packed Cholesky solve failed.");
}
MagmaPackedLowerInverse::MagmaPackedLowerInverse()
{
queue = Magma::Queue();
}
void MagmaPackedLowerInverse::Compute(
const TriPackLowerMatrix &A,
TriPackLowerMatrix &A_inv)
{
MFEM_VERIFY(queue != nullptr, "MAGMA queue is not set.");
n = A.GetNumRows();
batch_size = A.GetNumMatrices();
packed_size = A.GetPackedSize();
MFEM_VERIFY(n <= 64, "MAGMA packed inverse supports n <= 64.");
A_inv.SetSize(n, batch_size);
A_inv.UseDevice(true);
if (batch_size == 0) { return; }
A_inv.Data() = A.Data();
real_t *inv_data = A_inv.Data().ReadWrite();
real_t **d_inv_ptrs =
SetPackedPointerArray(inv_ptrs, inv_data, packed_size, batch_size, queue);
info.SetSize(batch_size, Device::GetDeviceMemoryType());
magma_int_t *d_info = info.Write();
magma_memset(d_info, 0, batch_size*sizeof(magma_int_t));
// MAGMA currently expects a valid pointer for device_lwork even when the
// required workspace is 0 bytes.
int64_t device_lwork[1] = {0};
const magma_int_t status =
MFEM_TRIPACK_MAGMA_PREFIX(ppinv_batched)(
MagmaLower, n, d_inv_ptrs,
/*device_work*/ nullptr, device_lwork,
d_info, batch_size, queue);
MFEM_VERIFY(status == MAGMA_SUCCESS, "MAGMA packed inverse failed.");
magma_queue_sync(queue);
const magma_int_t *h_info = info.HostRead();
for (int e = 0; e < batch_size; ++e)
{
MFEM_VERIFY(h_info[e] == 0,
"MAGMA packed inverse failed on matrix " << e << '.');
}
}
void MagmaPackedLowerInverse::ApplyInPlace(
const TriPackLowerMatrix &A_inv,
Vector &rhs_sol) const
{
MFEM_VERIFY(queue != nullptr, "MAGMA queue is not set.");
MFEM_VERIFY(A_inv.GetNumRows() > 0 || A_inv.GetNumMatrices() == 0,
"Invalid inverse dimensions.");
const int apply_n = A_inv.GetNumRows();
const int apply_batch = A_inv.GetNumMatrices();
const int apply_packed = A_inv.GetPackedSize();
MFEM_VERIFY(rhs_sol.Size() == apply_batch*apply_n,
"Right-hand side has the wrong size.");
if (apply_batch == 0) { return; }
// Prefer MAGMA's tuned packed-symmetric matvec when available (n <= 32).
// Fall back to an MFEM device kernel for larger n.
if (apply_n <= 32)
{
real_t *inv_data = const_cast<real_t *>(A_inv.Data().Read());
real_t **d_inv_ptrs =
SetPackedPointerArray(inv_ptrs, inv_data, apply_packed, apply_batch,
queue);
real_t *rhs_data = rhs_sol.ReadWrite();
real_t **d_rhs_ptrs =
SetPackedPointerArray(rhs_ptrs, rhs_data, apply_n, apply_batch, queue);
// Note: MAGMA's symv_packed_inplace_batched_small returns void; it will
// report argument errors via magma_xerbla.
MFEM_TRIPACK_MAGMA_PREFIX(symv_packed_inplace_batched_small)(
MagmaLower, apply_n, d_inv_ptrs, d_rhs_ptrs, apply_n, apply_batch,
queue);
return;
}
work.SetSize(apply_batch*apply_n);
work.UseDevice(true);
const real_t *AP = A_inv.Data().Read();
const real_t *X = rhs_sol.Read();
real_t *Y = work.Write();
mfem::forall(apply_batch*apply_n, [=] MFEM_HOST_DEVICE (int idx)
{
const int i = idx % apply_n;
const int e = idx / apply_n;
const real_t *APe = AP + e*apply_packed;
const real_t *Xe = X + e*apply_n;
real_t sum = 0.0;
for (int j = 0; j < apply_n; ++j)
{
const real_t aij =
(i >= j) ?
APe[TriPackLowerMatrix::LowerIndex(i, j, apply_n)] :
APe[TriPackLowerMatrix::LowerIndex(j, i, apply_n)];
sum += aij * Xe[j];
}
Y[idx] = sum;
});
const real_t *Y_in = work.Read();
real_t *X_out = rhs_sol.Write();
mfem::forall(apply_batch*apply_n, [=] MFEM_HOST_DEVICE (int idx)
{
X_out[idx] = Y_in[idx];
});
}
#undef MFEM_TRIPACK_MAGMA_SET_POINTER
#undef MFEM_TRIPACK_MAGMA_PREFIX
} // namespace mfem
#endif // MFEM_USE_MAGMA
+144
View File
@@ -0,0 +1,144 @@
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#ifndef MFEM_TRIPACK_MAGMA
#define MFEM_TRIPACK_MAGMA
#include "../config/config.hpp"
#include "tripack.hpp"
#ifdef MFEM_USE_MAGMA
#include "batched/magma.hpp"
namespace mfem
{
/// Workspace + operations for MAGMA packed-lower batched Cholesky and solve.
///
/// This class factors a batch of SPD matrices stored in packed lower-triangular
/// format (LAPACK/MAGMA column-major packed storage) and applies the inverse via
/// MAGMA batched triangular solves.
class MagmaPackedLowerCholesky
{
private:
int n = 0;
int batch_size = 0;
int packed_size = 0;
mutable Array<real_t *> factor_ptrs;
mutable Array<real_t *> rhs_ptrs;
Array<magma_int_t> info;
magma_queue_t queue = nullptr;
public:
MagmaPackedLowerCholesky();
void SetQueue(magma_queue_t q) { queue = q; }
int GetNumRows() const { return n; }
int GetNumMatrices() const { return batch_size; }
int GetPackedSize() const { return packed_size; }
/// Factor packed-lower matrices A into L (in-place copy then factor).
void Factor(const TriPackLowerMatrix &A,
TriPackLowerMatrix &L);
/// Solve A x = b using L from Factor(), overwriting rhs_sol with x.
void SolveInPlace(const TriPackLowerMatrix &L,
Vector &rhs_sol) const;
};
/// Workspace + operations for MAGMA packed-lower batched inverse and apply.
///
/// This class computes the inverse of a batch of SPD matrices stored in packed
/// lower-triangular format (LAPACK/MAGMA column-major packed storage) using
/// MAGMA's `ppinv_batched`. The resulting packed inverse can be applied to a
/// batch of vectors using MAGMA's packed-symmetric batched matvec when
/// available, falling back to an MFEM device kernel for larger sizes.
class MagmaPackedLowerInverse
{
private:
int n = 0;
int batch_size = 0;
int packed_size = 0;
mutable Array<real_t *> inv_ptrs;
mutable Array<real_t *> rhs_ptrs;
mutable Vector work;
Array<magma_int_t> info;
magma_queue_t queue = nullptr;
public:
MagmaPackedLowerInverse();
void SetQueue(magma_queue_t q) { queue = q; }
int GetNumRows() const { return n; }
int GetNumMatrices() const { return batch_size; }
int GetPackedSize() const { return packed_size; }
/// Compute packed inverse of A into A_inv (in-place copy then invert).
void Compute(const TriPackLowerMatrix &A,
TriPackLowerMatrix &A_inv);
/// Apply packed inverse to rhs_sol, overwriting rhs_sol with the result.
void ApplyInPlace(const TriPackLowerMatrix &A_inv,
Vector &rhs_sol) const;
};
namespace tripack
{
namespace magma
{
inline void ComputeCholeskyLower(
const TriPackLowerMatrix &packed_lower,
TriPackLowerMatrix &lower_factor,
MagmaPackedLowerCholesky &ws)
{
ws.Factor(packed_lower, lower_factor);
}
inline void SolveCholeskyLowerInPlace(
const TriPackLowerMatrix &lower_factor,
Vector &rhs_sol,
MagmaPackedLowerCholesky &ws)
{
ws.SolveInPlace(lower_factor, rhs_sol);
}
inline void ComputeInverseLower(
const TriPackLowerMatrix &packed_lower,
TriPackLowerMatrix &lower_inverse,
MagmaPackedLowerInverse &ws)
{
ws.Compute(packed_lower, lower_inverse);
}
inline void ApplyInverseLowerInPlace(
const TriPackLowerMatrix &lower_inverse,
Vector &rhs_sol,
MagmaPackedLowerInverse &ws)
{
ws.ApplyInPlace(lower_inverse, rhs_sol);
}
} // namespace magma
} // namespace tripack
} // namespace mfem
#endif // MFEM_USE_MAGMA
#endif // MFEM_TRIPACK_MAGMA
+3
View File
@@ -44,6 +44,9 @@
#endif // MFEM_USE_MPI
#include "linalg/linalg.hpp"
#ifdef MFEM_USE_MAGMA
#include "linalg/tripack_magma.hpp"
#endif
#include "mesh/mesh_headers.hpp"
+157 -1
View File
@@ -40,9 +40,15 @@
// mpirun -np 4 darcy -m ../../data/fichera-q2.mesh
#include "mfem.hpp"
#include <cstring>
#include <iostream>
#include <memory>
#ifdef MFEM_USE_UMPIRE
#include <umpire/Allocator.hpp>
#include <umpire/ResourceManager.hpp>
#endif
#include "discrete_divergence.hpp"
#include "hdiv_linear_solver.hpp"
@@ -53,6 +59,78 @@ using namespace mfem;
ParMesh LoadParMesh(const char *mesh_file, int ser_ref = 0, int par_ref = 0);
namespace
{
HdivSaddlePointSolver::L2InverseType ParseL2InverseType(const char *name)
{
if (!name || strcmp(name, "cg") == 0)
{
return HdivSaddlePointSolver::L2InverseType::CG;
}
if (strcmp(name, "magma-packed") == 0)
{
return HdivSaddlePointSolver::L2InverseType::MAGMA_PACKED;
}
if (strcmp(name, "magma-full") == 0)
{
return HdivSaddlePointSolver::L2InverseType::MAGMA_FULL;
}
MFEM_ABORT("Unknown -l2inv value: " << name
<< " (expected: cg | magma-packed | magma-full)");
return HdivSaddlePointSolver::L2InverseType::CG;
}
#ifdef MFEM_USE_UMPIRE
void ReportUmpireAllocator(const char *label, const char *alloc_name)
{
auto &rm = umpire::ResourceManager::getInstance();
if (!rm.isAllocator(alloc_name))
{
if (Mpi::Root())
{
cout << label << ": allocator '" << alloc_name
<< "' not found (no allocations yet?)\n";
}
return;
}
auto alloc = rm.getAllocator(alloc_name);
const unsigned long long cur = alloc.getCurrentSize();
const unsigned long long hwm = alloc.getHighWatermark();
unsigned long long cur_sum = 0, cur_max = 0;
unsigned long long hwm_sum = 0, hwm_max = 0;
MPI_Reduce(&cur, &cur_sum, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0,
MPI_COMM_WORLD);
MPI_Reduce(&cur, &cur_max, 1, MPI_UNSIGNED_LONG_LONG, MPI_MAX, 0,
MPI_COMM_WORLD);
MPI_Reduce(&hwm, &hwm_sum, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0,
MPI_COMM_WORLD);
MPI_Reduce(&hwm, &hwm_max, 1, MPI_UNSIGNED_LONG_LONG, MPI_MAX, 0,
MPI_COMM_WORLD);
if (Mpi::Root())
{
cout << label << " (Umpire '" << alloc_name << "'): "
<< "current(sum/max)=(" << cur_sum << "/" << cur_max << ") bytes, "
<< "hwm(sum/max)=(" << hwm_sum << "/" << hwm_max << ") bytes\n";
}
}
void ReportUmpireMemory(const char *label)
{
if (Mpi::Root()) { cout << label << '\n'; }
ReportUmpireAllocator(" host", MemoryManager::GetUmpireHostAllocatorName());
ReportUmpireAllocator(" device",
MemoryManager::GetUmpireDeviceAllocatorName());
}
#else
void ReportUmpireMemory(const char *) { }
#endif
} // namespace
int main(int argc, char *argv[])
{
Mpi::Init(argc, argv);
@@ -64,6 +142,11 @@ int main(int argc, char *argv[])
int par_ref = 1;
int order = 3;
real_t alpha = 0.0;
const char *l2inv = "cg";
bool bench_l2inv = false;
int l2inv_reps = 100;
bool use_umpire_pool = false;
bool report_umpire_mem = false;
OptionsParser args(argc, argv);
args.AddOption(&device_config, "-d", "--device",
@@ -75,8 +158,32 @@ int main(int argc, char *argv[])
"Number of times to refine the mesh in parallel.");
args.AddOption(&order, "-o", "--order", "Polynomial degree.");
args.AddOption(&alpha, "-a", "--alpha", "Value of alpha coefficient.");
args.AddOption(&l2inv, "-l2inv", "--l2-inverse",
"Local L2 mass inverse: cg | magma-packed | magma-full.");
args.AddOption(&bench_l2inv, "-l2bench", "--l2-bench",
"-no-l2bench", "--no-l2-bench",
"Benchmark the local L2 inverse apply.");
args.AddOption(&l2inv_reps, "-l2reps", "--l2-repetitions",
"Repetitions for -l2bench timing.");
args.AddOption(&use_umpire_pool, "-umpire-pool", "--umpire-pool",
"-no-umpire-pool", "--no-umpire-pool",
"Use Umpire QuickPool allocators for MFEM allocations.");
args.AddOption(&report_umpire_mem, "-mem", "--report-memory",
"-no-mem", "--no-report-memory",
"Report Umpire allocator memory usage.");
args.ParseCheck();
#ifdef MFEM_USE_UMPIRE
if (use_umpire_pool)
{
MemoryManager::SetUmpireHostAllocatorName("mfem_host_pool");
MemoryManager::SetUmpireDeviceAllocatorName("mfem_device_pool");
}
#else
MFEM_VERIFY(!use_umpire_pool, "MFEM was built without Umpire support.");
MFEM_VERIFY(!report_umpire_mem, "MFEM was built without Umpire support.");
#endif
Device device(device_config);
if (Mpi::Root()) { device.Print(); }
@@ -140,8 +247,26 @@ int main(int argc, char *argv[])
ConstantCoefficient one(1.0);
ConstantCoefficient alpha_coeff(alpha);
const auto solver_mode = HdivSaddlePointSolver::Mode::DARCY;
const auto l2inv_type = ParseL2InverseType(l2inv);
StopWatch setup_sw;
setup_sw.Start();
HdivSaddlePointSolver saddle_point_solver(
mesh, fes_rt, fes_l2, alpha_coeff, one, ess_rt_dofs, solver_mode);
mesh, fes_rt, fes_l2, alpha_coeff, one, ess_rt_dofs, solver_mode, l2inv_type);
MFEM_DEVICE_SYNC;
setup_sw.Stop();
if (Mpi::Root())
{
const int n = fes_l2.GetTypicalFE()->GetDof();
const int ne = mesh.GetNE();
const size_t full_bytes = static_cast<size_t>(ne)*n*n*sizeof(real_t);
const size_t packed_bytes = static_cast<size_t>(ne)*n*(n+1)/2*sizeof(real_t);
cout << "Setup time: " << setup_sw.RealTime() << " s\n"
<< "Local element matrices (theoretical): full=" << full_bytes
<< " bytes, packed=" << packed_bytes << " bytes\n";
}
if (report_umpire_mem) { ReportUmpireMemory("After setup"); }
const Array<int> &offsets = saddle_point_solver.GetOffsets();
BlockVector X_block(offsets), B_block(offsets);
@@ -160,6 +285,37 @@ int main(int argc, char *argv[])
<< saddle_point_solver.GetNumIterations()
<< "\nElapsed: " << tic_toc.RealTime() << endl;
}
if (report_umpire_mem) { ReportUmpireMemory("After solve"); }
if (bench_l2inv)
{
const int n_l2 = offsets[1];
Vector rhs(n_l2), x_l2(n_l2);
rhs.UseDevice(true);
x_l2.UseDevice(true);
rhs.Randomize(1);
// Warm up and time repeated applications.
saddle_point_solver.GetL2Inverse().Mult(rhs, x_l2);
MFEM_DEVICE_SYNC;
StopWatch sw;
sw.Start();
for (int r = 0; r < l2inv_reps; ++r)
{
saddle_point_solver.GetL2Inverse().Mult(rhs, x_l2);
}
MFEM_DEVICE_SYNC;
sw.Stop();
const double local_ms = 1000.0*sw.RealTime()/l2inv_reps;
double max_ms = 0.0;
MPI_Reduce(&local_ms, &max_ms, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
if (Mpi::Root())
{
cout << "L2 inverse apply (ms/apply, max over ranks): " << max_ms << '\n';
}
}
ParGridFunction x(&fes_l2);
x.SetFromTrueDofs(X_block.GetBlock(0));
+91 -1
View File
@@ -32,8 +32,13 @@
// mpirun -np 4 grad_div -sp -ams -lor -hb -m ../../data/fichera-q2.mesh -rp 0
#include "mfem.hpp"
#include <cstring>
#include <iostream>
#include <memory>
#ifdef MFEM_USE_UMPIRE
#include <umpire/Allocator.hpp>
#include <umpire/ResourceManager.hpp>
#endif
#include "hdiv_linear_solver.hpp"
#include "../solvers/lor_mms.hpp"
@@ -43,6 +48,67 @@ using namespace mfem;
ParMesh LoadParMesh(const char *mesh_file, int ser_ref = 0, int par_ref = 0);
void SolveCG(Operator &A, Solver &P, const Vector &B, Vector &X);
namespace
{
HdivSaddlePointSolver::L2InverseType ParseL2InverseType(const char *name)
{
if (!name || strcmp(name, "cg") == 0)
{
return HdivSaddlePointSolver::L2InverseType::CG;
}
if (strcmp(name, "magma-packed") == 0)
{
return HdivSaddlePointSolver::L2InverseType::MAGMA_PACKED;
}
if (strcmp(name, "magma-full") == 0)
{
return HdivSaddlePointSolver::L2InverseType::MAGMA_FULL;
}
MFEM_ABORT("Unknown -l2inv value: " << name
<< " (expected: cg | magma-packed | magma-full)");
return HdivSaddlePointSolver::L2InverseType::CG;
}
#ifdef MFEM_USE_UMPIRE
void ReportUmpireAllocator(const char *label, const char *alloc_name)
{
auto &rm = umpire::ResourceManager::getInstance();
if (!rm.isAllocator(alloc_name)) { return; }
auto alloc = rm.getAllocator(alloc_name);
const unsigned long long cur = alloc.getCurrentSize();
const unsigned long long hwm = alloc.getHighWatermark();
unsigned long long cur_sum = 0, cur_max = 0;
unsigned long long hwm_sum = 0, hwm_max = 0;
MPI_Reduce(&cur, &cur_sum, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0,
MPI_COMM_WORLD);
MPI_Reduce(&cur, &cur_max, 1, MPI_UNSIGNED_LONG_LONG, MPI_MAX, 0,
MPI_COMM_WORLD);
MPI_Reduce(&hwm, &hwm_sum, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0,
MPI_COMM_WORLD);
MPI_Reduce(&hwm, &hwm_max, 1, MPI_UNSIGNED_LONG_LONG, MPI_MAX, 0,
MPI_COMM_WORLD);
if (Mpi::Root())
{
cout << label << " (Umpire '" << alloc_name << "'): "
<< "current(sum/max)=(" << cur_sum << "/" << cur_max << ") bytes, "
<< "hwm(sum/max)=(" << hwm_sum << "/" << hwm_max << ") bytes\n";
}
}
void ReportUmpireMemory(const char *label)
{
if (Mpi::Root()) { cout << label << '\n'; }
ReportUmpireAllocator(" host", MemoryManager::GetUmpireHostAllocatorName());
ReportUmpireAllocator(" device",
MemoryManager::GetUmpireDeviceAllocatorName());
}
#else
void ReportUmpireMemory(const char *) { }
#endif
} // namespace
int main(int argc, char *argv[])
{
Mpi::Init(argc, argv);
@@ -57,6 +123,9 @@ int main(int argc, char *argv[])
bool use_ams = false;
bool use_lor_ams = false;
bool use_hybridization = false;
const char *l2inv = "cg";
bool use_umpire_pool = false;
bool report_umpire_mem = false;
OptionsParser args(argc, argv);
args.AddOption(&device_config, "-d", "--device",
@@ -77,8 +146,27 @@ int main(int argc, char *argv[])
args.AddOption(&use_hybridization,
"-hb", "--hybridization", "-no-hb", "--no-hybridization",
"Enable or disable hybridization solver.");
args.AddOption(&l2inv, "-l2inv", "--l2-inverse",
"Local L2 mass inverse: cg | magma-packed | magma-full.");
args.AddOption(&use_umpire_pool, "-umpire-pool", "--umpire-pool",
"-no-umpire-pool", "--no-umpire-pool",
"Use Umpire QuickPool allocators for MFEM allocations.");
args.AddOption(&report_umpire_mem, "-mem", "--report-memory",
"-no-mem", "--no-report-memory",
"Report Umpire allocator memory usage.");
args.ParseCheck();
#ifdef MFEM_USE_UMPIRE
if (use_umpire_pool)
{
MemoryManager::SetUmpireHostAllocatorName("mfem_host_pool");
MemoryManager::SetUmpireDeviceAllocatorName("mfem_device_pool");
}
#else
MFEM_VERIFY(!use_umpire_pool, "MFEM was built without Umpire support.");
MFEM_VERIFY(!report_umpire_mem, "MFEM was built without Umpire support.");
#endif
if (!use_saddle_point && !use_ams && !use_lor_ams && !use_hybridization)
{
if (Mpi::Root()) { cout << "No solver enabled. Exiting.\n"; }
@@ -125,9 +213,11 @@ int main(int argc, char *argv[])
L2_FECollection fec_l2(order-1, dim, b2, mt);
ParFiniteElementSpace fes_l2(&mesh, &fec_l2);
const auto l2inv_type = ParseL2InverseType(l2inv);
HdivSaddlePointSolver saddle_point_solver(
mesh, fes_rt, fes_l2, alpha_coeff, beta_coeff, ess_rt_dofs,
HdivSaddlePointSolver::Mode::GRAD_DIV);
HdivSaddlePointSolver::Mode::GRAD_DIV, l2inv_type);
if (report_umpire_mem) { ReportUmpireMemory("After saddle-point setup"); }
const Array<int> &offsets = saddle_point_solver.GetOffsets();
@@ -15,6 +15,199 @@
namespace mfem
{
namespace
{
#ifdef MFEM_USE_MAGMA
#ifdef MFEM_USE_SINGLE
#define MFEM_HDIV_MAGMA_PREFIX(stub) magma_s##stub
#define MFEM_HDIV_MAGMA_SET_POINTER magma_sset_pointer
#elif defined(MFEM_USE_DOUBLE)
#define MFEM_HDIV_MAGMA_PREFIX(stub) magma_d##stub
#define MFEM_HDIV_MAGMA_SET_POINTER magma_dset_pointer
#else
#error "Unsupported MFEM precision for MAGMA in hdiv-linear-solver."
#endif
real_t **SetMagmaPointerArray(Array<real_t *> &ptrs,
real_t *data,
const int stride,
const int batch_size,
const magma_queue_t queue)
{
if (ptrs.Size() != batch_size)
{
if (ptrs.Size() != 0) { magma_queue_sync(queue); }
ptrs.SetSize(batch_size, Device::GetDeviceMemoryType());
}
real_t **d_ptrs = ptrs.Write();
MFEM_HDIV_MAGMA_SET_POINTER(d_ptrs, data, 1, 0, 0, stride,
batch_size, queue);
return d_ptrs;
}
class MagmaPackedL2MassInverse final : public Solver
{
private:
const FiniteElementSpace &fes;
Coefficient &coeff;
const IntegrationRule &ir;
TriPackLowerMatrix L_factor;
mutable MagmaPackedLowerCholesky ws;
public:
MagmaPackedL2MassInverse(const FiniteElementSpace &fes_,
Coefficient &coeff_,
const IntegrationRule &ir_)
: Solver(fes_.GetTrueVSize()),
fes(fes_),
coeff(coeff_),
ir(ir_)
{
MFEM_VERIFY(fes.IsDGSpace(), "MagmaPackedL2MassInverse requires DG.");
MFEM_VERIFY(UsesTensorBasis(fes),
"MagmaPackedL2MassInverse requires a tensor basis.");
MFEM_VERIFY(Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK),
"MAGMA L2 inverse requires CUDA or HIP device backend.");
Update();
}
void Update()
{
MassIntegrator mass(coeff, &ir);
mass.AssembleEATriangular(fes, L_factor, false);
tripack::magma::ComputeCholeskyLower(L_factor, L_factor, ws);
}
void Mult(const Vector &b, Vector &u) const override
{
u = b;
u.UseDevice(true);
tripack::magma::SolveCholeskyLowerInPlace(L_factor, u, ws);
}
void SetOperator(const Operator &) override
{
MFEM_ABORT("SetOperator not supported with MagmaPackedL2MassInverse.");
}
};
class MagmaFullL2MassInverse final : public Solver
{
private:
const FiniteElementSpace &fes;
Coefficient &coeff;
const IntegrationRule &ir;
Vector A_factor;
int n = 0;
int batch_size = 0;
mutable Array<real_t *> mat_ptrs;
mutable Array<real_t *> rhs_ptrs;
Array<magma_int_t> info;
magma_queue_t queue = nullptr;
public:
MagmaFullL2MassInverse(const FiniteElementSpace &fes_,
Coefficient &coeff_,
const IntegrationRule &ir_)
: Solver(fes_.GetTrueVSize()),
fes(fes_),
coeff(coeff_),
ir(ir_)
{
MFEM_VERIFY(fes.IsDGSpace(), "MagmaFullL2MassInverse requires DG.");
MFEM_VERIFY(UsesTensorBasis(fes),
"MagmaFullL2MassInverse requires a tensor basis.");
MFEM_VERIFY(Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK),
"MAGMA L2 inverse requires CUDA or HIP device backend.");
queue = Magma::Queue();
Update();
}
void Update()
{
MassIntegrator mass(coeff, &ir);
n = fes.GetTypicalFE()->GetDof();
batch_size = fes.GetMesh()->GetNE();
// MassIntegrator::AssembleEA expects the output Vector to be sized by
// the caller (unlike AssembleEATriangular which sizes its output).
A_factor.SetSize(batch_size*n*n, Device::GetDeviceMemoryType());
A_factor.UseDevice(true);
mass.AssembleEA(fes, A_factor, false);
MFEM_VERIFY(A_factor.Size() == batch_size*n*n,
"Unexpected element matrix storage size.");
if (batch_size == 0) { return; }
real_t *A_data = A_factor.ReadWrite();
real_t **dA =
SetMagmaPointerArray(mat_ptrs, A_data, n*n, batch_size, queue);
info.SetSize(batch_size, Device::GetDeviceMemoryType());
magma_int_t *d_info = info.Write();
magma_memset(d_info, 0, batch_size*sizeof(magma_int_t));
const magma_int_t status =
MFEM_HDIV_MAGMA_PREFIX(potrf_batched)(
MagmaLower, n, dA, n, d_info, batch_size, queue);
MFEM_VERIFY(status == MAGMA_SUCCESS, "MAGMA full potrf batched failed.");
magma_queue_sync(queue);
const magma_int_t *h_info = info.HostRead();
for (int e = 0; e < batch_size; ++e)
{
MFEM_VERIFY(h_info[e] == 0,
"MAGMA full potrf failed on matrix " << e << '.');
}
}
void Mult(const Vector &b, Vector &u) const override
{
MFEM_VERIFY(queue != nullptr, "MAGMA queue is not set.");
MFEM_VERIFY(b.Size() == height, "Invalid RHS size.");
u = b;
u.UseDevice(true);
if (batch_size == 0) { return; }
real_t *A_data = const_cast<real_t *>(A_factor.Read());
real_t **dA =
SetMagmaPointerArray(mat_ptrs, A_data, n*n, batch_size, queue);
real_t *rhs_data = u.ReadWrite();
real_t **dB =
SetMagmaPointerArray(rhs_ptrs, rhs_data, n, batch_size, queue);
const magma_int_t status =
MFEM_HDIV_MAGMA_PREFIX(potrs_batched)(
MagmaLower, n, 1, dA, n, dB, n, batch_size, queue);
MFEM_VERIFY(status == MAGMA_SUCCESS, "MAGMA full potrs batched failed.");
}
void SetOperator(const Operator &) override
{
MFEM_ABORT("SetOperator not supported with MagmaFullL2MassInverse.");
}
};
#undef MFEM_HDIV_MAGMA_SET_POINTER
#undef MFEM_HDIV_MAGMA_PREFIX
#endif // MFEM_USE_MAGMA
} // namespace
/// Replace x[i] with 1.0/x[i] for all i.
void Reciprocal(Vector &x)
{
@@ -65,7 +258,7 @@ const IntegrationRule &GetMassIntRule(FiniteElementSpace &fes_l2)
HdivSaddlePointSolver::HdivSaddlePointSolver(
ParMesh &mesh, ParFiniteElementSpace &fes_rt_, ParFiniteElementSpace &fes_l2_,
Coefficient &L_coeff_, Coefficient &R_coeff_, const Array<int> &ess_rt_dofs_,
Mode mode_)
Mode mode_, L2InverseType l2_inv_type_)
: minres(mesh.GetComm()),
order(fes_rt_.GetMaxElementOrder()),
fec_l2(order - 1, mesh.Dimension(), b2, mt),
@@ -81,6 +274,7 @@ HdivSaddlePointSolver::HdivSaddlePointSolver(
L_coeff(L_coeff_),
R_coeff(R_coeff_),
mode(mode_),
l2_inv_type(l2_inv_type_),
qs(mesh, GetMassIntRule(fes_l2)),
W_coeff_qf(qs),
W_mix_coeff_qf(qs),
@@ -154,9 +348,10 @@ HdivSaddlePointSolver::HdivSaddlePointSolver(
HdivSaddlePointSolver::HdivSaddlePointSolver(
ParMesh &mesh, ParFiniteElementSpace &fes_rt_, ParFiniteElementSpace &fes_l2_,
Coefficient &R_coeff_, const Array<int> &ess_rt_dofs_)
Coefficient &R_coeff_, const Array<int> &ess_rt_dofs_,
L2InverseType l2_inv_type_)
: HdivSaddlePointSolver(mesh, fes_rt_, fes_l2_, zero, R_coeff_,
ess_rt_dofs_, Mode::DARCY)
ess_rt_dofs_, Mode::DARCY, l2_inv_type_)
{ }
void HdivSaddlePointSolver::Setup()
@@ -189,7 +384,28 @@ void HdivSaddlePointSolver::Setup()
});
}
L_inv.reset(new DGMassInverse(fes_l2, W_mix_coeff));
switch (l2_inv_type)
{
case L2InverseType::CG:
L_inv.reset(new DGMassInverse(fes_l2, W_mix_coeff));
break;
case L2InverseType::MAGMA_PACKED:
#ifdef MFEM_USE_MAGMA
L_inv.reset(new MagmaPackedL2MassInverse(fes_l2, W_mix_coeff,
qs.GetIntRule(0)));
#else
MFEM_ABORT("MFEM was built without MAGMA support.");
#endif
break;
case L2InverseType::MAGMA_FULL:
#ifdef MFEM_USE_MAGMA
L_inv.reset(new MagmaFullL2MassInverse(fes_l2, W_mix_coeff,
qs.GetIntRule(0)));
#else
MFEM_ABORT("MFEM was built without MAGMA support.");
#endif
break;
}
if (zero_l2_block)
{
@@ -33,6 +33,14 @@ public:
GRAD_DIV, ///< Grad-div problem.
DARCY ///< Darcy/mixed Poisson problem.
};
/// How to apply the local DG L2 mass inverse.
enum class L2InverseType
{
CG, ///< DGMassInverse (default).
MAGMA_PACKED, ///< MAGMA packed-lower batched Cholesky+solve.
MAGMA_FULL ///< MAGMA full dense batched Cholesky+solve.
};
private:
MINRESSolver minres;
@@ -63,7 +71,7 @@ private:
// Components needed for the block operator
OperatorHandle L, R, R_e; ///< Mass matrices.
std::unique_ptr<HypreParMatrix> D, Dt, D_e; ///< Divergence matrices.
std::shared_ptr<DGMassInverse> L_inv; ///< Inverse of the DG mass matrix.
std::shared_ptr<Solver> L_inv; ///< Inverse of the DG mass matrix.
std::shared_ptr<Operator> A_11; ///< (1,1)-block of the matrix
/// Diagonals of the mass matrices
@@ -85,6 +93,7 @@ private:
Coefficient &L_coeff, &R_coeff;
const Mode mode;
const L2InverseType l2_inv_type;
bool zero_l2_block = false;
QuadratureSpace qs;
QuadratureFunction W_coeff_qf, W_mix_coeff_qf;
@@ -135,7 +144,8 @@ public:
Coefficient &L_coeff_,
Coefficient &R_coeff_,
const Array<int> &ess_rt_dofs_,
Mode mode_);
Mode mode_,
L2InverseType l2_inv_type_ = L2InverseType::CG);
/// @brief Creates a linear solver for the case when the L2 diagonal block is
/// zero (for Darcy problems).
@@ -146,7 +156,8 @@ public:
ParFiniteElementSpace &fes_rt_,
ParFiniteElementSpace &fes_l2_,
Coefficient &R_coeff_,
const Array<int> &ess_rt_dofs_);
const Array<int> &ess_rt_dofs_,
L2InverseType l2_inv_type_ = L2InverseType::CG);
/// @brief Build the linear operator and solver. Must be called when the
/// coefficients change.
@@ -168,6 +179,12 @@ public:
const Array<int> &GetOffsets() const { return offsets; }
/// Returns the internal MINRES solver.
MINRESSolver &GetMINRES() { return minres; }
/// Returns the internal L2 DG mass inverse operator.
const Solver &GetL2Inverse() const { return *L_inv; }
/// Returns the selected local L2 inverse method.
L2InverseType GetL2InverseType() const { return l2_inv_type; }
};
} // namespace mfem
+16 -100
View File
@@ -33,7 +33,6 @@
//
// Sample runs: lor-transfer
// lor-transfer -h1
// lor-transfer -ea -w
// lor-transfer -t
// lor-transfer -m ../../data/star-q2.mesh -lref 5 -p 4
// lor-transfer -m ../../data/star-mixed.mesh -lref 3 -p 2
@@ -60,12 +59,11 @@ string direction;
// Exact functions to project
real_t RHO_exact(const Vector &x);
real_t W_exact(const Vector &x);
real_t weight(const Vector &x);
// Helper functions
void visualize(VisItDataCollection &, string, int, int, int visport = 19916);
real_t compute_mass(GridFunction &, real_t, string, CoefficientWithOrder);
real_t compute_mass(FiniteElementSpace *, real_t, VisItDataCollection &,
string);
int main(int argc, char *argv[])
{
@@ -78,7 +76,6 @@ int main(int argc, char *argv[])
bool useH1 = false;
int visport = 19916;
bool use_pointwise_transfer = false;
bool use_weighted_transfer = false;
const char *device_config = "cpu";
bool use_ea = false;
@@ -101,9 +98,6 @@ int main(int argc, char *argv[])
args.AddOption(&use_pointwise_transfer, "-t", "--use-pointwise-transfer",
"-no-t", "--dont-use-pointwise-transfer",
"Use pointwise transfer operators instead of L2 projection.");
args.AddOption(&use_weighted_transfer, "-w", "--use-weighted-transfer",
"-no-w", "--dont-use-weighted-transfer",
"Use coefficient-weighted L2 projection.");
args.AddOption(&device_config, "-d", "--device",
"Device configuration string, see Device::Configure().");
args.AddOption(&use_ea, "-ea", "--ea-version", "-no-ea",
@@ -113,15 +107,6 @@ int main(int argc, char *argv[])
// Configure device
Device device(device_config);
if (use_weighted_transfer && !use_pointwise_transfer)
{
if (problem != 5)
{
cout << "Switching to positive problem = 5 for weighted transfer.\n";
}
problem = 5;
}
// Read the mesh from the given mesh file.
Mesh mesh(mesh_file, 1, 1);
int dim = mesh.Dimension();
@@ -153,14 +138,6 @@ int main(int argc, char *argv[])
FiniteElementSpace fespace(&mesh, fec);
FiniteElementSpace fespace_lor(&mesh_lor, fec_lor);
FunctionCoefficient weight_fn_coeff(weight);
CoefficientWithOrder weight_coeff;
if (use_weighted_transfer)
{
weight_coeff.coeff = &weight_fn_coeff;
weight_coeff.order = 2;
}
GridFunction rho(&fespace);
GridFunction rho_lor(&fespace_lor);
@@ -188,7 +165,7 @@ int main(int argc, char *argv[])
rho.SetTrueVector();
rho.SetFromTrueVector();
real_t ho_mass = compute_mass(rho, -1.0, "HO ", weight_coeff);
real_t ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO ");
if (vis) { visualize(HO_dc, "HO", Wx, Wy, visport); Wx += offx; }
GridTransfer *gt;
@@ -198,8 +175,7 @@ int main(int argc, char *argv[])
}
else
{
gt = new L2ProjectionGridTransfer(fespace, fespace_lor, weight_coeff,
weight_coeff);
gt = new L2ProjectionGridTransfer(fespace, fespace_lor);
}
// Configure element assembly for device acceleration
@@ -210,44 +186,9 @@ int main(int argc, char *argv[])
// HO->LOR restriction
direction = "HO -> LOR @ LOR";
R.Mult(rho, rho_lor);
compute_mass(rho_lor, ho_mass, "R(HO) ", weight_coeff);
compute_mass(&fespace_lor, ho_mass, LOR_dc, "R(HO) ");
if (vis) { visualize(LOR_dc, "R(HO)", Wx, Wy, visport); Wx += offx; }
if (use_weighted_transfer && !use_pointwise_transfer)
{
// Transfer velocity while conserving rho-weighted momentum.
GridFunctionCoefficient rho_coeff(&rho);
GridFunctionCoefficient rho_lor_coeff(&rho_lor);
ProductCoefficient prod_coeff(weight_fn_coeff, rho_coeff);
ProductCoefficient prod_lor_coeff(weight_fn_coeff, rho_lor_coeff);
CoefficientWithOrder prod_weight(prod_coeff, order + 2);
CoefficientWithOrder prod_lor_weight(prod_lor_coeff, lorder + 2);
GridFunction w(&fespace), w_lor(&fespace_lor);
FunctionCoefficient W(W_exact);
w.ProjectCoefficient(W);
cout << '\n';
const real_t ho_momentum = compute_mass(w, -1.0, "rho w HO ", prod_weight);
L2ProjectionGridTransfer vel_gt(fespace, fespace_lor, prod_weight,
prod_lor_weight);
vel_gt.UseEA(use_ea);
vel_gt.ForwardOperator().Mult(w, w_lor);
compute_mass(w_lor, ho_momentum, "rho w LOR", prod_lor_weight);
if (vel_gt.SupportsBackwardsOperator())
{
GridFunction w_prev = w;
vel_gt.BackwardOperator().Mult(w_lor, w);
compute_mass(w, ho_momentum, "P(rho w) ", prod_weight);
w_prev -= w;
cout.precision(12);
cout << "|w - P(R(w))|_∞ = " << w_prev.Normlinf() << "\n\n";
}
}
if (gt->SupportsBackwardsOperator())
{
const Operator &P = gt->BackwardOperator();
@@ -255,7 +196,7 @@ int main(int argc, char *argv[])
direction = "HO -> LOR @ HO";
GridFunction rho_prev = rho;
P.Mult(rho_lor, rho);
compute_mass(rho, ho_mass, "P(R(HO)) ", weight_coeff);
compute_mass(&fespace, ho_mass, HO_dc, "P(R(HO)) ");
if (vis) { visualize(HO_dc, "P(R(HO))", Wx, Wy, visport); Wx = 0; Wy += offy; }
rho_prev -= rho;
@@ -277,7 +218,7 @@ int main(int argc, char *argv[])
direction = "LOR -> HO @ LOR";
rho_lor.ProjectCoefficient(RHO);
GridFunction rho_lor_prev = rho_lor;
real_t lor_mass = compute_mass(rho_lor, -1.0, "LOR ", weight_coeff);
real_t lor_mass = compute_mass(&fespace_lor, -1.0, LOR_dc, "LOR ");
if (vis) { visualize(LOR_dc, "LOR", Wx, Wy, visport); Wx += offx; }
if (gt->SupportsBackwardsOperator())
@@ -286,14 +227,14 @@ int main(int argc, char *argv[])
// Prolongate to HO space
direction = "LOR -> HO @ HO";
P.Mult(rho_lor, rho);
compute_mass(rho, lor_mass, "P(LOR) ", weight_coeff);
compute_mass(&fespace, lor_mass, HO_dc, "P(LOR) ");
if (vis) { visualize(HO_dc, "P(LOR)", Wx, Wy, visport); Wx += offx; }
// Restrict back to LOR space. This won't give the original function because
// the rho_lor doesn't necessarily live in the range of R.
direction = "LOR -> HO @ LOR";
R.Mult(rho, rho_lor);
compute_mass(rho_lor, lor_mass, "R(P(LOR))", weight_coeff);
compute_mass(&fespace_lor, lor_mass, LOR_dc, "R(P(LOR))");
if (vis) { visualize(LOR_dc, "R(P(LOR))", Wx, Wy, visport); }
rho_lor_prev -= rho_lor;
@@ -329,26 +270,12 @@ real_t RHO_exact(const Vector &x)
return M_PI/2-atan(5*(2*x.Norml2()-1));
case 4: // basis function
return (x.Norml2() < 0.1) ? 1 : 0;
case 5: // positive function
return 2.0 + 2*x(0)*x(0) + 3*x(1)*x(1) - x(0)*x(1) + 0.1*sin(x.Norml2());
default:
return 1.0;
}
}
real_t W_exact(const Vector &x)
{
return x(1) + 0.25*cos(2*M_PI*x.Norml2());
}
real_t weight(const Vector &x)
{
return x(0)*x(0) + x(1)*x(1) + 1.0;
}
void visualize(VisItDataCollection &dc, string prefix, int x, int y,
int visport)
{
@@ -365,32 +292,21 @@ void visualize(VisItDataCollection &dc, string prefix, int x, int y,
}
real_t compute_mass(GridFunction &gf, real_t oldmass, string prefix,
CoefficientWithOrder mass_coeff)
real_t compute_mass(FiniteElementSpace *L2, real_t massL2,
VisItDataCollection &dc, string prefix)
{
FiniteElementSpace &fes = *gf.FESpace();
Mesh &mesh = *fes.GetMesh();
// Integration order is a * (element order) + b.
const int a = 2;
const int b = mesh.GetTypicalElementTransformation()->OrderW() +
mass_coeff.order;
ConstantCoefficient one(1.0);
Coefficient &coeff = mass_coeff ? *mass_coeff.coeff : one;
DomainLFIntegrator *integ = new DomainLFIntegrator(coeff, a, b);
LinearForm lf(&fes);
lf.AddDomainIntegrator(integ);
LinearForm lf(L2);
lf.AddDomainIntegrator(new DomainLFIntegrator(one));
lf.Assemble();
const real_t newmass = lf(gf);
real_t newmass = lf(*dc.GetField("density"));
cout.precision(18);
cout << space << " " << prefix << " mass = " << newmass;
if (oldmass >= 0)
if (massL2 >= 0)
{
cout.precision(4);
cout << " (" << fabs(newmass-oldmass)*100/oldmass << "%)";
cout << " (" << fabs(newmass-massL2)*100/massL2 << "%)";
}
cout << endl;
return newmass;
+16 -106
View File
@@ -33,7 +33,6 @@
//
// Sample runs: plor-transfer
// plor-transfer -h1
// plor-transfer -ea -w
// plor-transfer -t
// plor-transfer -m ../../data/star-q2.mesh -lref 5 -p 4
// plor-transfer -m ../../data/star-mixed.mesh -lref 3 -p 2
@@ -60,12 +59,11 @@ string direction;
// Exact functions to project
real_t RHO_exact(const Vector &x);
real_t W_exact(const Vector &x);
real_t weight(const Vector &x);
// Helper functions
void visualize(VisItDataCollection &, string, int, int, int /* visport */);
real_t compute_mass(ParGridFunction &, real_t, string, CoefficientWithOrder);
real_t compute_mass(ParFiniteElementSpace *, real_t, VisItDataCollection &,
string);
int main(int argc, char *argv[])
{
@@ -82,7 +80,6 @@ int main(int argc, char *argv[])
bool useH1 = false;
int visport = 19916;
bool use_pointwise_transfer = false;
bool use_weighted_transfer = false;
const char *device_config = "cpu";
bool use_ea = false;
@@ -105,9 +102,6 @@ int main(int argc, char *argv[])
args.AddOption(&use_pointwise_transfer, "-t", "--use-pointwise-transfer",
"-no-t", "--dont-use-pointwise-transfer",
"Use pointwise transfer operators instead of L2 projection.");
args.AddOption(&use_weighted_transfer, "-w", "--use-weighted-transfer",
"-no-w", "--dont-use-weighted-transfer",
"Use coefficient-weighted L2 projection.");
args.AddOption(&device_config, "-d", "--device",
"Device configuration string, see Device::Configure().");
args.AddOption(&use_ea, "-ea", "--ea-version", "-no-ea",
@@ -118,15 +112,6 @@ int main(int argc, char *argv[])
Device device(device_config);
if (Mpi::Root()) { device.Print(); }
if (use_weighted_transfer && !use_pointwise_transfer)
{
if (problem != 5 && Mpi::Root())
{
cout << "Switching to positive problem = 5 for weighted transfer.\n";
}
problem = 5;
}
// Read the mesh from the given mesh file.
Mesh serial_mesh(mesh_file, 1, 1);
ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
@@ -169,14 +154,6 @@ int main(int argc, char *argv[])
ParFiniteElementSpace fespace(&mesh, fec);
ParFiniteElementSpace fespace_lor(&mesh_lor, fec_lor);
FunctionCoefficient weight_fn_coeff(weight);
CoefficientWithOrder weight_coeff;
if (use_weighted_transfer)
{
weight_coeff.coeff = &weight_fn_coeff;
weight_coeff.order = 2;
}
ParGridFunction rho(&fespace);
ParGridFunction rho_lor(&fespace_lor);
@@ -206,7 +183,7 @@ int main(int argc, char *argv[])
rho.SetTrueVector();
rho.SetFromTrueVector();
real_t ho_mass = compute_mass(rho, -1.0, "HO ", weight_coeff);
real_t ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO ");
if (vis) { visualize(HO_dc, "HO", Wx, Wy, visport); Wx += offx; }
GridTransfer *gt;
@@ -216,8 +193,7 @@ int main(int argc, char *argv[])
}
else
{
gt = new L2ProjectionGridTransfer(fespace, fespace_lor, weight_coeff,
weight_coeff);
gt = new L2ProjectionGridTransfer(fespace, fespace_lor);
}
// Configure element assembly for device acceleration
@@ -228,7 +204,7 @@ int main(int argc, char *argv[])
// HO->LOR restriction
direction = "HO -> LOR @ LOR";
R.Mult(rho, rho_lor);
compute_mass(rho_lor, ho_mass, "R(HO) ", weight_coeff);
compute_mass(&fespace_lor, ho_mass, LOR_dc, "R(HO) ");
if (vis) { visualize(LOR_dc, "R(HO)", Wx, Wy, visport); Wx += offx; }
auto global_max = [](const Vector& v)
{
@@ -238,47 +214,6 @@ int main(int argc, char *argv[])
return max;
};
if (use_weighted_transfer && !use_pointwise_transfer)
{
// Transfer velocity while conserving rho-weighted momentum.
GridFunctionCoefficient rho_coeff(&rho);
GridFunctionCoefficient rho_lor_coeff(&rho_lor);
ProductCoefficient prod_coeff(weight_fn_coeff, rho_coeff);
ProductCoefficient prod_lor_coeff(weight_fn_coeff, rho_lor_coeff);
CoefficientWithOrder prod_weight(prod_coeff, order + 2);
CoefficientWithOrder prod_lor_weight(prod_lor_coeff, lorder + 2);
ParGridFunction w(&fespace), w_lor(&fespace_lor);
FunctionCoefficient W(W_exact);
w.ProjectCoefficient(W);
if (Mpi::Root()) { cout << '\n'; }
const real_t ho_momentum = compute_mass(w, -1.0, "rho w HO ", prod_weight);
L2ProjectionGridTransfer vel_gt(fespace, fespace_lor, prod_weight,
prod_lor_weight);
vel_gt.UseEA(use_ea);
vel_gt.ForwardOperator().Mult(w, w_lor);
compute_mass(w_lor, ho_momentum, "rho w LOR", prod_lor_weight);
if (vel_gt.SupportsBackwardsOperator())
{
ParGridFunction w_prev = w;
vel_gt.BackwardOperator().Mult(w_lor, w);
compute_mass(w, ho_momentum, "P(rho w) ", prod_weight);
w_prev -= w;
Vector w_prev_true(fespace.GetTrueVSize());
w_prev.GetTrueDofs(w_prev_true);
const real_t l_inf = global_max(w_prev_true);
if (Mpi::Root())
{
cout.precision(12);
cout << "|w - P(R(w))|_∞ = " << l_inf << "\n\n";
}
}
}
if (gt->SupportsBackwardsOperator())
{
const Operator &P = gt->BackwardOperator();
@@ -286,7 +221,7 @@ int main(int argc, char *argv[])
direction = "HO -> LOR @ HO";
ParGridFunction rho_prev = rho;
P.Mult(rho_lor, rho);
compute_mass(rho, ho_mass, "P(R(HO)) ", weight_coeff);
compute_mass(&fespace, ho_mass, HO_dc, "P(R(HO)) ");
if (vis) { visualize(HO_dc, "P(R(HO))", Wx, Wy, visport); Wx = 0; Wy += offy; }
rho_prev -= rho;
@@ -328,7 +263,7 @@ int main(int argc, char *argv[])
direction = "LOR -> HO @ LOR";
rho_lor.ProjectCoefficient(RHO);
ParGridFunction rho_lor_prev = rho_lor;
real_t lor_mass = compute_mass(rho_lor, -1.0, "LOR ", weight_coeff);
real_t lor_mass = compute_mass(&fespace_lor, -1.0, LOR_dc, "LOR ");
if (vis) { visualize(LOR_dc, "LOR", Wx, Wy, visport); Wx += offx; }
if (gt->SupportsBackwardsOperator())
@@ -337,14 +272,14 @@ int main(int argc, char *argv[])
// Prolongate to HO space
direction = "LOR -> HO @ HO";
P.Mult(rho_lor, rho);
compute_mass(rho, lor_mass, "P(LOR) ", weight_coeff);
compute_mass(&fespace, lor_mass, HO_dc, "P(LOR) ");
if (vis) { visualize(HO_dc, "P(LOR)", Wx, Wy, visport); Wx += offx; }
// Restrict back to LOR space. This won't give the original function because
// the rho_lor doesn't necessarily live in the range of R.
direction = "LOR -> HO @ LOR";
R.Mult(rho, rho_lor);
compute_mass(rho_lor, lor_mass, "R(P(LOR))", weight_coeff);
compute_mass(&fespace_lor, lor_mass, LOR_dc, "R(P(LOR))");
if (vis) { visualize(LOR_dc, "R(P(LOR))", Wx, Wy, visport); }
rho_lor_prev -= rho_lor;
@@ -399,26 +334,12 @@ real_t RHO_exact(const Vector &x)
return M_PI/2-atan(5*(2*x.Norml2()-1));
case 4: // basis function
return (x.Norml2() < 0.1) ? 1 : 0;
case 5: // positive function
return 2.0 + 2*x(0)*x(0) + 3*x(1)*x(1) - x(0)*x(1) + 0.1*sin(x.Norml2());
default:
return 1.0;
}
}
real_t W_exact(const Vector &x)
{
return x(1) + 0.25*cos(2*M_PI*x.Norml2());
}
real_t weight(const Vector &x)
{
return x(0)*x(0) + x(1)*x(1) + 1.0;
}
void visualize(VisItDataCollection &dc, string prefix, int x, int y,
int visport)
{
@@ -437,34 +358,23 @@ void visualize(VisItDataCollection &dc, string prefix, int x, int y,
}
real_t compute_mass(ParGridFunction &gf, real_t oldmass, string prefix,
CoefficientWithOrder mass_coeff)
real_t compute_mass(ParFiniteElementSpace *L2, real_t massL2,
VisItDataCollection &dc, string prefix)
{
ParFiniteElementSpace &fes = *gf.ParFESpace();
Mesh &mesh = *fes.GetMesh();
// Integration order is a * (element order) + b.
const int a = 2;
const int b = mesh.GetTypicalElementTransformation()->OrderW() +
mass_coeff.order;
ConstantCoefficient one(1.0);
Coefficient &coeff = mass_coeff ? *mass_coeff.coeff : one;
DomainLFIntegrator *integ = new DomainLFIntegrator(coeff, a, b);
ParLinearForm lf(&fes);
lf.AddDomainIntegrator(integ);
ParLinearForm lf(L2);
lf.AddDomainIntegrator(new DomainLFIntegrator(one));
lf.Assemble();
const real_t newmass = lf(gf);
real_t newmass = lf(*dc.GetParField("density"));
if (Mpi::Root())
{
cout.precision(18);
cout << space << " " << prefix << " mass = " << newmass;
if (oldmass >= 0)
if (massL2 >= 0)
{
cout.precision(4);
cout << " (" << fabs(newmass-oldmass)*100/oldmass << "%)";
cout << " (" << fabs(newmass-massL2)*100/massL2 << "%)";
}
cout << endl;
}
+2 -1
View File
@@ -74,6 +74,8 @@ set(UNIT_TESTS_SRCS
linalg/test_particlevector.cpp
linalg/test_petsc_nonlinear.cpp
linalg/test_sparsesmoothers.cpp
linalg/test_tripack.cpp
linalg/test_tripack_magma.cpp
linalg/test_vector.cpp
mesh/mesh_test_utils.cpp
mesh/test_exodus_reader.cpp
@@ -97,7 +99,6 @@ set(UNIT_TESTS_SRCS
fem/test_2d_bilininteg.cpp
fem/test_3d_bilininteg.cpp
fem/test_assembly_levels.cpp
fem/test_bdr_edgedof.cpp
fem/test_bilinearform.cpp
fem/test_block_operators.cpp
fem/test_blocknonlinearform.cpp
-705
View File
@@ -1,705 +0,0 @@
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#include "unit_tests.hpp"
#include "mfem.hpp"
#include "../mesh/mesh_test_utils.hpp"
#include <set>
#include <unordered_set>
#include <vector>
using namespace mfem;
#ifdef MFEM_USE_MPI
TEST_CASE("BoundaryEdgeDOFsPartitionInvariant",
"[Parallel][ParMesh][BoundaryEdgeDOFs]")
{
constexpr int orientation = 3;
constexpr int order = 1;
// Use all available MPI processes for partitioning
const int test_num_procs = Mpi::WorldSize();
// Create base mesh
Mesh base_mesh = OrientedTriFaceMesh(orientation, true);
base_mesh.UniformRefinement();
const int n_elements = base_mesh.GetNE();
// Use a small set of representative partitionings
std::vector<std::vector<int>> all_partitionings;
// 1. All elements on rank 0
all_partitionings.push_back(std::vector<int>(n_elements, 0));
if (test_num_procs > 1)
{
// 2. Block partition: first half on rank 0, second half on last rank
std::vector<int> &block = all_partitionings.emplace_back(n_elements);
for (int i = 0; i < n_elements; i++)
{
block[i] = (i < n_elements/2) ? 0 : test_num_procs-1;
}
// 3. Round-robin partition: elements assigned cyclically to all ranks
std::vector<int> &round_robin = all_partitionings.emplace_back(n_elements);
for (int i = 0; i < n_elements; i++)
{
round_robin[i] = i % test_num_procs;
}
}
// Create reusable FEC
ND_FECollection fec(order, 3);
std::vector<int> all_results;
all_results.reserve(all_partitionings.size());
// Test each partitioning
for (const auto& partition : all_partitionings)
{
// Create parallel mesh with current partitioning
Mesh test_mesh = OrientedTriFaceMesh(orientation, true);
test_mesh.UniformRefinement();
// For single process, use default partitioning; for multiple, use custom partition
ParMesh pmesh = (test_num_procs == 1) ?
ParMesh(MPI_COMM_WORLD, test_mesh) :
ParMesh(MPI_COMM_WORLD, test_mesh, partition.data());
// Create finite element space
ParFiniteElementSpace fespace(&pmesh, &fec);
// Extract boundary edge DOFs
Array<int> ess_tdof_list;
Array<int> boundary_edge_ldofs;
std::vector<Array<int>> attr_to_elements;
// Select the shared face to be the tested boundary
int bdr_attr = pmesh.bdr_attributes.Max();
Array<int> bdr_attrs(1);
bdr_attrs[0] = bdr_attr;
fespace.GetBoundaryElementsByAttribute(bdr_attrs, attr_to_elements);
Array<int> boundary_elements = attr_to_elements[0];
Array<int> dof_edges, dof_boundary_elements, ess_edge_list;
fespace.GetBoundaryLoopEdgeDofs(boundary_elements, ess_tdof_list,
boundary_edge_ldofs, nullptr, &dof_edges,
&dof_boundary_elements, &ess_edge_list);
// Collect total boundary edge DOFs
int local_dofs = boundary_edge_ldofs.Size();
int total_dofs;
MPI_Allreduce(&local_dofs, &total_dofs, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
all_results.push_back(total_dofs);
}
// The set of boundary edge DOFs is a property of the mesh geometry and must
// not depend on how the elements are distributed across ranks. Each result
// is the global count of selected boundary edge DOFs for one partitioning, so
// if the method correctly removes the artificial edges introduced at
// processor boundaries, every partitioning yields the same total. A mismatch
// means some partition kept or dropped a DOF that another did not.
REQUIRE(!all_results.empty());
// One refinement splits the triangular face into four sub-triangles. Its
// perimeter has six loop edges (order-1 ND: one DOF per edge); the three
// interior edges of the middle sub-triangle are shared and correctly dropped.
constexpr int expected = 6;
for (int result : all_results)
{
REQUIRE(result == expected);
}
}
TEST_CASE("BoundaryEdgeDOFsBasicFunctionality",
"[Parallel][ParMesh][BoundaryEdgeDOFs]")
{
const int orientation = GENERATE(1, 3, 5);
const int order = GENERATE(1, 2);
CAPTURE(orientation, order);
// Create test mesh
Mesh mesh = OrientedTriFaceMesh(orientation, true);
mesh.UniformRefinement();
ParMesh pmesh(MPI_COMM_WORLD, mesh);
// Create finite element space
ND_FECollection fec(order, 3);
ParFiniteElementSpace fespace(&pmesh, &fec);
// Test boundary edge DOF extraction
Array<int> ess_tdof_list;
Array<int> boundary_edge_ldofs;
Array<int> ldof_marker;
std::vector<Array<int>> attr_to_elements;
// Get boundary elements for the shared face
int bdr_attr = pmesh.bdr_attributes.Max();
Array<int> bdr_attrs(1);
bdr_attrs[0] = bdr_attr;
fespace.GetBoundaryElementsByAttribute(bdr_attrs, attr_to_elements);
Array<int> boundary_elements = attr_to_elements[0];
Array<int> dof_edges, dof_boundary_elements, ess_edge_list;
fespace.GetBoundaryLoopEdgeDofs(boundary_elements, ess_tdof_list,
boundary_edge_ldofs, &ldof_marker, &dof_edges,
&dof_boundary_elements, &ess_edge_list);
// Basic validation
REQUIRE(ldof_marker.Size() == fespace.GetVSize());
REQUIRE(ess_tdof_list.Size() >= 0);
// The output arrays share a single indexing, so they must have equal size.
REQUIRE(boundary_edge_ldofs.Size() == dof_edges.Size());
REQUIRE(dof_edges.Size() == dof_boundary_elements.Size());
// Verify all boundary edge DOFs are marked in ldof_marker
for (int dof : boundary_edge_ldofs)
{
REQUIRE(ldof_marker[dof] == 1);
}
}
// Helper function to compute boundary loop length
real_t ComputeBoundaryLoopLength(ParMesh* pmesh, const Array<int>& dof_edges)
{
real_t local_length = 0.0;
std::unordered_set<int> processed_edges;
for (int i = 0; i < dof_edges.Size(); i++)
{
int edge_id = dof_edges[i];
if (!processed_edges.insert(edge_id).second) { continue; }
Array<int> edge_verts;
pmesh->GetEdgeVertices(edge_id, edge_verts);
const real_t* v0 = pmesh->GetVertex(edge_verts[0]);
const real_t* v1 = pmesh->GetVertex(edge_verts[1]);
real_t edge_length = 0.0;
for (int d = 0; d < pmesh->SpaceDimension(); d++)
{
real_t diff = v1[d] - v0[d];
edge_length += diff * diff;
}
local_length += sqrt(edge_length);
}
return local_length;
}
TEST_CASE("BoundaryEdgeDOFsNestedCubes",
"[Parallel][ParMesh][BoundaryEdgeDOFs]")
{
const int order = GENERATE(1, 2);
// Expected processor-invariant results for nested cubes mesh (1 refinement)
// order=1: 16 tdofs, sum=16.0, length=2.0
// order=2: 32 tdofs, sum=32.0, length=2.0
int exp_tdofs = (order == 1) ? 16 : 32;
real_t exp_sum = (order == 1) ? real_t(16.0) : real_t(32.0);
real_t exp_length = real_t(2.0);
struct BoundaryTest
{
int attr_value;
Vector normal;
std::string name;
};
std::vector<BoundaryTest> boundary_tests =
{
{7, Vector({0, 0, -1}), "-z"},
{8, Vector({0, 0, 1}), "+z"},
{9, Vector({0, -1, 0}), "-y"},
{10, Vector({1, 0, 0}), "+x"},
{11, Vector({0, 1, 0}), "+y"},
{12, Vector({-1, 0, 0}), "-x"}
};
const char* mesh_file = "../../data/nested_cubes.msh";
Mesh mesh(mesh_file, 1, 1);
mesh.UniformRefinement();
ParMesh pmesh(MPI_COMM_WORLD, mesh);
ND_FECollection fec(order, 3);
ParFiniteElementSpace fespace(&pmesh, &fec);
int num_procs;
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
for (const auto& test : boundary_tests)
{
CAPTURE(test.name, test.attr_value, order, num_procs);
std::vector<Array<int>> attr_to_elements;
Array<int> bdr_attrs(1);
bdr_attrs[0] = test.attr_value;
fespace.GetBoundaryElementsByAttribute(bdr_attrs, attr_to_elements);
Array<int> boundary_elements = attr_to_elements[0];
Array<int> ess_tdof_list;
Array<int> ldof_marker;
Array<int> boundary_edge_ldofs;
Array<int> dof_edges, dof_boundary_elements, ess_edge_list;
fespace.GetBoundaryLoopEdgeDofs(boundary_elements, ess_tdof_list,
boundary_edge_ldofs, &ldof_marker, &dof_edges,
&dof_boundary_elements, &ess_edge_list);
Array<int> dof_orientations;
fespace.ComputeLoopEdgeOrientations(dof_edges, dof_boundary_elements,
test.normal, dof_orientations);
ParGridFunction x(&fespace);
x = real_t(0.0);
for (int i = 0; i < boundary_edge_ldofs.Size(); i++)
{
x(boundary_edge_ldofs[i]) = real_t(1.0) * dof_orientations[i];
}
GroupCommunicator *gc = fespace.ScalarGroupComm();
Array<int> global_marker(ldof_marker);
gc->Reduce<int>(global_marker.GetData(), GroupCommunicator::BitOR<int>);
gc->Bcast(global_marker);
Array<real_t> values(x.GetData(), x.Size());
gc->ReduceBegin(values.GetData());
gc->ReduceMarked<real_t>(values.GetData(), global_marker, 0,
GroupCommunicator::MaxAbs<real_t>);
gc->Bcast(values.GetData());
delete gc;
Vector x_true;
x.GetTrueDofs(x_true);
int local_nonzero_tdofs = 0;
real_t local_tdof_sum = 0.0;
for (int tdof = 0; tdof < x_true.Size(); tdof++)
{
real_t tdof_value = x_true(tdof);
if (abs(tdof_value) > 1e-12)
{
local_nonzero_tdofs++;
local_tdof_sum += abs(tdof_value);
}
}
real_t local_length = ComputeBoundaryLoopLength(&pmesh, dof_edges);
int global_nonzero_tdofs;
real_t global_tdof_sum, total_length;
MPI_Allreduce(&local_nonzero_tdofs, &global_nonzero_tdofs, 1, MPI_INT, MPI_SUM,
MPI_COMM_WORLD);
MPI_Allreduce(&local_tdof_sum, &global_tdof_sum, 1,
MPITypeMap<real_t>::mpi_type, MPI_SUM,
MPI_COMM_WORLD);
MPI_Allreduce(&local_length, &total_length, 1,
MPITypeMap<real_t>::mpi_type, MPI_SUM,
MPI_COMM_WORLD);
// Verify processor-invariant results match expected values
REQUIRE(global_nonzero_tdofs == exp_tdofs);
REQUIRE(abs(global_tdof_sum - exp_sum) < real_t(1e-12));
REQUIRE(abs(total_length - exp_length) < real_t(1e-12));
}
}
TEST_CASE("BoundaryEdgeDOFs2DSquareInSquare",
"[Parallel][ParMesh][BoundaryEdgeDOFs]")
{
// Test 2D boundary edge DOF extraction using square-in-square mesh
constexpr int order = 2;
// Test multiple inner boundary attributes
std::vector<int> inner_attrs_to_test = {5, 6, 7, 8};
// Load 2D square-in-square mesh from file
const char* mesh_file = "../../data/square_in_square.msh";
Mesh serial_mesh(mesh_file, 1, 1);
serial_mesh.UniformRefinement();
int num_procs = Mpi::WorldSize();
// Test each boundary attribute
for (int inner_attr : inner_attrs_to_test)
{
CAPTURE(inner_attr); // Capture the attribute being tested for better test output
// Test that results are consistent across different mesh partitionings
const int n_elements = serial_mesh.GetNE();
// Generate multiple different partitionings
std::vector<std::vector<int>> all_partitionings;
// 1. All elements on rank 0
all_partitionings.push_back(std::vector<int>(n_elements, 0));
if (num_procs > 1)
{
// 2. Block partition: first half on rank 0, second half on last rank
std::vector<int> block(n_elements);
for (int i = 0; i < n_elements; i++)
{
block[i] = (i < n_elements/2) ? 0 : num_procs-1;
}
all_partitionings.push_back(block);
// 3. Round-robin partition: elements assigned cyclically to all ranks
std::vector<int> round_robin(n_elements);
for (int i = 0; i < n_elements; i++)
{
round_robin[i] = i % num_procs;
}
all_partitionings.push_back(round_robin);
}
ND_FECollection fec(order, 2);
std::vector<int> all_dof_results;
all_dof_results.reserve(all_partitionings.size());
// Test each partitioning
for (const auto& partition : all_partitionings)
{
// Create parallel mesh with current partitioning
Mesh test_mesh(mesh_file, 1, 1);
test_mesh.UniformRefinement();
ParMesh pmesh = (num_procs == 1) ?
ParMesh(MPI_COMM_WORLD, test_mesh) :
ParMesh(MPI_COMM_WORLD, test_mesh, partition.data());
ParFiniteElementSpace fespace(&pmesh, &fec);
// Find boundary elements with the inner attribute
std::vector<Array<int>> attr_to_elements;
Array<int> inner_attrs(1);
inner_attrs[0] = inner_attr;
fespace.GetBoundaryElementsByAttribute(inner_attrs, attr_to_elements);
Array<int> inner_boundary_elements = attr_to_elements[0];
Array<int> ess_tdofs, ess_edges;
Array<int> boundary_dofs;
Array<int> dof_edges, dof_boundary_elements;
fespace.GetBoundaryLoopEdgeDofs(inner_boundary_elements, ess_tdofs,
boundary_dofs, nullptr, &dof_edges,
&dof_boundary_elements, &ess_edges);
// The output arrays share one indexing, so their sizes must match.
REQUIRE(boundary_dofs.Size() == dof_edges.Size());
REQUIRE(dof_edges.Size() == dof_boundary_elements.Size());
// Gather global counts for this partitioning
int local_dof_count = boundary_dofs.Size();
int global_dof_count;
MPI_Allreduce(&local_dof_count, &global_dof_count, 1, MPI_INT, MPI_SUM,
MPI_COMM_WORLD);
all_dof_results.push_back(global_dof_count);
}
// Verify all partitionings give identical results
REQUIRE(!all_dof_results.empty());
int expected_dofs = all_dof_results[0];
for (int result : all_dof_results)
{
REQUIRE(result == expected_dofs);
}
} // End of inner_attr loop
}
TEST_CASE("BoundaryEdgeDOFsSharedDOFsAreOwnedBySomeRank",
"[Parallel][ParMesh][BoundaryEdgeDOFs]")
{
// Every selected shared DOF must appear in exactly one rank's ess_tdof_list.
// Only the group master owns the corresponding true DOF and returns a
// non-negative value from GetLocalTDofNumber(), so if the master holds none
// of the selected boundary elements the DOF would be emitted by no rank at
// all unless the local marker is synchronized across the sharing group.
const int nranks = Mpi::WorldSize();
if (nranks < 2) { return; }
constexpr int order = 1;
ND_FECollection fec(order, 3);
for (int orientation : {1, 3, 5})
{
Mesh probe = OrientedTriFaceMesh(orientation, true);
probe.UniformRefinement();
const int ne = probe.GetNE();
// Several partitionings, to vary which rank masters each shared group
std::vector<std::vector<int>> partitionings;
{
std::vector<int> round_robin(ne), block(ne), strided(ne);
for (int i = 0; i < ne; i++)
{
round_robin[i] = i % nranks;
block[i] = (i < ne/2) ? 0 : nranks-1;
strided[i] = (i * 7 + 3) % nranks;
}
partitionings = {round_robin, block, strided};
}
for (const auto &partition : partitionings)
{
Mesh mesh = OrientedTriFaceMesh(orientation, true);
mesh.UniformRefinement();
ParMesh pmesh(MPI_COMM_WORLD, mesh, partition.data());
ParFiniteElementSpace fes(&pmesh, &fec);
const int bdr_attr = pmesh.bdr_attributes.Max();
Array<int> bdr_attrs(1);
bdr_attrs[0] = bdr_attr;
std::vector<Array<int>> attr_to_elements;
fes.GetBoundaryElementsByAttribute(bdr_attrs, attr_to_elements);
Array<int> bdr_elements = attr_to_elements[0];
Array<int> ess_tdofs;
Array<int> boundary_dofs;
fes.GetBoundaryLoopEdgeDofs(bdr_elements, ess_tdofs, boundary_dofs);
// Identify DOFs by global true DOF number, which is agreed upon by all
// ranks sharing the DOF, then compare the set selected anywhere with
// the set actually emitted in ess_tdof_list.
std::set<HYPRE_BigInt> selected, emitted;
for (int dof : boundary_dofs)
{
selected.insert(fes.GetGlobalTDofNumber(dof));
}
for (int i = 0; i < ess_tdofs.Size(); i++)
{
emitted.insert(fes.GetMyTDofOffset() + ess_tdofs[i]);
}
auto all_gather = [nranks](const std::set<HYPRE_BigInt> &s)
{
std::vector<HYPRE_BigInt> local(s.begin(), s.end());
int n = static_cast<int>(local.size()), total = 0;
std::vector<int> counts(nranks), bytes(nranks), displs(nranks);
MPI_Allgather(&n, 1, MPI_INT, counts.data(), 1, MPI_INT,
MPI_COMM_WORLD);
constexpr int sz = sizeof(HYPRE_BigInt);
for (int r = 0; r < nranks; r++)
{
displs[r] = total * sz;
total += counts[r];
bytes[r] = counts[r] * sz;
}
std::vector<HYPRE_BigInt> all(total);
MPI_Allgatherv(local.data(), n * sz, MPI_BYTE, all.data(),
bytes.data(), displs.data(), MPI_BYTE,
MPI_COMM_WORLD);
return std::set<HYPRE_BigInt>(all.begin(), all.end());
};
// Gather both sets across all ranks. global_selected is every shared
// boundary DOF chosen on any rank; global_emitted is every true DOF
// actually placed in some rank's ess_tdof_list. A selected DOF missing
// from global_emitted is one that no rank owns and outputs, which is
// exactly the synchronization bug this test guards against.
const std::set<HYPRE_BigInt> global_selected = all_gather(selected);
const std::set<HYPRE_BigInt> global_emitted = all_gather(emitted);
int num_missing = 0;
for (auto gtdof : global_selected)
{
if (!global_emitted.count(gtdof)) { num_missing++; }
}
CAPTURE(orientation, nranks, global_selected.size(),
global_emitted.size(), num_missing);
REQUIRE(num_missing == 0);
}
}
}
TEST_CASE("BoundaryEdgeDOFs2DLoopVertexDOFsPartitionInvariant",
"[Parallel][ParMesh][BoundaryEdgeDOFs]")
{
// A closed boundary loop split between ranks must give the same result as
// the serial code. With a collection carrying vertex DOFs (ND_R2D), a vertex
// shared by two boundary segments is interior to the loop and must be
// dropped. When the two segments live on different ranks, each rank sees the
// vertex only once locally, so the occurrence parity has to be reconciled
// across the sharing group.
if (Mpi::WorldSize() < 2) { return; }
constexpr int order = 1;
ND_R2D_FECollection fec(order, 2);
// Serial reference result
Mesh serial_mesh = Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL, false,
1.0, 1.0);
FiniteElementSpace serial_fes(&serial_mesh, &fec);
Array<int> serial_bdr_elements(serial_mesh.GetNBE());
for (int i = 0; i < serial_bdr_elements.Size(); i++)
{
serial_bdr_elements[i] = i;
}
Array<int> serial_boundary_dofs;
serial_fes.GetBoundaryLoopEdgeDofs(serial_bdr_elements, serial_boundary_dofs);
const int serial_count = serial_boundary_dofs.Size();
// Compare against several partitionings of the same mesh
const int num_procs = Mpi::WorldSize();
std::vector<std::vector<int>> partitionings;
{
Mesh probe = Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL, false,
1.0, 1.0);
const int ne = probe.GetNE();
std::vector<int> block(ne), round_robin(ne);
for (int i = 0; i < ne; i++)
{
block[i] = (i < ne/2) ? 0 : num_procs-1;
round_robin[i] = i % num_procs;
}
partitionings.push_back(block);
partitionings.push_back(round_robin);
}
for (const auto &partition : partitionings)
{
Mesh mesh = Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL, false,
1.0, 1.0);
ParMesh pmesh(MPI_COMM_WORLD, mesh, partition.data());
ParFiniteElementSpace pfes(&pmesh, &fec);
Array<int> local_bdr_elements(pmesh.GetNBE());
for (int i = 0; i < local_bdr_elements.Size(); i++)
{
local_bdr_elements[i] = i;
}
Array<int> ess_tdofs;
Array<int> local_boundary_dofs;
pfes.GetBoundaryLoopEdgeDofs(local_bdr_elements, ess_tdofs,
local_boundary_dofs);
// The true DOFs are owned by exactly one rank each, so summing the local
// counts gives a partition-independent global count.
int local_tdofs = ess_tdofs.Size();
int global_tdofs = 0;
MPI_Allreduce(&local_tdofs, &global_tdofs, 1, MPI_INT, MPI_SUM,
MPI_COMM_WORLD);
CAPTURE(num_procs, serial_count, global_tdofs);
REQUIRE(global_tdofs == serial_count);
}
}
TEST_CASE("GroupCommunicatorReduceMarkedByGroupStride",
"[Parallel][GroupCommunicator]")
{
// Regression test for the neighbor-major stride of the byGroup receive
// buffer: with more than one DOF in a group, the contributions to DOF i are
// at buf[j*nldofs + i], so reducing a single marked DOF must gather the
// strided values rather than reading a contiguous run.
const int rank = Mpi::WorldRank();
const int nranks = Mpi::WorldSize();
if (nranks < 3) { return; }
ListOfIntegerSets groups;
IntegerSet local_group(1);
local_group[0] = rank;
groups.Insert(local_group);
IntegerSet shared_group(nranks);
for (int r = 0; r < nranks; r++)
{
shared_group[r] = r;
}
groups.Insert(shared_group);
GroupTopology topology(MPI_COMM_WORLD);
topology.Create(groups, 4983);
GroupCommunicator comm(topology, GroupCommunicator::byGroup);
// Two DOFs in the same shared group, so the buffer stride is 2.
Array<int> ldof_group(2);
ldof_group = 1;
comm.Create(ldof_group);
Array<real_t> values(2);
values[0] = real_t(10.0) * rank + real_t(1.0);
values[1] = real_t(100.0) * rank + real_t(2.0);
Array<int> marker(2);
marker = 1;
comm.ReduceBegin(values.GetData());
comm.ReduceMarked<real_t>(values.GetData(), marker, 0,
GroupCommunicator::Sum<real_t>);
comm.Bcast(values);
const real_t rank_sum = real_t(nranks) * real_t(nranks - 1) / real_t(2.0);
REQUIRE(values[0] == MFEM_Approx(real_t(10.0) * rank_sum + real_t(nranks)));
REQUIRE(values[1] == MFEM_Approx(real_t(100.0) * rank_sum +
real_t(2.0) * real_t(nranks)));
}
TEST_CASE("GroupCommunicatorMaxAbs", "[Parallel][GroupCommunicator]")
{
const int rank = Mpi::WorldRank();
const int nranks = Mpi::WorldSize();
if (nranks < 2) { return; }
ListOfIntegerSets groups;
IntegerSet local_group(1);
local_group[0] = rank;
groups.Insert(local_group);
IntegerSet shared_group(nranks);
for (int r = 0; r < nranks; r++)
{
shared_group[r] = r;
}
groups.Insert(shared_group);
GroupTopology topology(MPI_COMM_WORLD);
topology.Create(groups, 4983);
GroupCommunicator comm(topology, GroupCommunicator::byGroup);
Array<int> ldof_group(2);
ldof_group = 1;
comm.Create(ldof_group);
// The group master (rank 0) reduces the peers' contributions into its own.
Array<real_t> values(2);
// DOF 0: equal magnitude across ranks with opposite signs, with the negative
// value held by the master, so the opposite-sign tie must still resolve
// deterministically to the positive value.
values[0] = (rank == 0) ? real_t(-5.0) : real_t(5.0);
// DOF 1: the largest magnitude is negative and held by a peer, so it must
// win over the master's smaller positive value and keep its sign.
values[1] = (rank == nranks - 1) ? real_t(-10.0) : real_t(5.0);
comm.Reduce<real_t>(values.GetData(), GroupCommunicator::MaxAbs<real_t>);
comm.Bcast(values);
REQUIRE(values[0] == MFEM_Approx(real_t(5.0)));
REQUIRE(values[1] == MFEM_Approx(real_t(-10.0)));
}
#endif // MFEM_USE_MPI
+243
View File
@@ -0,0 +1,243 @@
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#include "mfem.hpp"
#include "unit_tests.hpp"
using namespace mfem;
namespace
{
void PackLower(const DenseMatrix &mat, real_t *packed)
{
const int n = mat.Height();
for (int j = 0; j < n; ++j)
{
for (int i = j; i < n; ++i)
{
packed[TriPackLowerMatrix::LowerIndex(i, j, n)] = mat(i, j);
}
}
}
void FillFullBatch(const DenseMatrix &mat, real_t *full)
{
const int n = mat.Height();
for (int j = 0; j < n; ++j)
{
for (int i = 0; i < n; ++i)
{
full[i + n*j] = mat(i, j);
}
}
}
void BuildLowerDense(const TriPackLowerMatrix &packed,
int e, DenseMatrix &mat)
{
const int n = packed.GetNumRows();
mat.SetSize(n);
mat = 0.0;
const real_t *data = packed.Data().HostRead() + e*packed.GetPackedSize();
for (int j = 0; j < n; ++j)
{
for (int i = j; i < n; ++i)
{
mat(i, j) = data[TriPackLowerMatrix::LowerIndex(i, j, n)];
}
}
}
}
TEST_CASE("TriPackLowerMatrix operations", "[TriPackLowerMatrix]")
{
constexpr int n = 3;
constexpr int batch_size = 2;
constexpr real_t tol = 1e-12;
DenseMatrix A0(n), A1(n);
A0 = 0.0;
A1 = 0.0;
A0(0,0) = 4.0; A0(0,1) = 1.0; A0(0,2) = 2.0;
A0(1,0) = 1.0; A0(1,1) = 5.0; A0(1,2) = 3.0;
A0(2,0) = 2.0; A0(2,1) = 3.0; A0(2,2) = 6.0;
A1(0,0) = 7.0; A1(0,1) = 2.0; A1(0,2) = 1.0;
A1(1,0) = 2.0; A1(1,1) = 8.0; A1(1,2) = 2.0;
A1(2,0) = 1.0; A1(2,1) = 2.0; A1(2,2) = 5.0;
TriPackLowerMatrix packed(n, batch_size);
packed = 0.0;
real_t *packed_data = packed.Data().HostWrite();
PackLower(A0, packed_data);
PackLower(A1, packed_data + packed.GetPackedSize());
Vector full(batch_size*n*n);
real_t *full_data = full.HostWrite();
FillFullBatch(A0, full_data);
FillFullBatch(A1, full_data + n*n);
SECTION("Compare with full symmetric matrices")
{
REQUIRE(tripack::CompareWithFull(packed, full, tol));
}
SECTION("Symmetric multiply and lumping")
{
Vector x({1.0, 2.0, 3.0, 4.0, 5.0, 6.0});
Vector y, lump;
tripack::Mult(packed, x, y);
tripack::Lump(packed, lump);
Vector y_expected(batch_size*n);
Vector lump_expected(batch_size*n);
y_expected = 0.0;
lump_expected = 0.0;
const DenseMatrix *mats[batch_size] = { &A0, &A1 };
for (int e = 0; e < batch_size; ++e)
{
const DenseMatrix &M = *mats[e];
for (int i = 0; i < n; ++i)
{
real_t rowsum = 0.0;
real_t val = 0.0;
for (int j = 0; j < n; ++j)
{
rowsum += M(i, j);
val += M(i, j) * x(e*n + j);
}
lump_expected(e*n + i) = rowsum;
y_expected(e*n + i) = val;
}
}
for (int i = 0; i < y.Size(); ++i)
{
REQUIRE(y(i) == MFEM_Approx(y_expected(i)));
REQUIRE(lump(i) == MFEM_Approx(lump_expected(i)));
}
}
SECTION("Lower Cholesky factor and solves")
{
TriPackLowerMatrix lfac;
Vector rhs({1.0, -1.0, 2.0, 0.5, 1.5, -2.0});
Vector y, t, x;
tripack::ComputeCholeskyLower(packed, lfac);
tripack::SolveLowerTranspose(lfac, rhs, t);
tripack::SolveLower(lfac, t, x);
tripack::SolveCholeskyLower(lfac, rhs, y);
const DenseMatrix *mats[batch_size] = { &A0, &A1 };
for (int e = 0; e < batch_size; ++e)
{
DenseMatrix L;
DenseMatrix recon(n);
BuildLowerDense(lfac, e, L);
MultABt(L, L, recon);
recon -= *mats[e];
REQUIRE(recon.MaxMaxNorm() == MFEM_Approx(0.0, tol, tol));
}
Vector x_expected(batch_size*n);
for (int e = 0; e < batch_size; ++e)
{
DenseMatrix inv(n);
CalcInverse(*mats[e], inv);
for (int i = 0; i < n; ++i)
{
real_t sum = 0.0;
for (int j = 0; j < n; ++j)
{
sum += inv(i, j) * rhs(e*n + j);
}
x_expected(e*n + i) = sum;
}
}
for (int i = 0; i < x.Size(); ++i)
{
REQUIRE(x(i) == MFEM_Approx(x_expected(i)).epsilon(tol));
REQUIRE(y(i) == MFEM_Approx(x_expected(i)).epsilon(tol));
}
}
SECTION("Cholesky lower inverse")
{
TriPackLowerMatrix linv;
Vector rhs({1.0, -1.0, 2.0, 0.5, 1.5, -2.0});
tripack::ComputeCholeskyLowerInverse(packed, linv);
const DenseMatrix *mats[batch_size] = { &A0, &A1 };
for (int e = 0; e < batch_size; ++e)
{
DenseMatrix L, recon(n), inv(n);
BuildLowerDense(linv, e, L);
MultAtB(L, L, recon);
CalcInverse(*mats[e], inv);
recon -= inv;
REQUIRE(recon.MaxMaxNorm() == MFEM_Approx(0.0, tol, tol));
}
Vector y(batch_size*n);
Vector y_expected(batch_size*n);
for (int e = 0; e < batch_size; ++e)
{
DenseMatrix Linv;
DenseMatrix inv(n);
CalcInverse(*mats[e], inv);
BuildLowerDense(linv, e, Linv);
Vector tmp(n);
tmp = 0.0;
for (int i = 0; i < n; ++i)
{
real_t sum = 0.0;
for (int j = 0; j <= i; ++j)
{
sum += Linv(i, j) * rhs(e*n + j);
}
tmp(i) = sum;
}
for (int i = 0; i < n; ++i)
{
real_t sum = 0.0;
for (int j = 0; j < n; ++j)
{
sum += inv(i, j) * rhs(e*n + j);
}
y_expected(e*n + i) = sum;
}
for (int i = 0; i < n; ++i)
{
real_t sum = 0.0;
for (int j = i; j < n; ++j)
{
sum += Linv(j, i) * tmp(j);
}
y(e*n + i) = sum;
}
}
for (int i = 0; i < y.Size(); ++i)
{
REQUIRE(y(i) == MFEM_Approx(y_expected(i)).epsilon(tol));
}
}
}
+166
View File
@@ -0,0 +1,166 @@
// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#include "mfem.hpp"
#include "unit_tests.hpp"
using namespace mfem;
#if defined(MFEM_USE_MAGMA) && (defined(MFEM_USE_HIP) || defined(MFEM_USE_CUDA))
#include <vector>
#if defined(MFEM_USE_HIP)
#include <hip/hip_runtime.h>
#elif defined(MFEM_USE_CUDA)
#include <cuda_runtime.h>
#endif
namespace
{
bool HasGpuDevice()
{
#if defined(MFEM_USE_HIP)
int count = 0;
const hipError_t err = hipGetDeviceCount(&count);
return err == hipSuccess && count > 0;
#elif defined(MFEM_USE_CUDA)
int count = 0;
const cudaError_t err = cudaGetDeviceCount(&count);
return err == cudaSuccess && count > 0;
#else
return false;
#endif
}
DenseMatrix MakeSPD(const int n, const int seed)
{
DenseMatrix B(n), A(n);
for (int j = 0; j < n; ++j)
{
for (int i = 0; i < n; ++i)
{
const int v = (17*(i + 1) + 31*(j + 1) + 7*seed) % 23;
B(i, j) = 0.05 * real_t(v);
}
}
MultAtB(B, B, A);
for (int i = 0; i < n; ++i) { A(i, i) += 1.0 + 0.1*i; }
return A;
}
void PackLower(const DenseMatrix &mat, real_t *packed)
{
const int n = mat.Height();
for (int j = 0; j < n; ++j)
{
for (int i = j; i < n; ++i)
{
packed[TriPackLowerMatrix::LowerIndex(i, j, n)] =
mat(i, j);
}
}
}
real_t MaxResidual(const DenseMatrix &A, const Vector &x, const Vector &b)
{
const int n = A.Height();
real_t max_abs = 0.0;
for (int i = 0; i < n; ++i)
{
real_t sum = 0.0;
for (int j = 0; j < n; ++j) { sum += A(i, j) * x(j); }
max_abs = std::max(max_abs, std::abs(sum - b(i)));
}
return max_abs;
}
} // namespace
TEST_CASE("MAGMA packed-lower Cholesky factor+solve",
"[MAGMA][TriPackLowerMatrix]")
{
if (!HasGpuDevice())
{
WARN("No GPU device visible; skipping MAGMA packed-lower tests.");
return;
}
Device device(
#if defined(MFEM_USE_HIP)
"hip"
#elif defined(MFEM_USE_CUDA)
"cuda"
#else
"cpu"
#endif
);
constexpr int n = 8;
constexpr int batch_size = 17;
constexpr double tol = 5e-9;
TriPackLowerMatrix A_packed(n, batch_size);
A_packed.UseDevice(true);
A_packed = 0.0;
std::vector<DenseMatrix> A_dense;
A_dense.reserve(batch_size);
real_t *h_packed = A_packed.Data().HostWrite();
const int ps = A_packed.GetPackedSize();
for (int e = 0; e < batch_size; ++e)
{
A_dense.emplace_back(MakeSPD(n, e + 1));
PackLower(A_dense.back(), h_packed + e*ps);
}
Vector b(batch_size*n);
Vector x(batch_size*n);
real_t *h_b = b.HostWrite();
for (int e = 0; e < batch_size; ++e)
{
for (int i = 0; i < n; ++i)
{
h_b[e*n + i] = 1.0 + real_t((13*(i + 1) + 7*(e + 1)) % 29)/real_t(29);
}
}
b.UseDevice(true);
x.UseDevice(true);
TriPackLowerMatrix L;
MagmaPackedLowerCholesky ws;
ws.Factor(A_packed, L);
x = b;
ws.SolveInPlace(L, x);
MFEM_DEVICE_SYNC;
const real_t *h_x = x.HostRead();
const real_t *h_b_ro = b.HostRead();
for (int e = 0; e < batch_size; ++e)
{
Vector xe(const_cast<real_t *>(h_x) + e*n, n);
Vector be(const_cast<real_t *>(h_b_ro) + e*n, n);
const real_t res = MaxResidual(A_dense[e], xe, be);
REQUIRE(res == MFEM_Approx(0.0, tol, tol));
}
}
#else
TEST_CASE("MAGMA packed-lower tests disabled", "[MAGMA][TriPackLowerMatrix]")
{
SUCCEED("MFEM was built without MAGMA+GPU support.");
}
#endif