Skip to content

Commit ad6ac76

Browse files
felsenhowerfelsenhower
committed
Remove perturbation_matrix; use same logic as in reference
implementation
1 parent d4f17a2 commit ad6ac76

6 files changed

Lines changed: 87 additions & 33 deletions

File tree

‎cython/calculate.pyx‎

Lines changed: 23 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77

88
from time import time
99

10-
from libc.math cimport fabs
10+
from libc.math cimport fabs, sin
1111

1212
import numpy as np
1313
cimport numpy as np
@@ -16,6 +16,7 @@ from partdiff_common.parse_args import (
1616
Options,
1717
CalculationMethod,
1818
TerminationCondition,
19+
PerturbationFunction,
1920
)
2021

2122
from partdiff_common import (
@@ -27,37 +28,51 @@ cdef int METHOD_GAUSS_SEIDEL = CalculationMethod.GAUSS_SEIDEL.value
2728
cdef int METHOD_JACOBI = CalculationMethod.JACOBI.value
2829
cdef int TERMINATION_ACCURACY = TerminationCondition.ACCURACY.value
2930
cdef int TERMINATION_ITERATIONS = TerminationCondition.ITERATIONS.value
31+
cdef int PERT_FUNC_F0 = PerturbationFunction.F0.value
32+
cdef int PERT_FUNC_FPISIN = PerturbationFunction.FPISIN.value
33+
34+
cdef double PI = 3.14159265358979323846
3035

3136
cdef tuple calculate_inner(
3237
int method,
38+
int pert_func,
3339
int termination,
3440
int term_iteration,
3541
double term_accuracy,
3642
int n,
43+
double h,
3744
double[:, :, :] tensor,
38-
double[:, :] perturbation_matrix,
3945
):
4046
cdef int i, j
41-
cdef double star, residuum, maxresiduum
47+
cdef double star, residuum, maxresiduum, pih, fpisin, fpisin_i
4248
cdef int stat_iteration = 0
4349
cdef double stat_accuracy = 0.0
4450
cdef double[:, :] matrix_out = tensor[0, :, :]
4551
cdef double[:, :] matrix_in = matrix_out
4652
cdef double[:, :] final_matrix
53+
pih = 0.0
54+
fpisin = 0.0
55+
fpisin_i = 0.0
4756
if method == METHOD_JACOBI:
4857
matrix_in = tensor[1, :, :]
58+
if pert_func == PERT_FUNC_FPISIN:
59+
pih = PI * h
60+
fpisin = 0.25 * (2.0 * PI * PI) * h * h
4961
while True:
5062
stat_iteration += 1
5163
maxresiduum = 0.0
5264
for i in range(1, n):
65+
if pert_func == PERT_FUNC_FPISIN:
66+
fpisin_i = fpisin * sin(pih * i)
5367
for j in range(1, n):
5468
star = 0.25 * (
5569
matrix_in[i - 1, j] +
5670
matrix_in[i, j - 1] +
5771
matrix_in[i, j + 1] +
5872
matrix_in[i + 1, j]
5973
)
60-
star += perturbation_matrix[i, j]
74+
if pert_func == PERT_FUNC_FPISIN:
75+
star += fpisin_i * sin(pih * j)
6176
if (
6277
termination == TERMINATION_ACCURACY
6378
or term_iteration == stat_iteration
@@ -79,21 +94,23 @@ cdef tuple calculate_inner(
7994

8095
def calculate(arguments: CalculationArguments, options: Options) -> CalculationResults:
8196
cdef int method = options.method.value
97+
cdef int pert_func = options.pert_func.value
8298
cdef int termination = options.termination.value
8399
cdef int term_iteration = options.term_iteration
84100
cdef double term_accuracy = options.term_accuracy
85101
cdef int n = arguments.n
102+
cdef double h = arguments.h
86103
cdef double[:, :, :] tensor = arguments.tensor
87-
cdef double[:, :] perturbation_matrix = arguments.perturbation_matrix
88104
start_time = time()
89105
final_matrix, stat_iteration, stat_accuracy = calculate_inner(
90106
method,
107+
pert_func,
91108
termination,
92109
term_iteration,
93110
term_accuracy,
94111
n,
112+
h,
95113
tensor,
96-
perturbation_matrix,
97114
)
98115
end_time = time()
99116
duration = end_time - start_time

‎np_vectorize/main.py‎

Lines changed: 21 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
from time import time
88

99
import numpy as np
10+
from numpy import sin
1011
from partdiff_common import (
1112
CalculationArguments,
1213
CalculationResults,
@@ -18,10 +19,13 @@
1819
from partdiff_common.parse_args import (
1920
CalculationMethod,
2021
Options,
22+
PerturbationFunction,
2123
TerminationCondition,
2224
parse_args,
2325
)
2426

27+
PI = 3.14159265358979323846
28+
2529

2630
def calculate_jacobi(
2731
arguments: CalculationArguments, options: Options
@@ -37,11 +41,19 @@ def calculate_jacobi(
3741
"""
3842
start_time = time()
3943
n = arguments.n
44+
h = arguments.h
4045
tensor = arguments.tensor
41-
perturbation_matrix = arguments.perturbation_matrix
4246
stat_accuracy = None
4347
matrix_out = tensor[0, :, :]
4448
matrix_in = tensor[1, :, :]
49+
perturbation_matrix = np.zeros((n + 1, n + 1), dtype=np.float64)
50+
if options.pert_func == PerturbationFunction.FPISIN:
51+
pih = PI * h
52+
fpisin = 0.25 * (2.0 * PI * PI) * h * h
53+
for i in range(1, n):
54+
fpisin_i = fpisin * sin(pih * i)
55+
for j in range(1, n):
56+
perturbation_matrix[i, j] = fpisin_i * sin(pih * j)
4557
for stat_iteration in count(start=1):
4658
maxresiduum = 0.0
4759
center = matrix_in[1:n, 1:n]
@@ -85,22 +97,28 @@ def calculate_gauss_seidel(
8597
"""
8698
start_time = time()
8799
n = arguments.n
100+
h = arguments.h
88101
tensor = arguments.tensor
89-
perturbation_matrix = arguments.perturbation_matrix
90102
stat_iteration = 0
91103
stat_accuracy = None
92104
matrix = tensor[0, :, :]
105+
if options.pert_func == PerturbationFunction.FPISIN:
106+
pih = PI * h
107+
fpisin = 0.25 * (2.0 * PI * PI) * h * h
93108
for stat_iteration in count(start=1):
94109
maxresiduum = 0.0
95110
for i in range(1, n):
111+
if options.pert_func == PerturbationFunction.FPISIN:
112+
fpisin_i = fpisin * sin(pih * float(i))
96113
for j in range(1, n):
97114
star = 0.25 * (
98115
matrix[i - 1, j]
99116
+ matrix[i, j - 1]
100117
+ matrix[i, j + 1]
101118
+ matrix[i + 1, j]
102119
)
103-
star += perturbation_matrix[i, j]
120+
if options.pert_func == PerturbationFunction.FPISIN:
121+
star += fpisin_i * sin(pih * float(j))
104122
if (
105123
options.termination == TerminationCondition.ACCURACY
106124
or stat_iteration == options.term_iteration

‎nuitka/main.py‎

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
from itertools import count
77
from time import time
88

9+
from numpy import sin
910
from partdiff_common import (
1011
CalculationArguments,
1112
CalculationResults,
@@ -17,10 +18,13 @@
1718
from partdiff_common.parse_args import (
1819
CalculationMethod,
1920
Options,
21+
PerturbationFunction,
2022
TerminationCondition,
2123
parse_args,
2224
)
2325

26+
PI = 3.14159265358979323846
27+
2428

2529
def calculate(arguments: CalculationArguments, options: Options) -> CalculationResults:
2630
"""Solve the Poisson equation iteratively using the Jacobi or Gauß-Seidel method.
@@ -34,25 +38,31 @@ def calculate(arguments: CalculationArguments, options: Options) -> CalculationR
3438
"""
3539
start_time = time()
3640
n = arguments.n
41+
h = arguments.h
3742
tensor = arguments.tensor
38-
perturbation_matrix = arguments.perturbation_matrix
3943
stat_iteration = 0
4044
stat_accuracy = None
4145
matrix_out = tensor[0, :, :]
4246
matrix_in = matrix_out
4347
if options.method == CalculationMethod.JACOBI:
4448
matrix_in = tensor[1, :, :]
49+
if options.pert_func == PerturbationFunction.FPISIN:
50+
pih = PI * h
51+
fpisin = 0.25 * (2.0 * PI * PI) * h * h
4552
for stat_iteration in count(start=1):
4653
maxresiduum = 0.0
4754
for i in range(1, n):
55+
if options.pert_func == PerturbationFunction.FPISIN:
56+
fpisin_i = fpisin * sin(pih * float(i))
4857
for j in range(1, n):
4958
star = 0.25 * (
5059
matrix_in[i - 1, j]
5160
+ matrix_in[i, j - 1]
5261
+ matrix_in[i, j + 1]
5362
+ matrix_in[i + 1, j]
5463
)
55-
star += perturbation_matrix[i, j]
64+
if options.pert_func == PerturbationFunction.FPISIN:
65+
star += fpisin_i * sin(pih * float(j))
5666
if (
5767
options.termination == TerminationCondition.ACCURACY
5868
or stat_iteration == options.term_iteration

‎numba/main.py‎

Lines changed: 16 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77

88
import numpy as np
99
from numba import njit
10+
from numpy import sin
1011
from partdiff_common import (
1112
CalculationArguments,
1213
CalculationResults,
@@ -18,34 +19,38 @@
1819
from partdiff_common.parse_args import (
1920
CalculationMethod,
2021
Options,
22+
PerturbationFunction,
2123
TermAccuracy,
2224
TerminationCondition,
2325
TermIterations,
2426
parse_args,
2527
)
2628

29+
PI = 3.14159265358979323846
30+
2731

2832
@njit
2933
def calculate_iterate(
3034
method: CalculationMethod,
35+
pert_func: PerturbationFunction,
3136
termination: TerminationCondition,
3237
term_iteration: TermIterations,
3338
term_accuracy: TermAccuracy,
3439
n: int,
40+
h: float,
3541
tensor: np.ndarray,
36-
perturbation_matrix: np.ndarray,
3742
) -> tuple[np.ndarray, int, float]:
3843
"""The inner calculation part of the calculate method which is just-in-time compiled
3944
with numba here.
4045
4146
Args:
4247
method (CalculationMethod): The method (Gauß-Seidel or Jacobi).
48+
pert_func (PerturbationFunction): The perturbation function.
4349
termination (TerminationCondition): Termination (Iterations or Accuracy).
4450
term_iteration (TermIterations): Max iterations.
4551
term_accuracy (TermAccuracy): Min accuracy.
4652
n (int): Problem size (matrix is (n+1)*(n+1)).
4753
tensor (np.ndarray): The problem matrices.
48-
perturbation_matrix (np.ndarray): Precomputed perturbation function values.
4954
5055
Returns:
5156
tuple[np.ndarray, int, float]: A tuple containing the final matrix,
@@ -57,18 +62,24 @@ def calculate_iterate(
5762
matrix_in = matrix_out
5863
if method == CalculationMethod.JACOBI:
5964
matrix_in = tensor[1, :, :]
65+
if pert_func == PerturbationFunction.FPISIN:
66+
pih = PI * h
67+
fpisin = 0.25 * (2.0 * PI * PI) * h * h
6068
while True:
6169
stat_iteration += 1
6270
maxresiduum = 0.0
6371
for i in range(1, n):
72+
if pert_func == PerturbationFunction.FPISIN:
73+
fpisin_i = fpisin * sin(pih * float(i))
6474
for j in range(1, n):
6575
star = 0.25 * (
6676
matrix_in[i - 1, j]
6777
+ matrix_in[i, j - 1]
6878
+ matrix_in[i, j + 1]
6979
+ matrix_in[i + 1, j]
7080
)
71-
star += perturbation_matrix[i, j]
81+
if pert_func == PerturbationFunction.FPISIN:
82+
star += fpisin_i * sin(pih * float(j))
7283
if (
7384
termination == TerminationCondition.ACCURACY
7485
or stat_iteration == term_iteration
@@ -101,12 +112,13 @@ def calculate(arguments: CalculationArguments, options: Options) -> CalculationR
101112
start_time = time()
102113
final_matrix, stat_iteration, stat_accuracy = calculate_iterate(
103114
options.method,
115+
options.pert_func,
104116
options.termination,
105117
options.term_iteration,
106118
options.term_accuracy,
107119
arguments.n,
120+
arguments.h,
108121
arguments.tensor,
109-
arguments.perturbation_matrix,
110122
)
111123
end_time = time()
112124
duration = end_time - start_time

‎partdiff_common/src/partdiff_common/__init__.py‎

Lines changed: 3 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,6 @@
55

66
import sys
77
from dataclasses import dataclass
8-
from math import sin
98

109
import numpy as np
1110
from pympler import asizeof
@@ -22,12 +21,10 @@ class CalculationArguments:
2221
"""This class contains the internal representation of the problem, i.e. the
2322
initialized matrices. All matrices have the size (n+1)*(n+1).
2423
The tensor may be 1*(n+1)*(n+1) or 2*(n+1)*(n+1), depending on the method.
25-
The perturbation matrix contains the precomputed values of the perturbation function.
2624
"""
27-
2825
n: int
26+
h: float
2927
tensor: np.ndarray
30-
perturbation_matrix: np.ndarray
3128

3229

3330
@dataclass(frozen=True)
@@ -52,8 +49,7 @@ def init_arguments(options: Options) -> CalculationArguments:
5249
n = (options.interlines * 8) + 9 - 1
5350
num_matrices = 2 if options.method == CalculationMethod.JACOBI else 1
5451
h = 1.0 / n
55-
matrix_shape = (n + 1, n + 1)
56-
tensor_shape = (num_matrices, *matrix_shape)
52+
tensor_shape = (num_matrices, n + 1, n + 1)
5753
tensor = np.zeros(tensor_shape, dtype=np.float64)
5854
if options.pert_func == PerturbationFunction.F0:
5955
for g in range(num_matrices):
@@ -66,16 +62,7 @@ def init_arguments(options: Options) -> CalculationArguments:
6662
tensor[g, n, i] = c2
6763
tensor[g, n, 0] = 0.0
6864
tensor[g, 0, n] = 0.0
69-
perturbation_matrix = np.zeros(matrix_shape, dtype=np.float64)
70-
if options.pert_func == PerturbationFunction.FPISIN:
71-
pi = 3.14159265358979323846
72-
pih = pi * h
73-
fpisin = 0.25 * (2.0 * pi * pi) * h * h
74-
for i in range(1, n):
75-
fpisin_i = fpisin * sin(pih * i)
76-
for j in range(1, n):
77-
perturbation_matrix[i, j] = fpisin_i * sin(pih * j)
78-
return CalculationArguments(n, tensor, perturbation_matrix)
65+
return CalculationArguments(n, h, tensor)
7966

8067

8168
def calculate_memory_usage(

0 commit comments

Comments
 (0)