PLaSK library

Loading...
Searching...
No Matches
diffusion2.py
Go to the documentation of this file.
1 import sys
2 import re
3
4 import numpy as np
5 import matplotlib.pyplot as plt
6
7 import sympy as sym
8 from sympy import latex
9 from sympy.printing import cxxcode
10
11 try:
12 from sympy.printing.cxx import CXX11CodePrinter
13 except ImportError:
14 from sympy.printing.cxxcode import CXX11CodePrinter
15
16 from mpi4py.MPI import COMM_WORLD as mpi
17
18msize = mpi.Get_size()
19mrank = mpi.Get_rank()
20
21
22 def evaluate(expr, args):
23 largs = np.array_split(args, msize)[mrank]
24 if len(largs.shape) == 1:
25 lres = [expr(arg) for arg in largs]
26 else:
27 lres = [expr(*arg) for arg in largs]
28 res = mpi.gather(lres, root=0)
29 if mrank == 0:
30 return np.concatenate(res)
31
32
33 def evaluate_matrix(expr):
34 args = [(i,j) for j in range(12) for i in range(j+1)]
35 vals = evaluate(expr, args)
36 if mrank == 0:
37 matrix = sym.Matrix(np.zeros((12, 12)))
38 for (i,j), val in zip(args, vals):
39 matrix[i,j] = val
40 if i != j:
41 matrix[j,i] = val
42 return matrix
43
44
45 def evaluate_vector(expr):
46 args = range(12)
47 vals = evaluate(expr, args)
48 if mrank == 0:
49 vector = sym.Matrix(np.zeros(12))
50 for i, val in zip(args, vals):
51 vector[i] = val
52 return vector
53
54
55 def _print_Indexed(self, expr):
56 indices = expr.indices
57 if len(indices) == 1:
58 # return f"{self._print(expr.base.label)}[idx(e,{self._print(indices[0])})]"
59 return f"{self._print(expr.base.label)}{self._print(indices[0])}"
60 else:
61 # return f"{self._print(expr.base.label)}[idx(e,{','.join(self._print(i) for i in indices)})]"
62 return f"{self._print(expr.base.label)}{''.join(self._print(i) for i in indices)}"
63CXX11CodePrinter._print_Indexed = _print_Indexed
64
65
66 def cpp(v):
67 s = cxxcode(sym.simplify(v), standard='C++11')
68 #for i in range(1, 5):
69 # s = s.replace(f"a[{i}]", f"a{i}")
70 s = re.sub(r'std::pow\‍(([^)]+), (\d+)\‍)', r'std::pow(\1,\2)', s)
71 s = re.sub(r'std::pow\‍(([^)]+),2\‍)', r'(\1*\1)', s)
72 s = re.sub(r'std::pow\‍(([^)]+),3\‍)', r'(\1*\1*\1)', s)
73 s = re.sub(r'\b(X|Y)\b', r'e.\1', s)
74 s = re.sub(r'\bU(\d\d)\b', r'U[e.i\1]', s)
75 s = re.sub(r'\bJ(\d\d)\b', r'J[e.n\1]', s)
76 s = re.sub(r'\bP(\d\d)(\d)\b', r'P[e.n\1].c\2\2', s)
77 s = re.sub(r'(d?G)(\d)\b', r'\1.c\2\2', s)
78 return s
79
80
81sidx = lambda i: "e.i{}{}".format(*idx[i])
82
83 def _print_cpp(kij, kvals, fvals, file=sys.stdout):
84 for (i,j), kval in zip(kij, kvals):
85 print(f"K({sidx(i)}, {sidx(j)}) += {kval};", file=file)
86 for i, fval in enumerate(fvals):
87 print(f"F[{sidx(i)}] += {fval};", file=file)
88
89
90 def print_cpp(KK, FF, fname=None):
91 kij = [(i,j) for j in range(12) for i in range(j+1)]
92 kvals = evaluate(lambda i,j: ' + '.join(cpp(K[i,j]) for K in KK), kij)
93 fvals = evaluate(lambda i: ' + '.join(cpp(F[i]) for F in FF), range(12))
94 if mrank == 0:
95 if fname is None:
96 _print_cpp(kij, kvals, fvals)
97 else:
98 with open(fname, 'w') as file:
99 _print_cpp(kij, kvals, fvals, file=file)
100
101
102# Podstawowe obliczenia
103
104A, B, C, D, = sym.symbols('A B C D')
105
106x, y = sym.symbols('x y')
107X, Y = sym.symbols('X Y')
108
109ξ = x / X
110ζ = y / Y
111
112Φ = np.linalg.inv(np.array([[1, 0, 0, 0], [0, 1, 0, 0], [1, 1, 1, 1], [0, 1, 2, 3]])).T
113
114φx = [sum(int(c) * ξ**i for i, c in enumerate(p)) for p in Φ]
115φx[1] *= X
116φx[3] *= X
117φx = np.array([sym.expand(a) for a in φx])
118
119φy = [sum(int(c) * ζ**i for i, c in enumerate(p)) for p in Φ]
120φy[1] *= Y
121φy[3] *= Y
122φy = np.array([sym.expand(a) for a in φy])
123
124idx = [(0,0), (0,1), (1,0), (0,2), (0,3), (1,2), (2,0), (2,1), (3,0), (2,2), (2,3), (3,2)]
125φ = np.array([φx[idx[i][0]] * φy[idx[i][1]] for i in range(12)])
126
127dφx = np.array([sym.expand(sym.diff(a, x)) for a in φ])
128dφy = np.array([sym.expand(sym.diff(a, y)) for a in φ])
129
130U = sym.IndexedBase('U', shape=(4, 4))
131u0 = sum(U[idx[k]] * φ[k] for k in range(12))
132
133 if mrank == 0:
134 print("U:")
135 print(" return ", cpp(u0).replace('U', 'active.U'), ';', sep='')
136
137zx = np.array([0, X])
138zy = np.array([0, Y])
139
140J = sym.IndexedBase('J', shape=(2, 2))
141
142
143 def La(i, x, z):
144 res = 1
145 for j in range(len(z)):
146 if j == i: continue
147 res *= (x - z[j]) / (z[i] - z[j])
148 return res
149
150
151Jxy = sym.simplify(sum(J[i,j] * La(i, x, zx) * La(j, y, zy) for i in range(len(zx)) for j in range(len(zy))))
152
153
154KD = evaluate_matrix(lambda i,j:
155 sym.simplify(sym.integrate(sym.integrate(D * (dφx[i] * dφx[j] + dφy[i] * dφy[j]), (x, 0, X)), (y, 0, Y)))
156)
157 if mrank == 0: print('KD = sym.', KD, end='\n\n')
158
159KA = evaluate_matrix(lambda i,j:
160 sym.simplify(sym.integrate(sym.integrate(A * φ[i] * φ[j], (x, 0, X)), (y, 0, Y)))
161)
162 if mrank == 0: print('KA = sym.', KA, end='\n\n')
163
164KB = evaluate_matrix(lambda i,j:
165 sym.simplify(sym.integrate(sym.integrate(2 * B * u0 * φ[i] * φ[j], (x, 0, X)), (y, 0, Y)))
166)
167 if mrank == 0: print('KB = sym.', KB, end='\n\n')
168
169u02 = sym.expand(u0**2)
170u03 = sym.expand(u0**3)
171
172KC = evaluate_matrix(lambda i,j:
173 sym.simplify(sym.integrate(sym.simplify(sym.integrate(3 * C * u02 * φ[i] * φ[j], (x, 0, X))), (y, 0, Y)))
174)
175 if mrank == 0: print('KC = sym.', KC, end='\n\n')
176
177 if mrank == 0:
178 KK = KD, KA, KB, KC
179 else:
180 KK = None
181KK = mpi.bcast(KK, root=0)
182
183FB = evaluate_vector(lambda i:
184 sym.simplify(sym.integrate(sym.simplify(sym.integrate(B * u02 * φ[i], (x, 0, X))), (y, 0, Y)))
185)
186 if mrank == 0: print('FB = sym.', FB, end='\n\n')
187
188FC = evaluate_vector(lambda i:
189 sym.simplify(sym.integrate(sym.simplify(sym.integrate(2 * C * u03 * φ[i], (x, 0, X))), (y, 0, Y)))
190)
191 if mrank == 0: print('FC = sym.', FC, end='\n\n')
192
193F0 = evaluate_vector(lambda i:
194 sym.simplify(sym.integrate(sym.integrate(Jxy * φ[i], (x, 0, X)), (y, 0, Y)))
195)
196 if mrank == 0: print('F0 = sym.', F0, end='\n\n')
197
198 if mrank == 0:
199 FF = FB, FC, F0
200 else:
201 FF = None
202FF = mpi.bcast(FF, root=0)
203
204 print_cpp(KK, FF, 'diffusion3d-eval.ipp')
205
206
207# Wypalanie nośników
208
209P = sym.IndexedBase('P', shape=(2,2,2))
210G = sym.IndexedBase('G', shape=(2,))
211dG = sym.IndexedBase('dG', shape=(2,))
212ug = sym.symbols('Ug')
213
214Pxy = [sym.simplify(sum(P[i, j, c] * La(i, x, zx) * La(j, y, zy) for i in range(len(zx)) for j in range(len(zy)))) for c in range(2)]
215
216KL = evaluate_matrix(lambda i,j:
217 sym.simplify(sym.integrate(sym.integrate(sum(Pxy[i] * dG[i] for i in range(2)) * φ[i] * φ[j], (x, 0, X)), (y, 0, Y)))
218)
219 if mrank == 0: print('KL = sym.', KL, end='\n\n')
220
221FL = evaluate_vector(lambda i:
222 sym.simplify(sym.integrate(sym.integrate(sum(Pxy[i] * (dG[i] * ug - G[i]) for i in range(2)) * φ[i], (x, 0, X)), (y, 0, Y)))
223)
224 if mrank == 0: print('FL = sym.', FL, end='\n\n')
225
226KL = mpi.bcast(KL, root=0)
227FL = mpi.bcast(FL, root=0)
228
229 print_cpp((KL,), (FL,), 'diffusion3d-eval-shb.ipp')