Note
Go to the end to download the full example code.
Symbolic SPH kernels & c++ tests#
7 from __future__ import division
8
9 import matplotlib.pyplot as plt
10 import mpmath
11 import numpy as np
12 from sympy import *
13
14 import shamrock
15
16 # Set precision to 32 digits
17 mpmath.mp.dps = 32
18
19 q = symbols("q")
Utilities#
26 def getnorm(w, R):
27 c1D = sympify(1) / (2 * integrate(w, (q, 0, R)))
28 c2D = sympify(1) / (integrate(2 * pi * q * w, (q, 0, R)))
29 c3D = sympify(1) / (integrate(4 * pi * q * q * w, (q, 0, R)))
30 return (c1D, c2D, c3D)
31
32
33 def replace_pi_constants(cpp_code):
34 """Replace M_PI and M_1_PI with shambase::constants::pi<Tscal>"""
35 cpp_code = cpp_code.replace("M_1_PI", "1/shambase::constants::pi<Tscal>")
36 cpp_code = cpp_code.replace("M_PI", "shambase::constants::pi<Tscal>")
37 return cpp_code
38
39
40 def sympy_to_cpp(expr):
41 """Convert a sympy expression to C++ code with proper formatting"""
42
43 import re
44
45 # Expand the expression first
46 # expr = expand(expr)
47 cpp_code = ccode(expr)
48
49 # Replace pow(q, n) with q*q*q... for powers up to 11
50 # for power in range(11, 1, -1): # From 11 down to 2
51 # q_mult = '*'.join(['q'] * power)
52 # cpp_code = cpp_code.replace(f'pow(q, {power})', q_mult)
53
54 # Replace pow(base, n) with sham::pow_constexpr<n>(base)
55 # Use regex to match pow(anything, number) pattern
56 def replace_pow_with_constexpr(match):
57 base = match.group(1)
58 exponent = match.group(2)
59 return f"sham::pow_constexpr<{exponent}>({base})"
60
61 cpp_code = re.sub(r"pow\(([^,]+),\s*(\d+)\)", replace_pow_with_constexpr, cpp_code)
62
63 # Wrap divisions in parentheses for clarity
64 # Pattern: number/number followed by * (but not already in parentheses)
65 # Match something like: 1155.0/512.0*q -> (1155.0/512.0)*q
66 cpp_code = re.sub(r"(\d+\.?\d*)/(\d+\.?\d*)(\*)", r"(\1/\2)\3", cpp_code)
67
68 # Replace pi constants
69 cpp_code = replace_pi_constants(cpp_code)
70
71 return cpp_code
72
73
74 def process_power_variables(cpp_exprs, indent=8):
75 """
76 Self-contained power variable processing:
77 1. Extracts sham::pow_constexpr calls from expressions
78 2. Replaces them with temporary variable names (t1, t2, ...)
79 3. Generates header declarations
80
81 Returns: (modified_exprs, power_header_string)
82 """
83 spaces = " " * indent
84
85 # Find all sham::pow_constexpr<n>(...) calls
86 # Need to handle nested parentheses properly
87 def find_pow_constexpr_calls(text):
88 """Find all sham::pow_constexpr calls with proper parenthesis matching"""
89 calls = []
90 pattern = "sham::pow_constexpr<"
91 pos = 0
92
93 while True:
94 pos = text.find(pattern, pos)
95 if pos == -1:
96 break
97
98 # Find the power number
99 power_start = pos + len(pattern)
100 power_end = text.find(">", power_start)
101 power = text[power_start:power_end]
102
103 # Find the matching closing parenthesis
104 paren_start = text.find("(", power_end)
105 if paren_start == -1:
106 break
107
108 # Count parentheses to find the matching closing one
109 paren_count = 1
110 i = paren_start + 1
111 while i < len(text) and paren_count > 0:
112 if text[i] == "(":
113 paren_count += 1
114 elif text[i] == ")":
115 paren_count -= 1
116 i += 1
117
118 if paren_count == 0:
119 # Found matching closing parenthesis
120 full_match = text[pos:i]
121 base = text[paren_start + 1 : i - 1]
122 calls.append((full_match, power, base))
123 pos = i
124 else:
125 pos += 1
126
127 return calls
128
129 # Dictionary to store unique pow_constexpr calls
130 # Key: full match string, Value: (var_name, power, base)
131 pow_calls = {}
132 var_counter = 1
133
134 vars = {}
135
136 for expr in cpp_exprs:
137 for full_match, power, base in find_pow_constexpr_calls(expr):
138 vars[full_match] = (power, base)
139
140 # print(vars)
141
142 # sort them lexicographically by match (from the end of the key)
143 sorted_vars = sorted(vars.items(), key=lambda x: x[0])
144 for full_match, (power, base) in sorted_vars:
145 if full_match not in pow_calls:
146 var_name = f"t{var_counter}"
147 pow_calls[full_match] = (var_name, power, base)
148 var_counter += 1
149
150 # Replace pow_constexpr calls with variable names in expressions
151 # Sort by length (longest first) to avoid partial replacements
152 sorted_pow_calls = sorted(pow_calls.items(), key=lambda x: len(x[0]), reverse=True)
153 for full_match, (var_name, power, base) in sorted_pow_calls:
154 cpp_exprs = [expr.replace(full_match, var_name) for expr in cpp_exprs]
155
156 # Generate header declarations
157 power_header = ""
158 for full_match, (var_name, power, base) in pow_calls.items():
159 power_header += f"{spaces} Tscal {var_name} = sham::pow_constexpr<{power}>({base});\n"
160
161 return cpp_exprs, power_header
162
163
164 def generate_function_body(pieces, cpp_exprs, indent=8):
165 """Generate the body section with if-else conditions and return statements"""
166 spaces = " " * indent
167 body = ""
168
169 for i, ((expr, cond), cpp_expr) in enumerate(zip(pieces, cpp_exprs)):
170 if cond == True:
171 body += f"{spaces} else\n"
172 body += f"{spaces} return {cpp_expr};\n"
173 else:
174 if_keyword = "if" if i == 0 else "else if"
175 # Convert condition to string
176 cond_str = ccode(cond)
177
178 body += f"{spaces} {if_keyword} ({cond_str}) {{\n"
179 body += f"{spaces} return {cpp_expr};\n"
180 body += f"{spaces} }}"
181 if i < len(pieces) - 1 and pieces[i + 1][1] != True:
182 body += " "
183 else:
184 body += "\n"
185
186 return body
187
188
189 def generate_piecewise_function(func_expr, func_name, indent=8):
190 """
191 Generate C++ code for a piecewise function.
192 Works with separate header (declarations) and body (logic) sections.
193 """
194 spaces = " " * indent
195
196 # Extract pieces from Piecewise object
197 pieces = []
198 for arg in func_expr.args:
199 # Each arg is an ExprCondPair (expr, cond) tuple
200 if hasattr(arg, "expr") and hasattr(arg, "cond"):
201 pieces.append((arg.expr, arg.cond))
202 elif isinstance(arg, tuple) and len(arg) == 2:
203 pieces.append(arg)
204
205 # Phase 1: Generate raw C++ expressions for the body
206 cpp_exprs = []
207 for expr, _ in pieces:
208 cpp_exprs.append(sympy_to_cpp(expr))
209
210 # Phase 2: Process power variables (self-contained: detect, replace, generate header)
211 cpp_exprs, power_header = process_power_variables(cpp_exprs, indent)
212
213 # Phase 3: Extract constants and replace in body
214 # constants = extract_constants(cpp_exprs)
215 # cpp_exprs = replace_constants_in_body(cpp_exprs, constants)
216
217 # Phase 4: Generate constant declarations header
218 # const_header = generate_constant_header(constants, indent)
219
220 # Phase 5: Generate body (if-else structure)
221 body = generate_function_body(pieces, cpp_exprs, indent)
222
223 # Phase 6: Combine everything (constants first, then powers, then body)
224 code = f"{spaces}inline static Tscal {func_name}(Tscal q) {{\n"
225 # if const_header:
226 # code += const_header
227 if power_header:
228 code += power_header
229 code += body
230 code += f"{spaces}}}\n"
231
232 return code
233
234
235 def kernel_to_shamrock(kernel_gen):
236 """
237 Generate complete C++ kernel definition from SymPy expression
238
239 Args:
240 f_func: Function that returns (R, f_expr, name) tuple
241 generate_phi: Whether to generate phi_tilde_3d (requires manual derivation)
242 """
243 R, f_expr, name = kernel_gen()
244
245 # Compute normalization constants
246 c_norm_1d, c_norm_2d, c_norm_3d = getnorm(f_expr, R)
247
248 class text_body:
249 def __init__(self):
250 self.text = ""
251
252 def __call__(self, text=""):
253 self.text += text + "\n"
254
255 text = text_body()
256
257 # Generate class header
258 text("template<class Tscal>")
259 text(f"class KernelDef{name} {{")
260 text(" public:")
261 text(
262 f" inline static constexpr Tscal Rkern = {ccode(R)}; ///< Compact support radius of the kernel"
263 )
264 text(" /// default hfact to be used for this kernel")
265 text(" inline static constexpr Tscal hfactd = 1.2;")
266 text()
267
268 # Normalize constants with pi handling
269 norm_1d_str = replace_pi_constants(ccode(c_norm_1d))
270 norm_2d_str = replace_pi_constants(ccode(c_norm_2d))
271 norm_3d_str = replace_pi_constants(ccode(c_norm_3d))
272
273 text(" /// 1D norm of the kernel")
274 text(f" inline static constexpr Tscal norm_1d = {norm_1d_str};")
275 text(" /// 2D norm of the kernel")
276 text(f" inline static constexpr Tscal norm_2d = {norm_2d_str};")
277 text(" /// 3D norm of the kernel")
278 text(f" inline static constexpr Tscal norm_3d = {norm_3d_str};")
279 text()
280
281 from sympy.polys.polyfuncs import horner
282
283 # Generate f function
284 text(generate_piecewise_function(f_expr, "f", indent=4))
285 text()
286
287 # Generate df function
288 df_expr = diff(f_expr, q)
289 text(generate_piecewise_function(df_expr, "df", indent=4))
290
291 # Generate ddf function
292 ddf_expr = diff(df_expr, q)
293 text(generate_piecewise_function(ddf_expr, "ddf", indent=4))
294
295 # phi_expr
296 phi_tilde_3d_expr = integrate(integrate(f_expr * 4 * pi * q * q, q) * (1 / q**2), q).simplify()
297
298 # correct phi to have 0 at infinity
299 filter_cond = next(expr for expr, cond in phi_tilde_3d_expr.args if cond == True)
300 lim_phi = limit(filter_cond, q, oo)
301 phi_tilde_3d_expr -= lim_phi
302 phi_tilde_3d_expr = phi_tilde_3d_expr.simplify()
303
304 text(generate_piecewise_function(phi_tilde_3d_expr, "phi_tilde_3d", indent=4))
305
306 # phi_tilde_3d_prime
307 phi_tilde_3d_prime_expr = diff(phi_tilde_3d_expr, q)
308
309 text(generate_piecewise_function(phi_tilde_3d_prime_expr, "phi_tilde_3d_prime", indent=4))
310
311 text("};")
312 text()
313
314 return {
315 "R": R,
316 "name": name,
317 "norm_1d": c_norm_1d.evalf(),
318 "norm_2d": c_norm_2d.evalf(),
319 "norm_3d": c_norm_3d.evalf(),
320 "Rkern": R,
321 "f": f_expr,
322 "df": df_expr,
323 "ddf": ddf_expr,
324 "phi_tilde_3d": phi_tilde_3d_expr,
325 "phi_tilde_3d_prime": phi_tilde_3d_prime_expr,
326 "text": text.text,
327 }
Testing the kernels#
335 def test_kernel(ret, tolerance=1e-12):
336 R = ret["R"]
337 name = ret["name"]
338 norm_1d = ret["norm_1d"]
339 norm_2d = ret["norm_2d"]
340 norm_3d = ret["norm_3d"]
341 Rkern = ret["Rkern"]
342 f_expr = ret["f"]
343 df_expr = ret["df"]
344 ddf_expr = ret["ddf"]
345 phi_tilde_3d_expr = ret["phi_tilde_3d"]
346 phi_tilde_3d_prime_expr = ret["phi_tilde_3d_prime"]
347 text = ret["text"]
348
349 print("------------------------------------------")
350 print(f"Testing kernel {name} matching with Shamrock code")
351 print("------------------------------------------")
352
353 f = lambdify((q), f_expr, modules=["mpmath"])
354 df = lambdify((q), df_expr, modules=["mpmath"])
355 ddf = lambdify((q), ddf_expr, modules=["mpmath"])
356 phi_tilde_3d = lambdify((q), phi_tilde_3d_expr, modules=["mpmath"])
357 phi_tilde_3d_prime = lambdify((q), phi_tilde_3d_prime_expr, modules=["mpmath"])
358
359 print("Testing norms:")
360 shamrock_norm_1d = getattr(shamrock.math.sphkernel, f"{name}_norm_1d")()
361 shamrock_norm_2d = getattr(shamrock.math.sphkernel, f"{name}_norm_2d")()
362 shamrock_norm_3d = getattr(shamrock.math.sphkernel, f"{name}_norm_3d")()
363
364 print(
365 f"Shamrock norm_1d = {shamrock_norm_1d}, expected={norm_1d}, delta={abs(shamrock_norm_1d - norm_1d)}"
366 )
367 print(
368 f"Shamrock norm_2d = {shamrock_norm_2d}, expected={norm_2d}, delta={abs(shamrock_norm_2d - norm_2d)}"
369 )
370 print(
371 f"Shamrock norm_3d = {shamrock_norm_3d}, expected={norm_3d}, delta={abs(shamrock_norm_3d - norm_3d)}"
372 )
373
374 # test the norm constants down to 1e-12
375 assert abs(shamrock_norm_1d - norm_1d) < tolerance
376 assert abs(shamrock_norm_2d - norm_2d) < tolerance
377 assert abs(shamrock_norm_3d - norm_3d) < tolerance
378
379 print()
380 print("Testing kernel radius:")
381 shamrock_Rkern = getattr(shamrock.math.sphkernel, f"{name}_Rkern")()
382
383 print(f"Shamrock Rkern = {shamrock_Rkern}, delta={abs(shamrock_Rkern - Rkern)}")
384 assert abs(shamrock_Rkern - Rkern) < tolerance
385
386 print()
387 print("Testing kernel functions:")
388 shamrock_f = getattr(shamrock.math.sphkernel, f"{name}_f")
389 shamrock_df = getattr(shamrock.math.sphkernel, f"{name}_df")
390 shamrock_ddf = getattr(shamrock.math.sphkernel, f"{name}_ddf")
391 shamrock_phi_tilde_3d = getattr(shamrock.math.sphkernel, f"{name}_phi_tilde_3d")
392 shamrock_phi_tilde_3d_prime = getattr(shamrock.math.sphkernel, f"{name}_phi_tilde_3d_prime")
393
394 print(f"Shamrock f(q) = {shamrock_f}")
395 print(f"Shamrock df(q) = {shamrock_df}")
396 print(f"Shamrock ddf(q) = {shamrock_ddf}")
397 print(f"Shamrock phi_tilde_3d(q) = {shamrock_phi_tilde_3d}")
398 print(f"Shamrock phi_tilde_3d_prime(q) = {shamrock_phi_tilde_3d_prime}")
399
400 q_arr = np.linspace(0, max(1.1 * float(Rkern), 5.0), 1000)
401 shamrock_f = [shamrock_f(x) for x in q_arr]
402 shamrock_df = [shamrock_df(x) for x in q_arr]
403 shamrock_ddf = [shamrock_ddf(x) for x in q_arr]
404 shamrock_phi_tilde_3d = [shamrock_phi_tilde_3d(x) for x in q_arr]
405 shamrock_phi_tilde_3d_prime = [shamrock_phi_tilde_3d_prime(x) for x in q_arr]
406
407 try:
408 sympy_f = [f(x) for x in q_arr]
409 sympy_df = [df(x) for x in q_arr]
410 sympy_ddf = [ddf(x) for x in q_arr]
411 sympy_phi_tilde_3d = [phi_tilde_3d(x) for x in q_arr]
412 sympy_phi_tilde_3d_prime = [phi_tilde_3d_prime(x) for x in q_arr]
413 except ZeroDivisionError:
414 # This can happen for some kernels at q=0
415 print("ZeroDivisionError: using limits towards 0")
416 q_arr_no_zero = q_arr[1:]
417
418 # Compute the limits as q approaches 0
419 limit_f = float(limit(f_expr, q, 0))
420 limit_df = float(limit(df_expr, q, 0))
421 limit_ddf = float(limit(ddf_expr, q, 0))
422 limit_phi_tilde_3d = float(limit(phi_tilde_3d_expr, q, 0))
423 limit_phi_tilde_3d_prime = float(limit(phi_tilde_3d_prime_expr, q, 0))
424
425 print(f"limit_f = {limit_f} | shamrock_f[0] = {shamrock_f[0]}")
426 print(f"limit_df = {limit_df} | shamrock_df[0] = {shamrock_df[0]}")
427 print(f"limit_ddf = {limit_ddf} | shamrock_ddf[0] = {shamrock_ddf[0]}")
428 print(
429 f"limit_phi_tilde_3d = {limit_phi_tilde_3d} | shamrock_phi_tilde_3d[0] = {shamrock_phi_tilde_3d[0]}"
430 )
431 print(
432 f"limit_phi_tilde_3d_prime = {limit_phi_tilde_3d_prime} | shamrock_phi_tilde_3d_prime[0] = {shamrock_phi_tilde_3d_prime[0]}"
433 )
434
435 assert abs(limit_f - shamrock_f[0]) < tolerance
436 assert abs(limit_df - shamrock_df[0]) < tolerance
437 assert abs(limit_ddf - shamrock_ddf[0]) < tolerance
438 assert abs(limit_phi_tilde_3d - shamrock_phi_tilde_3d[0]) < tolerance
439 assert abs(limit_phi_tilde_3d_prime - shamrock_phi_tilde_3d_prime[0]) < tolerance
440
441 # Compute values for q > 0 and insert limits at the beginning
442 sympy_f = [limit_f] + [f(x) for x in q_arr_no_zero]
443 sympy_df = [limit_df] + [df(x) for x in q_arr_no_zero]
444 sympy_ddf = [limit_ddf] + [ddf(x) for x in q_arr_no_zero]
445 sympy_phi_tilde_3d = [limit_phi_tilde_3d] + [phi_tilde_3d(x) for x in q_arr_no_zero]
446 sympy_phi_tilde_3d_prime = [limit_phi_tilde_3d_prime] + [
447 phi_tilde_3d_prime(x) for x in q_arr_no_zero
448 ]
449
450 # compute the absolute error
451 abs_err_f = np.max(np.abs(np.array(shamrock_f) - np.array(sympy_f)))
452 abs_err_df = np.max(np.abs(np.array(shamrock_df) - np.array(sympy_df)))
453 abs_err_ddf = np.max(np.abs(np.array(shamrock_ddf) - np.array(sympy_ddf)))
454 abs_err_phi_tilde_3d = np.max(
455 np.abs(np.array(shamrock_phi_tilde_3d) - np.array(sympy_phi_tilde_3d))
456 )
457 abs_err_phi_tilde_3d_prime = np.max(
458 np.abs(np.array(shamrock_phi_tilde_3d_prime) - np.array(sympy_phi_tilde_3d_prime))
459 )
460
461 print(f"Absolute error f(q) = {abs_err_f}")
462 print(f"Absolute error df(q) = {abs_err_df}")
463 print(f"Absolute error ddf(q) = {abs_err_ddf}")
464 print(f"Absolute error phi_tilde_3d(q) = {abs_err_phi_tilde_3d}")
465 print(f"Absolute error phi_tilde_3d_prime(q) = {abs_err_phi_tilde_3d_prime}")
466
467 assert abs_err_f < tolerance
468 assert abs_err_df < tolerance * 10
469 assert abs_err_ddf < tolerance * 100
470 assert abs_err_phi_tilde_3d < tolerance * 100
471 assert abs_err_phi_tilde_3d_prime < tolerance * 1000
472
473 print("------------------------------------------")
474 print("")
475
476 return {
477 "q_arr": q_arr,
478 "shamrock_Cf": np.array(shamrock_f) * norm_3d,
479 "shamrock_Cdf": np.array(shamrock_df) * norm_3d,
480 "shamrock_Cddf": np.array(shamrock_ddf) * norm_3d,
481 "shamrock_Cphi_tilde_3d": np.array(shamrock_phi_tilde_3d) * norm_3d,
482 "shamrock_Cphi_tilde_3d_prime": np.array(shamrock_phi_tilde_3d_prime) * norm_3d,
483 }
484
485
486 def print_kernel_info(ret):
487 print("------------------------------------------")
488 print(f"Sympy expression for kernel {ret['name']}")
489 print("------------------------------------------")
490 print(f"f(q) = {ret['f']}")
491 print(f"df(q) = {ret['df']}")
492 print(f"ddf(q) = {ret['ddf']}")
493 print(f"phi_tilde_3d(q) = {ret['phi_tilde_3d']}")
494 print(f"phi_tilde_3d_prime(q) = {ret['phi_tilde_3d_prime']}")
495 print("------------------------------------------")
496 print("")
497
498
499 def print_kernel_cpp_code(ret):
500 print("------------------------------------------")
501 print(f"C++ generated code for kernel {ret['name']}")
502 print("------------------------------------------")
503 print(ret["text"])
504 print("------------------------------------------")
505 print("")
All the SPH kernels#
Cubic splines kernels (M-series)
515 def m4():
516 Rkern = 2
517 R = sympify(Rkern)
518 f = Piecewise(
519 (sympify(1) / 4 * (R - q) ** 3 - (R / 2 - q) ** 3, q < R / 2),
520 (sympify(1) / 4 * (R - q) ** 3, q < R),
521 (0, True),
522 )
523 return (R, f, "M4")
524
525
526 def m5():
527 Rkern = sympify(5) / 2
528 R = Rkern
529 term1 = (R - q) ** 4
530 term2 = -5 * (sympify(3) / 5 * R - q) ** 4
531 term3 = 10 * (sympify(1) / 5 * R - q) ** 4
532 f = Piecewise(
533 (term1 + term2 + term3, q < sympify(1) / 5 * R),
534 (term1 + term2, q < sympify(3) / 5 * R),
535 (term1, q < R),
536 (0, True),
537 )
538 return (R, f, "M5")
539
540
541 def m6():
542 Rkern = 3
543 R = sympify(Rkern)
544 term1 = (R - q) ** 5
545 term2 = -6 * (sympify(2) / 3 * R - q) ** 5
546 term3 = 15 * (sympify(1) / 3 * R - q) ** 5
547 f = Piecewise(
548 (term1 + term2 + term3, q < sympify(1) / 3 * R),
549 (term1 + term2, q < sympify(2) / 3 * R),
550 (term1, q < R),
551 (0, True),
552 )
553 return (R, f, "M6")
554
555
556 def m7():
557 Rkern = sympify(7) / 2
558 R = Rkern
559 term1 = (R - q) ** 6
560 term2 = -7 * (sympify(5) / 7 * R - q) ** 6
561 term3 = 21 * (sympify(3) / 7 * R - q) ** 6
562 term4 = -35 * (sympify(1) / 7 * R - q) ** 6
563 f = Piecewise(
564 (term1 + term2 + term3 + term4, q < sympify(1) / 7 * R),
565 (term1 + term2 + term3, q < sympify(3) / 7 * R),
566 (term1 + term2, q < sympify(5) / 7 * R),
567 (term1, q < R),
568 (0, True),
569 )
570 return (R, f, "M7")
571
572
573 def m8():
574 Rkern = 4
575 R = sympify(Rkern)
576 term1 = (4 - q) ** 7
577 term2 = -8 * (3 - q) ** 7
578 term3 = 28 * (2 - q) ** 7
579 term4 = -56 * (1 - q) ** 7
580 f = Piecewise(
581 (term1 + term2 + term3 + term4, q < 1),
582 (term1 + term2 + term3, q < 2),
583 (term1 + term2, q < 3),
584 (term1, q < R),
585 (0, True),
586 )
587 return (R, f, "M8")
588
589
590 def m9():
591 Rkern = sympify(9) / 2
592 R = Rkern
593 term1 = (R - q) ** 8
594 term2 = -9 * (sympify(7) / 9 * R - q) ** 8
595 term3 = 36 * (sympify(5) / 9 * R - q) ** 8
596 term4 = -84 * (sympify(3) / 9 * R - q) ** 8
597 term5 = 126 * (sympify(1) / 9 * R - q) ** 8
598 f = Piecewise(
599 (term1 + term2 + term3 + term4 + term5, q < sympify(1) / 9 * R),
600 (term1 + term2 + term3 + term4, q < sympify(3) / 9 * R),
601 (term1 + term2 + term3, q < sympify(5) / 9 * R),
602 (term1 + term2, q < sympify(7) / 9 * R),
603 (term1, q < R),
604 (0, True),
605 )
606 return (R, f, "M9")
607
608
609 def m10():
610 Rkern = 5
611 R = sympify(Rkern)
612 term1 = (R - q) ** 9
613 term2 = -10 * (sympify(4) / 5 * R - q) ** 9
614 term3 = 45 * (sympify(3) / 5 * R - q) ** 9
615 term4 = -120 * (sympify(2) / 5 * R - q) ** 9
616 term5 = 210 * (sympify(1) / 5 * R - q) ** 9
617 f = Piecewise(
618 (term1 + term2 + term3 + term4 + term5, q < sympify(1) / 5 * R),
619 (term1 + term2 + term3 + term4, q < sympify(2) / 5 * R),
620 (term1 + term2 + term3, q < sympify(3) / 5 * R),
621 (term1 + term2, q < sympify(4) / 5 * R),
622 (term1, q < R),
623 (0, True),
624 )
625 return (R, f, "M10")
Wendland kernels (C-series)
630 def c2():
631 Rkern = 2
632 R = sympify(Rkern)
633 f = Piecewise(((1 - q / 2) ** 4 * (1 + 2 * q), q < R), (0, True))
634 return (R, f, "C2")
635
636
637 def c4():
638 Rkern = 2
639 R = sympify(Rkern)
640 f = Piecewise(((1 - q / 2) ** 6 * (1 + 3 * q + sympify(35) / 12 * q**2), q < R), (0, True))
641 return (R, f, "C4")
642
643
644 def c6():
645 Rkern = 2
646 R = sympify(Rkern)
647 f = Piecewise(
648 ((1 - q / 2) ** 8 * (1 + 4 * q + sympify(25) / 4 * q**2 + 4 * q**3), q < R), (0, True)
649 )
650 return (R, f, "C6")
Helper function to multiply into piecewise without expanding
655 def multiply_piecewise(piecewise_expr, factor):
656 """Multiply a factor into each piece of a Piecewise expression without expanding"""
657 new_args = []
658 for arg in piecewise_expr.args:
659 if hasattr(arg, "expr") and hasattr(arg, "cond"):
660 # Multiply the expression part, keep condition unchanged
661 new_expr = factor * arg.expr
662 new_args.append((new_expr, arg.cond))
663 elif isinstance(arg, tuple) and len(arg) == 2:
664 new_expr = factor * arg[0]
665 new_args.append((new_expr, arg[1]))
666
667 return Piecewise(*new_args)
Helper function to shift and scale a kernel
672 def shift_scale_kernel(piecewise_expr, shift_val, scale_val):
673 return piecewise_fold(
674 Piecewise(
675 (1, q < shift_val),
676 (piecewise_expr.subs({q: (q - shift_val) * scale_val}), q >= shift_val),
677 )
678 )
Double hump
683 def m4dh():
684 R, f, _ = m4()
685 f = multiply_piecewise(f, q * q)
686 return (R, f, "M4DH")
687
688
689 def m4dh3():
690 R, f, _ = m4()
691 f = multiply_piecewise(f, q * q * q)
692 return (R, f, "M4DH3")
693
694
695 def m4dh5():
696 R, f, _ = m4()
697 f = multiply_piecewise(f, q * q * q * q * q)
698 return (R, f, "M4DH5")
699
700
701 def m4dh7():
702 R, f, _ = m4()
703 f = multiply_piecewise(f, q * q * q * q * q * q * q)
704 return (R, f, "M4DH7")
M4Shift kernels
709 def m4shift2():
710 R, f, _ = m4()
711 # For q < 1: return 1
712 # For q >= 1: return M4((q - 1) * 2)
713 f_shifted = shift_scale_kernel(f, shift_val=1, scale_val=2)
714 return (R, f_shifted, "M4Shift2")
715
716
717 def m4shift4():
718 R, f, _ = m4()
719 # For q < 1.5: return 1
720 # For q >= 1.5: return M4((q - 1.5) * 4)
721 f_shifted = shift_scale_kernel(f, shift_val=sympify(3) / 2, scale_val=4)
722 return (R, f_shifted, "M4Shift4")
723
724
725 def m4shift8():
726 R, f, _ = m4()
727 # For q < 1.75: return 1
728 # For q >= 1.75: return M4((q - 1.75) * 8)
729 f_shifted = shift_scale_kernel(f, shift_val=sympify(7) / 4, scale_val=8)
730 return (R, f_shifted, "M4Shift8")
731
732
733 def m4shift16():
734 R, f, _ = m4()
735 # For q < 1.875: return 1
736 # For q >= 1.875: return M4((q - 1.875) * 16)
737 f_shifted = shift_scale_kernel(f, shift_val=sympify(15) / 8, scale_val=16)
738 return (R, f_shifted, "M4Shift16")
TGauss3 kernel (gaussian until q=3)
743 def tgauss3():
744 R = 3
745 f = Piecewise(
746 (exp(-q * q), q < R),
747 (0, True),
748 )
749 return (R, f, "TGauss3")
TGauss5 kernel
754 def tgauss5():
755 R = 5
756 f = Piecewise(
757 (exp(-q * q), q < R),
758 (0, True),
759 )
760 return (R, f, "TGauss5")
M4 Kernel#
Generate c++ code for the kernel
770 ret = kernel_to_shamrock(m4)
773 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M4
------------------------------------------
f(q) = Piecewise((-(1 - q)**3 + (2 - q)**3/4, q < 1), ((2 - q)**3/4, q < 2), (0, True))
df(q) = Piecewise((3*(1 - q)**2 - 3*(2 - q)**2/4, q < 1), (-3*(2 - q)**2/4, q < 2), (0, True))
ddf(q) = Piecewise((9*q/2 - 3, q < 1), (3 - 3*q/2, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**2*(3*q**3 - 9*q**2 + 20) - 42)/30, q < 1), (pi*(-q**6 + 9*q**5 - 30*q**4 + 40*q**3 - 48*q + 2)/(30*q), q < 2), (-pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**2*(9*q**2 - 18*q) + 2*q*(3*q**3 - 9*q**2 + 20))/30, q < 1), (pi*(-6*q**5 + 45*q**4 - 120*q**3 + 120*q**2 - 48)/(30*q) - pi*(-q**6 + 9*q**5 - 30*q**4 + 40*q**3 - 48*q + 2)/(30*q**2), q < 2), (pi/q**2, True))
------------------------------------------
776 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M4
------------------------------------------
template<class Tscal>
class KernelDefM4 {
public:
inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 2.0/3.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (10.0/7.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = 1/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<3>(1 - q);
Tscal t2 = sham::pow_constexpr<3>(2 - q);
if (q < 1) {
return -t1 + (1.0/4.0)*t2;
} else if (q < 2) {
return (1.0/4.0)*t2;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(1 - q);
Tscal t2 = sham::pow_constexpr<2>(2 - q);
if (q < 1) {
return 3*t1 - (3.0/4.0)*t2;
} else if (q < 2) {
return -(3.0/4.0)*t2;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
if (q < 1) {
return (9.0/2.0)*q - 3;
} else if (q < 2) {
return 3 - (3.0/2.0)*q;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
if (q < 1) {
return (1.0/30.0)*shambase::constants::pi<Tscal>*(t1*(3*t2 - 9*t1 + 20) - 42);
} else if (q < 2) {
return (1.0/30.0)*shambase::constants::pi<Tscal>*(-t5 + 9*t4 - 30*t3 + 40*t2 - 48*q + 2)/q;
}
else
return -shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
if (q < 1) {
return (1.0/30.0)*shambase::constants::pi<Tscal>*(t1*(9*t1 - 18*q) + 2*q*(3*t2 - 9*t1 + 20));
} else if (q < 2) {
return (1.0/30.0)*shambase::constants::pi<Tscal>*(-6*t4 + 45*t3 - 120*t2 + 120*t1 - 48)/q - (1.0/30.0)*shambase::constants::pi<Tscal>*(-t5 + 9*t4 - 30*t3 + 40*t2 - 48*q + 2)/t1;
}
else
return shambase::constants::pi<Tscal>/t1;
}
};
------------------------------------------
Test the kernel
780 test_result_m4 = test_kernel(ret)
------------------------------------------
Testing kernel M4 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.6666666666666666, expected=0.666666666666667, delta=0
Shamrock norm_2d = 0.4547284088339867, expected=0.454728408833987, delta=0
Shamrock norm_3d = 0.3183098861837907, expected=0.318309886183791, delta=0
Testing kernel radius:
Shamrock Rkern = 2.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M4_f of PyCapsule object at 0x7fb9ee6b5830>
Shamrock df(q) = <built-in method M4_df of PyCapsule object at 0x7fb9ee6b58c0>
Shamrock ddf(q) = <built-in method M4_ddf of PyCapsule object at 0x7fb9ef375920>
Shamrock phi_tilde_3d(q) = <built-in method M4_phi_tilde_3d of PyCapsule object at 0x7fb9ee6b5980>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6b59b0>
Absolute error f(q) = 3.3306690738754696212708950042725e-16
Absolute error df(q) = 6.6613381477509392425417900085449e-16
Absolute error ddf(q) = 2.7755575615628913510590791702271e-16
Absolute error phi_tilde_3d(q) = 6.1447719867801382157510209786384e-15
Absolute error phi_tilde_3d_prime(q) = 9.4052927422763077193482903413079e-15
------------------------------------------
M5 Kernel#
Generate c++ code for the kernel
789 ret = kernel_to_shamrock(m5)
792 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M5
------------------------------------------
f(q) = Piecewise((10*(1/2 - q)**4 - 5*(3/2 - q)**4 + (5/2 - q)**4, q < 1/2), (-5*(3/2 - q)**4 + (5/2 - q)**4, q < 3/2), ((5/2 - q)**4, q < 5/2), (0, True))
df(q) = Piecewise((-40*(1/2 - q)**3 + 20*(3/2 - q)**3 - 4*(5/2 - q)**3, q < 1/2), (20*(3/2 - q)**3 - 4*(5/2 - q)**3, q < 3/2), (-4*(5/2 - q)**3, q < 5/2), (0, True))
ddf(q) = Piecewise((120*(1/2 - q)**2 - 60*(3/2 - q)**2 + 12*(5/2 - q)**2, q < 1/2), (-60*(3/2 - q)**2 + 12*(5/2 - q)**2, q < 3/2), (12*(5/2 - q)**2, q < 5/2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(192*q**6 - 1008*q**4 + 3220*q**2 - 8393)/336, q < 1/2), (pi*(-128*q**7 + 896*q**6 - 2016*q**5 + 560*q**4 + 3080*q**3 - 8386*q - 1)/(336*q), q < 3/2), (pi*(64*q**7 - 896*q**6 + 5040*q**5 - 14000*q**4 + 17500*q**3 - 21875*q + 2185)/(672*q), q < 5/2), (-20*pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((pi*(1152*q**5 - 4032*q**3 + 6440*q)/336, q < 1/2), (pi*(-896*q**6 + 5376*q**5 - 10080*q**4 + 2240*q**3 + 9240*q**2 - 8386)/(336*q) - pi*(-128*q**7 + 896*q**6 - 2016*q**5 + 560*q**4 + 3080*q**3 - 8386*q - 1)/(336*q**2), q < 3/2), (pi*(448*q**6 - 5376*q**5 + 25200*q**4 - 56000*q**3 + 52500*q**2 - 21875)/(672*q) - pi*(64*q**7 - 896*q**6 + 5040*q**5 - 14000*q**4 + 17500*q**3 - 21875*q + 2185)/(672*q**2), q < 5/2), (20*pi/q**2, True))
------------------------------------------
795 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M5
------------------------------------------
template<class Tscal>
class KernelDefM5 {
public:
inline static constexpr Tscal Rkern = 5.0/2.0; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 1.0/24.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (96.0/1199.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (1.0/20.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<4>(1.0/2.0 - q);
Tscal t2 = sham::pow_constexpr<4>(3.0/2.0 - q);
Tscal t3 = sham::pow_constexpr<4>(5.0/2.0 - q);
if (q < 1.0/2.0) {
return 10*t1 - 5*t2 + t3;
} else if (q < 3.0/2.0) {
return -5*t2 + t3;
} else if (q < 5.0/2.0) {
return t3;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<3>(1.0/2.0 - q);
Tscal t2 = sham::pow_constexpr<3>(3.0/2.0 - q);
Tscal t3 = sham::pow_constexpr<3>(5.0/2.0 - q);
if (q < 1.0/2.0) {
return -40*t1 + 20*t2 - 4*t3;
} else if (q < 3.0/2.0) {
return 20*t2 - 4*t3;
} else if (q < 5.0/2.0) {
return -4*t3;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(1.0/2.0 - q);
Tscal t2 = sham::pow_constexpr<2>(3.0/2.0 - q);
Tscal t3 = sham::pow_constexpr<2>(5.0/2.0 - q);
if (q < 1.0/2.0) {
return 120*t1 - 60*t2 + 12*t3;
} else if (q < 3.0/2.0) {
return -60*t2 + 12*t3;
} else if (q < 5.0/2.0) {
return 12*t3;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
Tscal t6 = sham::pow_constexpr<7>(q);
if (q < 1.0/2.0) {
return (1.0/336.0)*shambase::constants::pi<Tscal>*(192*t5 - 1008*t3 + 3220*t1 - 8393);
} else if (q < 3.0/2.0) {
return (1.0/336.0)*shambase::constants::pi<Tscal>*(-128*t6 + 896*t5 - 2016*t4 + 560*t3 + 3080*t2 - 8386*q - 1)/q;
} else if (q < 5.0/2.0) {
return (1.0/672.0)*shambase::constants::pi<Tscal>*(64*t6 - 896*t5 + 5040*t4 - 14000*t3 + 17500*t2 - 21875*q + 2185)/q;
}
else
return -20*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
Tscal t6 = sham::pow_constexpr<7>(q);
if (q < 1.0/2.0) {
return (1.0/336.0)*shambase::constants::pi<Tscal>*(1152*t4 - 4032*t2 + 6440*q);
} else if (q < 3.0/2.0) {
return (1.0/336.0)*shambase::constants::pi<Tscal>*(-896*t5 + 5376*t4 - 10080*t3 + 2240*t2 + 9240*t1 - 8386)/q - (1.0/336.0)*shambase::constants::pi<Tscal>*(-128*t6 + 896*t5 - 2016*t4 + 560*t3 + 3080*t2 - 8386*q - 1)/t1;
} else if (q < 5.0/2.0) {
return (1.0/672.0)*shambase::constants::pi<Tscal>*(448*t5 - 5376*t4 + 25200*t3 - 56000*t2 + 52500*t1 - 21875)/q - (1.0/672.0)*shambase::constants::pi<Tscal>*(64*t6 - 896*t5 + 5040*t4 - 14000*t3 + 17500*t2 - 21875*q + 2185)/t1;
}
else
return 20*shambase::constants::pi<Tscal>/t1;
}
};
------------------------------------------
Test the kernel
799 test_result_m5 = test_kernel(ret)
------------------------------------------
Testing kernel M5 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.041666666666666664, expected=0.0416666666666667, delta=0
Shamrock norm_2d = 0.025486029252413597, expected=0.0254860292524136, delta=0
Shamrock norm_3d = 0.015915494309189534, expected=0.0159154943091895, delta=0
Testing kernel radius:
Shamrock Rkern = 2.5, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M5_f of PyCapsule object at 0x7fb9ee6b5d70>
Shamrock df(q) = <built-in method M5_df of PyCapsule object at 0x7fb9ee6b5dd0>
Shamrock ddf(q) = <built-in method M5_ddf of PyCapsule object at 0x7fb9ee6b5e60>
Shamrock phi_tilde_3d(q) = <built-in method M5_phi_tilde_3d of PyCapsule object at 0x7fb9ee6b5ec0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M5_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6b5ef0>
Absolute error f(q) = 1.1020960353947708389243550387522e-14
Absolute error df(q) = 1.6987646535410197641121131341036e-14
Absolute error ddf(q) = 2.8918153501785600372766502424727e-14
Absolute error phi_tilde_3d(q) = 2.3924264058170384135427132230976e-13
Absolute error phi_tilde_3d_prime(q) = 5.1144045415903564046637237745082e-13
------------------------------------------
M6 Kernel#
Generate c++ code for the kernel
809 ret = kernel_to_shamrock(m6)
812 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M6
------------------------------------------
f(q) = Piecewise((15*(1 - q)**5 - 6*(2 - q)**5 + (3 - q)**5, q < 1), (-6*(2 - q)**5 + (3 - q)**5, q < 2), ((3 - q)**5, q < 3), (0, True))
df(q) = Piecewise((-75*(1 - q)**4 + 30*(2 - q)**4 - 5*(3 - q)**4, q < 1), (30*(2 - q)**4 - 5*(3 - q)**4, q < 2), (-5*(3 - q)**4, q < 3), (0, True))
ddf(q) = Piecewise((300*(1 - q)**3 - 120*(2 - q)**3 + 20*(3 - q)**3, q < 1), (-120*(2 - q)**3 + 20*(3 - q)**3, q < 2), (20*(3 - q)**3, q < 3), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(-5*q**7 + 20*q**6 - 84*q**4 + 308*q**2 - 956)/7, q < 1), (pi*(5*q**8 - 60*q**7 + 280*q**6 - 588*q**5 + 350*q**4 + 476*q**3 - 1892*q - 5)/(14*q), q < 2), (pi*(-q**8 + 20*q**7 - 168*q**6 + 756*q**5 - 1890*q**4 + 2268*q**3 - 2916*q + 507)/(14*q), q < 3), (-120*pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((pi*(-35*q**6 + 120*q**5 - 336*q**3 + 616*q)/7, q < 1), (pi*(40*q**7 - 420*q**6 + 1680*q**5 - 2940*q**4 + 1400*q**3 + 1428*q**2 - 1892)/(14*q) - pi*(5*q**8 - 60*q**7 + 280*q**6 - 588*q**5 + 350*q**4 + 476*q**3 - 1892*q - 5)/(14*q**2), q < 2), (pi*(-8*q**7 + 140*q**6 - 1008*q**5 + 3780*q**4 - 7560*q**3 + 6804*q**2 - 2916)/(14*q) - pi*(-q**8 + 20*q**7 - 168*q**6 + 756*q**5 - 1890*q**4 + 2268*q**3 - 2916*q + 507)/(14*q**2), q < 3), (120*pi/q**2, True))
------------------------------------------
815 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M6
------------------------------------------
template<class Tscal>
class KernelDefM6 {
public:
inline static constexpr Tscal Rkern = 3; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 1.0/120.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (7.0/478.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (1.0/120.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<5>(1 - q);
Tscal t2 = sham::pow_constexpr<5>(2 - q);
Tscal t3 = sham::pow_constexpr<5>(3 - q);
if (q < 1) {
return 15*t1 - 6*t2 + t3;
} else if (q < 2) {
return -6*t2 + t3;
} else if (q < 3) {
return t3;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<4>(1 - q);
Tscal t2 = sham::pow_constexpr<4>(2 - q);
Tscal t3 = sham::pow_constexpr<4>(3 - q);
if (q < 1) {
return -75*t1 + 30*t2 - 5*t3;
} else if (q < 2) {
return 30*t2 - 5*t3;
} else if (q < 3) {
return -5*t3;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<3>(1 - q);
Tscal t2 = sham::pow_constexpr<3>(2 - q);
Tscal t3 = sham::pow_constexpr<3>(3 - q);
if (q < 1) {
return 300*t1 - 120*t2 + 20*t3;
} else if (q < 2) {
return -120*t2 + 20*t3;
} else if (q < 3) {
return 20*t3;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
Tscal t6 = sham::pow_constexpr<7>(q);
Tscal t7 = sham::pow_constexpr<8>(q);
if (q < 1) {
return (1.0/7.0)*shambase::constants::pi<Tscal>*(-5*t6 + 20*t5 - 84*t3 + 308*t1 - 956);
} else if (q < 2) {
return (1.0/14.0)*shambase::constants::pi<Tscal>*(5*t7 - 60*t6 + 280*t5 - 588*t4 + 350*t3 + 476*t2 - 1892*q - 5)/q;
} else if (q < 3) {
return (1.0/14.0)*shambase::constants::pi<Tscal>*(-t7 + 20*t6 - 168*t5 + 756*t4 - 1890*t3 + 2268*t2 - 2916*q + 507)/q;
}
else
return -120*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
Tscal t6 = sham::pow_constexpr<7>(q);
Tscal t7 = sham::pow_constexpr<8>(q);
if (q < 1) {
return (1.0/7.0)*shambase::constants::pi<Tscal>*(-35*t5 + 120*t4 - 336*t2 + 616*q);
} else if (q < 2) {
return (1.0/14.0)*shambase::constants::pi<Tscal>*(40*t6 - 420*t5 + 1680*t4 - 2940*t3 + 1400*t2 + 1428*t1 - 1892)/q - (1.0/14.0)*shambase::constants::pi<Tscal>*(5*t7 - 60*t6 + 280*t5 - 588*t4 + 350*t3 + 476*t2 - 1892*q - 5)/t1;
} else if (q < 3) {
return (1.0/14.0)*shambase::constants::pi<Tscal>*(-8*t6 + 140*t5 - 1008*t4 + 3780*t3 - 7560*t2 + 6804*t1 - 2916)/q - (1.0/14.0)*shambase::constants::pi<Tscal>*(-t7 + 20*t6 - 168*t5 + 756*t4 - 1890*t3 + 2268*t2 - 2916*q + 507)/t1;
}
else
return 120*shambase::constants::pi<Tscal>/t1;
}
};
------------------------------------------
Test the kernel
819 test_result_m6 = test_kernel(ret)
------------------------------------------
Testing kernel M6 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.008333333333333333, expected=0.00833333333333333, delta=0
Shamrock norm_2d = 0.00466144184787978, expected=0.00466144184787978, delta=0
Shamrock norm_3d = 0.0026525823848649226, expected=0.00265258238486492, delta=4.33680868994202E-19
Testing kernel radius:
Shamrock Rkern = 3.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M6_f of PyCapsule object at 0x7fb9ee6b61f0>
Shamrock df(q) = <built-in method M6_df of PyCapsule object at 0x7fb9ee6b6250>
Shamrock ddf(q) = <built-in method M6_ddf of PyCapsule object at 0x7fb9ee6b62e0>
Shamrock phi_tilde_3d(q) = <built-in method M6_phi_tilde_3d of PyCapsule object at 0x7fb9ee6b6310>
Shamrock phi_tilde_3d_prime(q) = <built-in method M6_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6b6340>
Absolute error f(q) = 7.105427357601002e-14
Absolute error df(q) = 1.4921397450962104e-13
Absolute error ddf(q) = 2.2737367544323206e-13
Absolute error phi_tilde_3d(q) = 3.4086883168235832315012898375205e-12
Absolute error phi_tilde_3d_prime(q) = 4.3695323848306376217663142727157e-12
------------------------------------------
M7 Kernel#
Generate c++ code for the kernel
829 ret = kernel_to_shamrock(m7)
832 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M7
------------------------------------------
f(q) = Piecewise((-35*(1/2 - q)**6 + 21*(3/2 - q)**6 - 7*(5/2 - q)**6 + (7/2 - q)**6, q < 1/2), (21*(3/2 - q)**6 - 7*(5/2 - q)**6 + (7/2 - q)**6, q < 3/2), (-7*(5/2 - q)**6 + (7/2 - q)**6, q < 5/2), ((7/2 - q)**6, q < 7/2), (0, True))
df(q) = Piecewise((210*(1/2 - q)**5 - 126*(3/2 - q)**5 + 42*(5/2 - q)**5 - 6*(7/2 - q)**5, q < 1/2), (-126*(3/2 - q)**5 + 42*(5/2 - q)**5 - 6*(7/2 - q)**5, q < 3/2), (42*(5/2 - q)**5 - 6*(7/2 - q)**5, q < 5/2), (-6*(7/2 - q)**5, q < 7/2), (0, True))
ddf(q) = Piecewise((-1050*(1/2 - q)**4 + 630*(3/2 - q)**4 - 210*(5/2 - q)**4 + 30*(7/2 - q)**4, q < 1/2), (630*(3/2 - q)**4 - 210*(5/2 - q)**4 + 30*(7/2 - q)**4, q < 3/2), (-210*(5/2 - q)**4 + 30*(7/2 - q)**4, q < 5/2), (30*(7/2 - q)**4, q < 7/2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(-1280*q**8 + 11520*q**6 - 66528*q**4 + 282576*q**2 - 1018341)/1152, q < 1/2), (pi*(3840*q**9 - 34560*q**8 + 103680*q**7 - 53760*q**6 - 235872*q**5 - 10080*q**4 + 1131984*q**3 - 4073409*q + 5)/(4608*q), q < 3/2), (pi*(-768*q**9 + 13824*q**8 - 103680*q**7 + 408576*q**6 - 852768*q**5 + 729792*q**4 + 198576*q**3 - 1948131*q - 29522)/(2304*q), q < 5/2), (pi*(256*q**9 - 6912*q**8 + 80640*q**7 - 526848*q**6 + 2074464*q**5 - 4840416*q**4 + 5647152*q**3 - 7411887*q + 1894081)/(4608*q), q < 7/2), (-840*pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((pi*(-10240*q**7 + 69120*q**5 - 266112*q**3 + 565152*q)/1152, q < 1/2), (pi*(34560*q**8 - 276480*q**7 + 725760*q**6 - 322560*q**5 - 1179360*q**4 - 40320*q**3 + 3395952*q**2 - 4073409)/(4608*q) - pi*(3840*q**9 - 34560*q**8 + 103680*q**7 - 53760*q**6 - 235872*q**5 - 10080*q**4 + 1131984*q**3 - 4073409*q + 5)/(4608*q**2), q < 3/2), (pi*(-6912*q**8 + 110592*q**7 - 725760*q**6 + 2451456*q**5 - 4263840*q**4 + 2919168*q**3 + 595728*q**2 - 1948131)/(2304*q) - pi*(-768*q**9 + 13824*q**8 - 103680*q**7 + 408576*q**6 - 852768*q**5 + 729792*q**4 + 198576*q**3 - 1948131*q - 29522)/(2304*q**2), q < 5/2), (pi*(2304*q**8 - 55296*q**7 + 564480*q**6 - 3161088*q**5 + 10372320*q**4 - 19361664*q**3 + 16941456*q**2 - 7411887)/(4608*q) - pi*(256*q**9 - 6912*q**8 + 80640*q**7 - 526848*q**6 + 2074464*q**5 - 4840416*q**4 + 5647152*q**3 - 7411887*q + 1894081)/(4608*q**2), q < 7/2), (840*pi/q**2, True))
------------------------------------------
835 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M7
------------------------------------------
template<class Tscal>
class KernelDefM7 {
public:
inline static constexpr Tscal Rkern = 7.0/2.0; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 1.0/720.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (256.0/113149.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (1.0/840.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<6>(1.0/2.0 - q);
Tscal t2 = sham::pow_constexpr<6>(3.0/2.0 - q);
Tscal t3 = sham::pow_constexpr<6>(5.0/2.0 - q);
Tscal t4 = sham::pow_constexpr<6>(7.0/2.0 - q);
if (q < 1.0/2.0) {
return -35*t1 + 21*t2 - 7*t3 + t4;
} else if (q < 3.0/2.0) {
return 21*t2 - 7*t3 + t4;
} else if (q < 5.0/2.0) {
return -7*t3 + t4;
} else if (q < 7.0/2.0) {
return t4;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<5>(1.0/2.0 - q);
Tscal t2 = sham::pow_constexpr<5>(3.0/2.0 - q);
Tscal t3 = sham::pow_constexpr<5>(5.0/2.0 - q);
Tscal t4 = sham::pow_constexpr<5>(7.0/2.0 - q);
if (q < 1.0/2.0) {
return 210*t1 - 126*t2 + 42*t3 - 6*t4;
} else if (q < 3.0/2.0) {
return -126*t2 + 42*t3 - 6*t4;
} else if (q < 5.0/2.0) {
return 42*t3 - 6*t4;
} else if (q < 7.0/2.0) {
return -6*t4;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<4>(1.0/2.0 - q);
Tscal t2 = sham::pow_constexpr<4>(3.0/2.0 - q);
Tscal t3 = sham::pow_constexpr<4>(5.0/2.0 - q);
Tscal t4 = sham::pow_constexpr<4>(7.0/2.0 - q);
if (q < 1.0/2.0) {
return -1050*t1 + 630*t2 - 210*t3 + 30*t4;
} else if (q < 3.0/2.0) {
return 630*t2 - 210*t3 + 30*t4;
} else if (q < 5.0/2.0) {
return -210*t3 + 30*t4;
} else if (q < 7.0/2.0) {
return 30*t4;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
Tscal t6 = sham::pow_constexpr<7>(q);
Tscal t7 = sham::pow_constexpr<8>(q);
Tscal t8 = sham::pow_constexpr<9>(q);
if (q < 1.0/2.0) {
return (1.0/1152.0)*shambase::constants::pi<Tscal>*(-1280*t7 + 11520*t5 - 66528*t3 + 282576*t1 - 1018341);
} else if (q < 3.0/2.0) {
return (1.0/4608.0)*shambase::constants::pi<Tscal>*(3840*t8 - 34560*t7 + 103680*t6 - 53760*t5 - 235872*t4 - 10080*t3 + 1131984*t2 - 4073409*q + 5)/q;
} else if (q < 5.0/2.0) {
return (1.0/2304.0)*shambase::constants::pi<Tscal>*(-768*t8 + 13824*t7 - 103680*t6 + 408576*t5 - 852768*t4 + 729792*t3 + 198576*t2 - 1948131*q - 29522)/q;
} else if (q < 7.0/2.0) {
return (1.0/4608.0)*shambase::constants::pi<Tscal>*(256*t8 - 6912*t7 + 80640*t6 - 526848*t5 + 2074464*t4 - 4840416*t3 + 5647152*t2 - 7411887*q + 1894081)/q;
}
else
return -840*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
Tscal t6 = sham::pow_constexpr<7>(q);
Tscal t7 = sham::pow_constexpr<8>(q);
Tscal t8 = sham::pow_constexpr<9>(q);
if (q < 1.0/2.0) {
return (1.0/1152.0)*shambase::constants::pi<Tscal>*(-10240*t6 + 69120*t4 - 266112*t2 + 565152*q);
} else if (q < 3.0/2.0) {
return (1.0/4608.0)*shambase::constants::pi<Tscal>*(34560*t7 - 276480*t6 + 725760*t5 - 322560*t4 - 1179360*t3 - 40320*t2 + 3395952*t1 - 4073409)/q - (1.0/4608.0)*shambase::constants::pi<Tscal>*(3840*t8 - 34560*t7 + 103680*t6 - 53760*t5 - 235872*t4 - 10080*t3 + 1131984*t2 - 4073409*q + 5)/t1;
} else if (q < 5.0/2.0) {
return (1.0/2304.0)*shambase::constants::pi<Tscal>*(-6912*t7 + 110592*t6 - 725760*t5 + 2451456*t4 - 4263840*t3 + 2919168*t2 + 595728*t1 - 1948131)/q - (1.0/2304.0)*shambase::constants::pi<Tscal>*(-768*t8 + 13824*t7 - 103680*t6 + 408576*t5 - 852768*t4 + 729792*t3 + 198576*t2 - 1948131*q - 29522)/t1;
} else if (q < 7.0/2.0) {
return (1.0/4608.0)*shambase::constants::pi<Tscal>*(2304*t7 - 55296*t6 + 564480*t5 - 3161088*t4 + 10372320*t3 - 19361664*t2 + 16941456*t1 - 7411887)/q - (1.0/4608.0)*shambase::constants::pi<Tscal>*(256*t8 - 6912*t7 + 80640*t6 - 526848*t5 + 2074464*t4 - 4840416*t3 + 5647152*t2 - 7411887*q + 1894081)/t1;
}
else
return 840*shambase::constants::pi<Tscal>/t1;
}
};
------------------------------------------
Test the kernel
839 test_result_m7 = test_kernel(ret, tolerance=1e-11)
------------------------------------------
Testing kernel M7 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.001388888888888889, expected=0.00138888888888889, delta=0
Shamrock norm_2d = 0.0007201772076028106, expected=0.000720177207602811, delta=0
Shamrock norm_3d = 0.0003789403406949889, expected=0.000378940340694989, delta=0
Testing kernel radius:
Shamrock Rkern = 3.5, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M7_f of PyCapsule object at 0x7fb9ee6b66d0>
Shamrock df(q) = <built-in method M7_df of PyCapsule object at 0x7fb9ee6b6760>
Shamrock ddf(q) = <built-in method M7_ddf of PyCapsule object at 0x7fb9ee6b67c0>
Shamrock phi_tilde_3d(q) = <built-in method M7_phi_tilde_3d of PyCapsule object at 0x7fb9ee6b6820>
Shamrock phi_tilde_3d_prime(q) = <built-in method M7_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6b6850>
Absolute error f(q) = 7.1169918592659551989348027572597e-13
Absolute error df(q) = 1.5356443651156504336746027078175e-12
Absolute error ddf(q) = 2.5422845301193941132174144420942e-12
Absolute error phi_tilde_3d(q) = 5.033416387056549054525265560333e-11
Absolute error phi_tilde_3d_prime(q) = 9.1132909954461041619386612291864e-11
------------------------------------------
M8 Kernel#
Generate c++ code for the kernel
848 ret = kernel_to_shamrock(m8)
851 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M8
------------------------------------------
f(q) = Piecewise((-56*(1 - q)**7 + 28*(2 - q)**7 - 8*(3 - q)**7 + (4 - q)**7, q < 1), (28*(2 - q)**7 - 8*(3 - q)**7 + (4 - q)**7, q < 2), (-8*(3 - q)**7 + (4 - q)**7, q < 3), ((4 - q)**7, q < 4), (0, True))
df(q) = Piecewise((392*(1 - q)**6 - 196*(2 - q)**6 + 56*(3 - q)**6 - 7*(4 - q)**6, q < 1), (-196*(2 - q)**6 + 56*(3 - q)**6 - 7*(4 - q)**6, q < 2), (56*(3 - q)**6 - 7*(4 - q)**6, q < 3), (-7*(4 - q)**6, q < 4), (0, True))
ddf(q) = Piecewise((-2352*(1 - q)**5 + 1176*(2 - q)**5 - 336*(3 - q)**5 + 42*(4 - q)**5, q < 1), (1176*(2 - q)**5 - 336*(3 - q)**5 + 42*(4 - q)**5, q < 2), (-336*(3 - q)**5 + 42*(4 - q)**5, q < 3), (42*(4 - q)**5, q < 4), (0, True))
phi_tilde_3d(q) = Piecewise((2*pi*(q**2*(7*q**7 - 35*q**6 + 240*q**4 - 1512*q**2 + 7248) - 29740)/9, q < 1), (2*pi*(-21*q**10 + 315*q**9 - 1890*q**8 + 5400*q**7 - 5880*q**6 - 2268*q**5 - 2940*q**4 + 37080*q**3 - 148770*q + 14)/(45*q), q < 2), (2*pi*(7*q**10 - 175*q**9 + 1890*q**8 - 11400*q**7 + 41160*q**6 - 86940*q**5 + 91140*q**4 - 16680*q**3 - 130850*q - 7154)/(45*q), q < 3), (2*pi*(-q**10 + 35*q**9 - 540*q**8 + 4800*q**7 - 26880*q**6 + 96768*q**5 - 215040*q**4 + 245760*q**3 - 327680*q + 110944)/(45*q), q < 4), (-6720*pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((2*pi*(q**2*(49*q**6 - 210*q**5 + 960*q**3 - 3024*q) + 2*q*(7*q**7 - 35*q**6 + 240*q**4 - 1512*q**2 + 7248))/9, q < 1), (2*pi*(-210*q**9 + 2835*q**8 - 15120*q**7 + 37800*q**6 - 35280*q**5 - 11340*q**4 - 11760*q**3 + 111240*q**2 - 148770)/(45*q) - 2*pi*(-21*q**10 + 315*q**9 - 1890*q**8 + 5400*q**7 - 5880*q**6 - 2268*q**5 - 2940*q**4 + 37080*q**3 - 148770*q + 14)/(45*q**2), q < 2), (2*pi*(70*q**9 - 1575*q**8 + 15120*q**7 - 79800*q**6 + 246960*q**5 - 434700*q**4 + 364560*q**3 - 50040*q**2 - 130850)/(45*q) - 2*pi*(7*q**10 - 175*q**9 + 1890*q**8 - 11400*q**7 + 41160*q**6 - 86940*q**5 + 91140*q**4 - 16680*q**3 - 130850*q - 7154)/(45*q**2), q < 3), (2*pi*(-10*q**9 + 315*q**8 - 4320*q**7 + 33600*q**6 - 161280*q**5 + 483840*q**4 - 860160*q**3 + 737280*q**2 - 327680)/(45*q) - 2*pi*(-q**10 + 35*q**9 - 540*q**8 + 4800*q**7 - 26880*q**6 + 96768*q**5 - 215040*q**4 + 245760*q**3 - 327680*q + 110944)/(45*q**2), q < 4), (6720*pi/q**2, True))
------------------------------------------
854 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M8
------------------------------------------
template<class Tscal>
class KernelDefM8 {
public:
inline static constexpr Tscal Rkern = 4; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 1.0/5040.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (9.0/29740.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (1.0/6720.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<7>(1 - q);
Tscal t2 = sham::pow_constexpr<7>(2 - q);
Tscal t3 = sham::pow_constexpr<7>(3 - q);
Tscal t4 = sham::pow_constexpr<7>(4 - q);
if (q < 1) {
return -56*t1 + 28*t2 - 8*t3 + t4;
} else if (q < 2) {
return 28*t2 - 8*t3 + t4;
} else if (q < 3) {
return -8*t3 + t4;
} else if (q < 4) {
return t4;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<6>(1 - q);
Tscal t2 = sham::pow_constexpr<6>(2 - q);
Tscal t3 = sham::pow_constexpr<6>(3 - q);
Tscal t4 = sham::pow_constexpr<6>(4 - q);
if (q < 1) {
return 392*t1 - 196*t2 + 56*t3 - 7*t4;
} else if (q < 2) {
return -196*t2 + 56*t3 - 7*t4;
} else if (q < 3) {
return 56*t3 - 7*t4;
} else if (q < 4) {
return -7*t4;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<5>(1 - q);
Tscal t2 = sham::pow_constexpr<5>(2 - q);
Tscal t3 = sham::pow_constexpr<5>(3 - q);
Tscal t4 = sham::pow_constexpr<5>(4 - q);
if (q < 1) {
return -2352*t1 + 1176*t2 - 336*t3 + 42*t4;
} else if (q < 2) {
return 1176*t2 - 336*t3 + 42*t4;
} else if (q < 3) {
return -336*t3 + 42*t4;
} else if (q < 4) {
return 42*t4;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<10>(q);
Tscal t2 = sham::pow_constexpr<2>(q);
Tscal t3 = sham::pow_constexpr<3>(q);
Tscal t4 = sham::pow_constexpr<4>(q);
Tscal t5 = sham::pow_constexpr<5>(q);
Tscal t6 = sham::pow_constexpr<6>(q);
Tscal t7 = sham::pow_constexpr<7>(q);
Tscal t8 = sham::pow_constexpr<8>(q);
Tscal t9 = sham::pow_constexpr<9>(q);
if (q < 1) {
return (2.0/9.0)*shambase::constants::pi<Tscal>*(t2*(7*t7 - 35*t6 + 240*t4 - 1512*t2 + 7248) - 29740);
} else if (q < 2) {
return (2.0/45.0)*shambase::constants::pi<Tscal>*(-21*t1 + 315*t9 - 1890*t8 + 5400*t7 - 5880*t6 - 2268*t5 - 2940*t4 + 37080*t3 - 148770*q + 14)/q;
} else if (q < 3) {
return (2.0/45.0)*shambase::constants::pi<Tscal>*(7*t1 - 175*t9 + 1890*t8 - 11400*t7 + 41160*t6 - 86940*t5 + 91140*t4 - 16680*t3 - 130850*q - 7154)/q;
} else if (q < 4) {
return (2.0/45.0)*shambase::constants::pi<Tscal>*(-t1 + 35*t9 - 540*t8 + 4800*t7 - 26880*t6 + 96768*t5 - 215040*t4 + 245760*t3 - 327680*q + 110944)/q;
}
else
return -6720*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<10>(q);
Tscal t2 = sham::pow_constexpr<2>(q);
Tscal t3 = sham::pow_constexpr<3>(q);
Tscal t4 = sham::pow_constexpr<4>(q);
Tscal t5 = sham::pow_constexpr<5>(q);
Tscal t6 = sham::pow_constexpr<6>(q);
Tscal t7 = sham::pow_constexpr<7>(q);
Tscal t8 = sham::pow_constexpr<8>(q);
Tscal t9 = sham::pow_constexpr<9>(q);
if (q < 1) {
return (2.0/9.0)*shambase::constants::pi<Tscal>*(t2*(49*t6 - 210*t5 + 960*t3 - 3024*q) + 2*q*(7*t7 - 35*t6 + 240*t4 - 1512*t2 + 7248));
} else if (q < 2) {
return (2.0/45.0)*shambase::constants::pi<Tscal>*(-210*t9 + 2835*t8 - 15120*t7 + 37800*t6 - 35280*t5 - 11340*t4 - 11760*t3 + 111240*t2 - 148770)/q - (2.0/45.0)*shambase::constants::pi<Tscal>*(-21*t1 + 315*t9 - 1890*t8 + 5400*t7 - 5880*t6 - 2268*t5 - 2940*t4 + 37080*t3 - 148770*q + 14)/t2;
} else if (q < 3) {
return (2.0/45.0)*shambase::constants::pi<Tscal>*(70*t9 - 1575*t8 + 15120*t7 - 79800*t6 + 246960*t5 - 434700*t4 + 364560*t3 - 50040*t2 - 130850)/q - (2.0/45.0)*shambase::constants::pi<Tscal>*(7*t1 - 175*t9 + 1890*t8 - 11400*t7 + 41160*t6 - 86940*t5 + 91140*t4 - 16680*t3 - 130850*q - 7154)/t2;
} else if (q < 4) {
return (2.0/45.0)*shambase::constants::pi<Tscal>*(-10*t9 + 315*t8 - 4320*t7 + 33600*t6 - 161280*t5 + 483840*t4 - 860160*t3 + 737280*t2 - 327680)/q - (2.0/45.0)*shambase::constants::pi<Tscal>*(-t1 + 35*t9 - 540*t8 + 4800*t7 - 26880*t6 + 96768*t5 - 215040*t4 + 245760*t3 - 327680*q + 110944)/t2;
}
else
return 6720*shambase::constants::pi<Tscal>/t2;
}
};
------------------------------------------
Test the kernel
858 test_result_m8 = test_kernel(ret, tolerance=1e-10)
------------------------------------------
Testing kernel M8 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.0001984126984126984, expected=0.000198412698412698, delta=0
Shamrock norm_2d = 9.63278068478183e-05, expected=0.0000963278068478183, delta=0
Shamrock norm_3d = 4.736754258687362e-05, expected=0.0000473675425868736, delta=6.77626357803440E-21
Testing kernel radius:
Shamrock Rkern = 4.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M8_f of PyCapsule object at 0x7fb9ee6b6c10>
Shamrock df(q) = <built-in method M8_df of PyCapsule object at 0x7fb9ee6b6ca0>
Shamrock ddf(q) = <built-in method M8_ddf of PyCapsule object at 0x7fb9ee6b6d00>
Shamrock phi_tilde_3d(q) = <built-in method M8_phi_tilde_3d of PyCapsule object at 0x7fb9ee6b6d60>
Shamrock phi_tilde_3d_prime(q) = <built-in method M8_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6b6d90>
Absolute error f(q) = 5.002220859751105e-12
Absolute error df(q) = 1.1368683772161603e-11
Absolute error ddf(q) = 2.9103830456733704e-11
Absolute error phi_tilde_3d(q) = 0.000000001322561202158742498480693263775
Absolute error phi_tilde_3d_prime(q) = 0.0000000018193685467461683891335635391754
------------------------------------------
M9 Kernel#
Generate c++ code for the kernel
867 ret = kernel_to_shamrock(m9)
870 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M9
------------------------------------------
f(q) = Piecewise((126*(1/2 - q)**8 - 84*(3/2 - q)**8 + 36*(5/2 - q)**8 - 9*(7/2 - q)**8 + (9/2 - q)**8, q < 1/2), (-84*(3/2 - q)**8 + 36*(5/2 - q)**8 - 9*(7/2 - q)**8 + (9/2 - q)**8, q < 3/2), (36*(5/2 - q)**8 - 9*(7/2 - q)**8 + (9/2 - q)**8, q < 5/2), (-9*(7/2 - q)**8 + (9/2 - q)**8, q < 7/2), ((9/2 - q)**8, q < 9/2), (0, True))
df(q) = Piecewise((-1008*(1/2 - q)**7 + 672*(3/2 - q)**7 - 288*(5/2 - q)**7 + 72*(7/2 - q)**7 - 8*(9/2 - q)**7, q < 1/2), (672*(3/2 - q)**7 - 288*(5/2 - q)**7 + 72*(7/2 - q)**7 - 8*(9/2 - q)**7, q < 3/2), (-288*(5/2 - q)**7 + 72*(7/2 - q)**7 - 8*(9/2 - q)**7, q < 5/2), (72*(7/2 - q)**7 - 8*(9/2 - q)**7, q < 7/2), (-8*(9/2 - q)**7, q < 9/2), (0, True))
ddf(q) = Piecewise((7056*(1/2 - q)**6 - 4704*(3/2 - q)**6 + 2016*(5/2 - q)**6 - 504*(7/2 - q)**6 + 56*(9/2 - q)**6, q < 1/2), (-4704*(3/2 - q)**6 + 2016*(5/2 - q)**6 - 504*(7/2 - q)**6 + 56*(9/2 - q)**6, q < 3/2), (2016*(5/2 - q)**6 - 504*(7/2 - q)**6 + 56*(9/2 - q)**6, q < 5/2), (-504*(7/2 - q)**6 + 56*(9/2 - q)**6, q < 7/2), (56*(9/2 - q)**6, q < 9/2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(7168*q**10 - 98560*q**8 + 908160*q**6 - 6408864*q**4 + 34283436*q**2 - 157802293)/2816, q < 1/2), (pi*(-28672*q**11 + 315392*q**10 - 1182720*q**9 + 887040*q**8 + 3801600*q**7 + 413952*q**6 - 32199552*q**5 + 36960*q**4 + 171412560*q**3 - 789011388*q - 7)/(14080*q), q < 3/2), (pi*(14336*q**11 - 315392*q**10 + 2956800*q**9 - 15079680*q**8 + 43718400*q**7 - 66646272*q**6 + 43243200*q**5 - 53850720*q**4 + 191620440*q**3 - 792042570*q + 826679)/(14080*q), q < 5/2), (pi*(-4096*q**11 + 135168*q**10 - 1971200*q**9 + 16600320*q**8 - 88281600*q**7 + 302953728*q**6 - 649756800*q**5 + 771149280*q**4 - 324004560*q**3 - 577198820*q - 96829571)/(14080*q), q < 7/2), (pi*(1024*q**11 - 45056*q**10 + 887040*q**9 - 10264320*q**8 + 76982400*q**7 - 387991296*q**6 + 1309470624*q**5 - 2806008480*q**4 + 3156759540*q**3 - 4261625379*q + 1783667601)/(28160*q), q < 9/2), (-60480*pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((pi*(71680*q**9 - 788480*q**7 + 5448960*q**5 - 25635456*q**3 + 68566872*q)/2816, q < 1/2), (pi*(-315392*q**10 + 3153920*q**9 - 10644480*q**8 + 7096320*q**7 + 26611200*q**6 + 2483712*q**5 - 160997760*q**4 + 147840*q**3 + 514237680*q**2 - 789011388)/(14080*q) - pi*(-28672*q**11 + 315392*q**10 - 1182720*q**9 + 887040*q**8 + 3801600*q**7 + 413952*q**6 - 32199552*q**5 + 36960*q**4 + 171412560*q**3 - 789011388*q - 7)/(14080*q**2), q < 3/2), (pi*(157696*q**10 - 3153920*q**9 + 26611200*q**8 - 120637440*q**7 + 306028800*q**6 - 399877632*q**5 + 216216000*q**4 - 215402880*q**3 + 574861320*q**2 - 792042570)/(14080*q) - pi*(14336*q**11 - 315392*q**10 + 2956800*q**9 - 15079680*q**8 + 43718400*q**7 - 66646272*q**6 + 43243200*q**5 - 53850720*q**4 + 191620440*q**3 - 792042570*q + 826679)/(14080*q**2), q < 5/2), (pi*(-45056*q**10 + 1351680*q**9 - 17740800*q**8 + 132802560*q**7 - 617971200*q**6 + 1817722368*q**5 - 3248784000*q**4 + 3084597120*q**3 - 972013680*q**2 - 577198820)/(14080*q) - pi*(-4096*q**11 + 135168*q**10 - 1971200*q**9 + 16600320*q**8 - 88281600*q**7 + 302953728*q**6 - 649756800*q**5 + 771149280*q**4 - 324004560*q**3 - 577198820*q - 96829571)/(14080*q**2), q < 7/2), (pi*(11264*q**10 - 450560*q**9 + 7983360*q**8 - 82114560*q**7 + 538876800*q**6 - 2327947776*q**5 + 6547353120*q**4 - 11224033920*q**3 + 9470278620*q**2 - 4261625379)/(28160*q) - pi*(1024*q**11 - 45056*q**10 + 887040*q**9 - 10264320*q**8 + 76982400*q**7 - 387991296*q**6 + 1309470624*q**5 - 2806008480*q**4 + 3156759540*q**3 - 4261625379*q + 1783667601)/(28160*q**2), q < 9/2), (60480*pi/q**2, True))
------------------------------------------
873 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M9
------------------------------------------
template<class Tscal>
class KernelDefM9 {
public:
inline static constexpr Tscal Rkern = 9.0/2.0; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 1.0/40320.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (512.0/14345663.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (1.0/60480.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<8>(1.0/2.0 - q);
Tscal t2 = sham::pow_constexpr<8>(3.0/2.0 - q);
Tscal t3 = sham::pow_constexpr<8>(5.0/2.0 - q);
Tscal t4 = sham::pow_constexpr<8>(7.0/2.0 - q);
Tscal t5 = sham::pow_constexpr<8>(9.0/2.0 - q);
if (q < 1.0/2.0) {
return 126*t1 - 84*t2 + 36*t3 - 9*t4 + t5;
} else if (q < 3.0/2.0) {
return -84*t2 + 36*t3 - 9*t4 + t5;
} else if (q < 5.0/2.0) {
return 36*t3 - 9*t4 + t5;
} else if (q < 7.0/2.0) {
return -9*t4 + t5;
} else if (q < 9.0/2.0) {
return t5;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<7>(1.0/2.0 - q);
Tscal t2 = sham::pow_constexpr<7>(3.0/2.0 - q);
Tscal t3 = sham::pow_constexpr<7>(5.0/2.0 - q);
Tscal t4 = sham::pow_constexpr<7>(7.0/2.0 - q);
Tscal t5 = sham::pow_constexpr<7>(9.0/2.0 - q);
if (q < 1.0/2.0) {
return -1008*t1 + 672*t2 - 288*t3 + 72*t4 - 8*t5;
} else if (q < 3.0/2.0) {
return 672*t2 - 288*t3 + 72*t4 - 8*t5;
} else if (q < 5.0/2.0) {
return -288*t3 + 72*t4 - 8*t5;
} else if (q < 7.0/2.0) {
return 72*t4 - 8*t5;
} else if (q < 9.0/2.0) {
return -8*t5;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<6>(1.0/2.0 - q);
Tscal t2 = sham::pow_constexpr<6>(3.0/2.0 - q);
Tscal t3 = sham::pow_constexpr<6>(5.0/2.0 - q);
Tscal t4 = sham::pow_constexpr<6>(7.0/2.0 - q);
Tscal t5 = sham::pow_constexpr<6>(9.0/2.0 - q);
if (q < 1.0/2.0) {
return 7056*t1 - 4704*t2 + 2016*t3 - 504*t4 + 56*t5;
} else if (q < 3.0/2.0) {
return -4704*t2 + 2016*t3 - 504*t4 + 56*t5;
} else if (q < 5.0/2.0) {
return 2016*t3 - 504*t4 + 56*t5;
} else if (q < 7.0/2.0) {
return -504*t4 + 56*t5;
} else if (q < 9.0/2.0) {
return 56*t5;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<10>(q);
Tscal t2 = sham::pow_constexpr<11>(q);
Tscal t3 = sham::pow_constexpr<2>(q);
Tscal t4 = sham::pow_constexpr<3>(q);
Tscal t5 = sham::pow_constexpr<4>(q);
Tscal t6 = sham::pow_constexpr<5>(q);
Tscal t7 = sham::pow_constexpr<6>(q);
Tscal t8 = sham::pow_constexpr<7>(q);
Tscal t9 = sham::pow_constexpr<8>(q);
Tscal t10 = sham::pow_constexpr<9>(q);
if (q < 1.0/2.0) {
return (1.0/2816.0)*shambase::constants::pi<Tscal>*(7168*t1 - 98560*t9 + 908160*t7 - 6408864*t5 + 34283436*t3 - 157802293);
} else if (q < 3.0/2.0) {
return (1.0/14080.0)*shambase::constants::pi<Tscal>*(-28672*t2 + 315392*t1 - 1182720*t10 + 887040*t9 + 3801600*t8 + 413952*t7 - 32199552*t6 + 36960*t5 + 171412560*t4 - 789011388*q - 7)/q;
} else if (q < 5.0/2.0) {
return (1.0/14080.0)*shambase::constants::pi<Tscal>*(14336*t2 - 315392*t1 + 2956800*t10 - 15079680*t9 + 43718400*t8 - 66646272*t7 + 43243200*t6 - 53850720*t5 + 191620440*t4 - 792042570*q + 826679)/q;
} else if (q < 7.0/2.0) {
return (1.0/14080.0)*shambase::constants::pi<Tscal>*(-4096*t2 + 135168*t1 - 1971200*t10 + 16600320*t9 - 88281600*t8 + 302953728*t7 - 649756800*t6 + 771149280*t5 - 324004560*t4 - 577198820*q - 96829571)/q;
} else if (q < 9.0/2.0) {
return (1.0/28160.0)*shambase::constants::pi<Tscal>*(1024*t2 - 45056*t1 + 887040*t10 - 10264320*t9 + 76982400*t8 - 387991296*t7 + 1309470624*t6 - 2806008480*t5 + 3156759540*t4 - 4261625379*q + 1783667601)/q;
}
else
return -60480*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<10>(q);
Tscal t2 = sham::pow_constexpr<11>(q);
Tscal t3 = sham::pow_constexpr<2>(q);
Tscal t4 = sham::pow_constexpr<3>(q);
Tscal t5 = sham::pow_constexpr<4>(q);
Tscal t6 = sham::pow_constexpr<5>(q);
Tscal t7 = sham::pow_constexpr<6>(q);
Tscal t8 = sham::pow_constexpr<7>(q);
Tscal t9 = sham::pow_constexpr<8>(q);
Tscal t10 = sham::pow_constexpr<9>(q);
if (q < 1.0/2.0) {
return (1.0/2816.0)*shambase::constants::pi<Tscal>*(71680*t10 - 788480*t8 + 5448960*t6 - 25635456*t4 + 68566872*q);
} else if (q < 3.0/2.0) {
return (1.0/14080.0)*shambase::constants::pi<Tscal>*(-315392*t1 + 3153920*t10 - 10644480*t9 + 7096320*t8 + 26611200*t7 + 2483712*t6 - 160997760*t5 + 147840*t4 + 514237680*t3 - 789011388)/q - (1.0/14080.0)*shambase::constants::pi<Tscal>*(-28672*t2 + 315392*t1 - 1182720*t10 + 887040*t9 + 3801600*t8 + 413952*t7 - 32199552*t6 + 36960*t5 + 171412560*t4 - 789011388*q - 7)/t3;
} else if (q < 5.0/2.0) {
return (1.0/14080.0)*shambase::constants::pi<Tscal>*(157696*t1 - 3153920*t10 + 26611200*t9 - 120637440*t8 + 306028800*t7 - 399877632*t6 + 216216000*t5 - 215402880*t4 + 574861320*t3 - 792042570)/q - (1.0/14080.0)*shambase::constants::pi<Tscal>*(14336*t2 - 315392*t1 + 2956800*t10 - 15079680*t9 + 43718400*t8 - 66646272*t7 + 43243200*t6 - 53850720*t5 + 191620440*t4 - 792042570*q + 826679)/t3;
} else if (q < 7.0/2.0) {
return (1.0/14080.0)*shambase::constants::pi<Tscal>*(-45056*t1 + 1351680*t10 - 17740800*t9 + 132802560*t8 - 617971200*t7 + 1817722368*t6 - 3248784000*t5 + 3084597120*t4 - 972013680*t3 - 577198820)/q - (1.0/14080.0)*shambase::constants::pi<Tscal>*(-4096*t2 + 135168*t1 - 1971200*t10 + 16600320*t9 - 88281600*t8 + 302953728*t7 - 649756800*t6 + 771149280*t5 - 324004560*t4 - 577198820*q - 96829571)/t3;
} else if (q < 9.0/2.0) {
return (1.0/28160.0)*shambase::constants::pi<Tscal>*(11264*t1 - 450560*t10 + 7983360*t9 - 82114560*t8 + 538876800*t7 - 2327947776*t6 + 6547353120*t5 - 11224033920*t4 + 9470278620*t3 - 4261625379)/q - (1.0/28160.0)*shambase::constants::pi<Tscal>*(1024*t2 - 45056*t1 + 887040*t10 - 10264320*t9 + 76982400*t8 - 387991296*t7 + 1309470624*t6 - 2806008480*t5 + 3156759540*t4 - 4261625379*q + 1783667601)/t3;
}
else
return 60480*shambase::constants::pi<Tscal>/t3;
}
};
------------------------------------------
Test the kernel
877 test_result_m9 = test_kernel(ret, tolerance=1e-9)
------------------------------------------
Testing kernel M9 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 2.48015873015873e-05, expected=0.0000248015873015873, delta=0
Shamrock norm_2d = 1.1360552783520763e-05, expected=0.0000113605527835208, delta=0
Shamrock norm_3d = 5.263060287430401e-06, expected=0.00000526306028743040, delta=0
Testing kernel radius:
Shamrock Rkern = 4.5, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M9_f of PyCapsule object at 0x7fb9ee6b7120>
Shamrock df(q) = <built-in method M9_df of PyCapsule object at 0x7fb9ee6b7180>
Shamrock ddf(q) = <built-in method M9_ddf of PyCapsule object at 0x7fb9ee6b7210>
Shamrock phi_tilde_3d(q) = <built-in method M9_phi_tilde_3d of PyCapsule object at 0x7fb9ee6b7240>
Shamrock phi_tilde_3d_prime(q) = <built-in method M9_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6b7270>
Absolute error f(q) = 1.3116584471113417810615033274653e-10
Absolute error df(q) = 1.6072615378338547344603540902465e-10
Absolute error ddf(q) = 2.873446049366578385242747112273e-10
Absolute error phi_tilde_3d(q) = 0.000000020706113693881974910981485323116
Absolute error phi_tilde_3d_prime(q) = 0.000000039433226222822340380178224869235
------------------------------------------
M10 Kernel#
Generate c++ code for the kernel
886 ret = kernel_to_shamrock(m10)
889 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M10
------------------------------------------
f(q) = Piecewise((210*(1 - q)**9 - 120*(2 - q)**9 + 45*(3 - q)**9 - 10*(4 - q)**9 + (5 - q)**9, q < 1), (-120*(2 - q)**9 + 45*(3 - q)**9 - 10*(4 - q)**9 + (5 - q)**9, q < 2), (45*(3 - q)**9 - 10*(4 - q)**9 + (5 - q)**9, q < 3), (-10*(4 - q)**9 + (5 - q)**9, q < 4), ((5 - q)**9, q < 5), (0, True))
df(q) = Piecewise((-1890*(1 - q)**8 + 1080*(2 - q)**8 - 405*(3 - q)**8 + 90*(4 - q)**8 - 9*(5 - q)**8, q < 1), (1080*(2 - q)**8 - 405*(3 - q)**8 + 90*(4 - q)**8 - 9*(5 - q)**8, q < 2), (-405*(3 - q)**8 + 90*(4 - q)**8 - 9*(5 - q)**8, q < 3), (90*(4 - q)**8 - 9*(5 - q)**8, q < 4), (-9*(5 - q)**8, q < 5), (0, True))
ddf(q) = Piecewise((15120*(1 - q)**7 - 8640*(2 - q)**7 + 3240*(3 - q)**7 - 720*(4 - q)**7 + 72*(5 - q)**7, q < 1), (-8640*(2 - q)**7 + 3240*(3 - q)**7 - 720*(4 - q)**7 + 72*(5 - q)**7, q < 2), (3240*(3 - q)**7 - 720*(4 - q)**7 + 72*(5 - q)**7, q < 3), (-720*(4 - q)**7 + 72*(5 - q)**7, q < 4), (72*(5 - q)**7, q < 5), (0, True))
phi_tilde_3d(q) = Piecewise((2*pi*(-63*q**11 + 378*q**10 - 3850*q**8 + 37620*q**6 - 291060*q**4 + 1718090*q**2 - 8766690)/33, q < 1), (2*pi*(42*q**12 - 756*q**11 + 5544*q**10 - 20020*q**9 + 31185*q**8 - 3960*q**7 + 38808*q**6 - 316008*q**5 + 10395*q**4 + 1715780*q**3 - 8766564*q - 21)/(33*q), q < 2), (2*pi*(-18*q**12 + 540*q**11 - 7128*q**10 + 53900*q**9 - 253935*q**8 + 756360*q**7 - 1380456*q**6 + 1508760*q**5 - 1510245*q**4 + 2391620*q**3 - 8914020*q + 49131)/(33*q), q < 3), (pi*(9*q**12 - 378*q**11 + 7128*q**10 - 79310*q**9 + 574695*q**8 - 2817540*q**7 + 9363816*q**6 - 20365884*q**5 + 26208765*q**4 - 14702930*q**3 - 8262102*q - 4684707)/(33*q), q < 4), (pi*(-q**12 + 54*q**11 - 1320*q**10 + 19250*q**9 - 185625*q**8 + 1237500*q**7 - 5775000*q**6 + 18562500*q**5 - 38671875*q**4 + 42968750*q**3 - 58593750*q + 28869725)/(33*q), q < 5), (-604800*pi/q, True))
phi_tilde_3d_prime(q) = Piecewise((2*pi*(-693*q**10 + 3780*q**9 - 30800*q**7 + 225720*q**5 - 1164240*q**3 + 3436180*q)/33, q < 1), (2*pi*(504*q**11 - 8316*q**10 + 55440*q**9 - 180180*q**8 + 249480*q**7 - 27720*q**6 + 232848*q**5 - 1580040*q**4 + 41580*q**3 + 5147340*q**2 - 8766564)/(33*q) - 2*pi*(42*q**12 - 756*q**11 + 5544*q**10 - 20020*q**9 + 31185*q**8 - 3960*q**7 + 38808*q**6 - 316008*q**5 + 10395*q**4 + 1715780*q**3 - 8766564*q - 21)/(33*q**2), q < 2), (2*pi*(-216*q**11 + 5940*q**10 - 71280*q**9 + 485100*q**8 - 2031480*q**7 + 5294520*q**6 - 8282736*q**5 + 7543800*q**4 - 6040980*q**3 + 7174860*q**2 - 8914020)/(33*q) - 2*pi*(-18*q**12 + 540*q**11 - 7128*q**10 + 53900*q**9 - 253935*q**8 + 756360*q**7 - 1380456*q**6 + 1508760*q**5 - 1510245*q**4 + 2391620*q**3 - 8914020*q + 49131)/(33*q**2), q < 3), (pi*(108*q**11 - 4158*q**10 + 71280*q**9 - 713790*q**8 + 4597560*q**7 - 19722780*q**6 + 56182896*q**5 - 101829420*q**4 + 104835060*q**3 - 44108790*q**2 - 8262102)/(33*q) - pi*(9*q**12 - 378*q**11 + 7128*q**10 - 79310*q**9 + 574695*q**8 - 2817540*q**7 + 9363816*q**6 - 20365884*q**5 + 26208765*q**4 - 14702930*q**3 - 8262102*q - 4684707)/(33*q**2), q < 4), (pi*(-12*q**11 + 594*q**10 - 13200*q**9 + 173250*q**8 - 1485000*q**7 + 8662500*q**6 - 34650000*q**5 + 92812500*q**4 - 154687500*q**3 + 128906250*q**2 - 58593750)/(33*q) - pi*(-q**12 + 54*q**11 - 1320*q**10 + 19250*q**9 - 185625*q**8 + 1237500*q**7 - 5775000*q**6 + 18562500*q**5 - 38671875*q**4 + 42968750*q**3 - 58593750*q + 28869725)/(33*q**2), q < 5), (604800*pi/q**2, True))
------------------------------------------
892 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M10
------------------------------------------
template<class Tscal>
class KernelDefM10 {
public:
inline static constexpr Tscal Rkern = 5; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 1.0/362880.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (11.0/2922230.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (1.0/604800.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<9>(1 - q);
Tscal t2 = sham::pow_constexpr<9>(2 - q);
Tscal t3 = sham::pow_constexpr<9>(3 - q);
Tscal t4 = sham::pow_constexpr<9>(4 - q);
Tscal t5 = sham::pow_constexpr<9>(5 - q);
if (q < 1) {
return 210*t1 - 120*t2 + 45*t3 - 10*t4 + t5;
} else if (q < 2) {
return -120*t2 + 45*t3 - 10*t4 + t5;
} else if (q < 3) {
return 45*t3 - 10*t4 + t5;
} else if (q < 4) {
return -10*t4 + t5;
} else if (q < 5) {
return t5;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<8>(1 - q);
Tscal t2 = sham::pow_constexpr<8>(2 - q);
Tscal t3 = sham::pow_constexpr<8>(3 - q);
Tscal t4 = sham::pow_constexpr<8>(4 - q);
Tscal t5 = sham::pow_constexpr<8>(5 - q);
if (q < 1) {
return -1890*t1 + 1080*t2 - 405*t3 + 90*t4 - 9*t5;
} else if (q < 2) {
return 1080*t2 - 405*t3 + 90*t4 - 9*t5;
} else if (q < 3) {
return -405*t3 + 90*t4 - 9*t5;
} else if (q < 4) {
return 90*t4 - 9*t5;
} else if (q < 5) {
return -9*t5;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<7>(1 - q);
Tscal t2 = sham::pow_constexpr<7>(2 - q);
Tscal t3 = sham::pow_constexpr<7>(3 - q);
Tscal t4 = sham::pow_constexpr<7>(4 - q);
Tscal t5 = sham::pow_constexpr<7>(5 - q);
if (q < 1) {
return 15120*t1 - 8640*t2 + 3240*t3 - 720*t4 + 72*t5;
} else if (q < 2) {
return -8640*t2 + 3240*t3 - 720*t4 + 72*t5;
} else if (q < 3) {
return 3240*t3 - 720*t4 + 72*t5;
} else if (q < 4) {
return -720*t4 + 72*t5;
} else if (q < 5) {
return 72*t5;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<10>(q);
Tscal t2 = sham::pow_constexpr<11>(q);
Tscal t3 = sham::pow_constexpr<12>(q);
Tscal t4 = sham::pow_constexpr<2>(q);
Tscal t5 = sham::pow_constexpr<3>(q);
Tscal t6 = sham::pow_constexpr<4>(q);
Tscal t7 = sham::pow_constexpr<5>(q);
Tscal t8 = sham::pow_constexpr<6>(q);
Tscal t9 = sham::pow_constexpr<7>(q);
Tscal t10 = sham::pow_constexpr<8>(q);
Tscal t11 = sham::pow_constexpr<9>(q);
if (q < 1) {
return (2.0/33.0)*shambase::constants::pi<Tscal>*(-63*t2 + 378*t1 - 3850*t10 + 37620*t8 - 291060*t6 + 1718090*t4 - 8766690);
} else if (q < 2) {
return (2.0/33.0)*shambase::constants::pi<Tscal>*(42*t3 - 756*t2 + 5544*t1 - 20020*t11 + 31185*t10 - 3960*t9 + 38808*t8 - 316008*t7 + 10395*t6 + 1715780*t5 - 8766564*q - 21)/q;
} else if (q < 3) {
return (2.0/33.0)*shambase::constants::pi<Tscal>*(-18*t3 + 540*t2 - 7128*t1 + 53900*t11 - 253935*t10 + 756360*t9 - 1380456*t8 + 1508760*t7 - 1510245*t6 + 2391620*t5 - 8914020*q + 49131)/q;
} else if (q < 4) {
return (1.0/33.0)*shambase::constants::pi<Tscal>*(9*t3 - 378*t2 + 7128*t1 - 79310*t11 + 574695*t10 - 2817540*t9 + 9363816*t8 - 20365884*t7 + 26208765*t6 - 14702930*t5 - 8262102*q - 4684707)/q;
} else if (q < 5) {
return (1.0/33.0)*shambase::constants::pi<Tscal>*(-t3 + 54*t2 - 1320*t1 + 19250*t11 - 185625*t10 + 1237500*t9 - 5775000*t8 + 18562500*t7 - 38671875*t6 + 42968750*t5 - 58593750*q + 28869725)/q;
}
else
return -604800*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<10>(q);
Tscal t2 = sham::pow_constexpr<11>(q);
Tscal t3 = sham::pow_constexpr<12>(q);
Tscal t4 = sham::pow_constexpr<2>(q);
Tscal t5 = sham::pow_constexpr<3>(q);
Tscal t6 = sham::pow_constexpr<4>(q);
Tscal t7 = sham::pow_constexpr<5>(q);
Tscal t8 = sham::pow_constexpr<6>(q);
Tscal t9 = sham::pow_constexpr<7>(q);
Tscal t10 = sham::pow_constexpr<8>(q);
Tscal t11 = sham::pow_constexpr<9>(q);
if (q < 1) {
return (2.0/33.0)*shambase::constants::pi<Tscal>*(-693*t1 + 3780*t11 - 30800*t9 + 225720*t7 - 1164240*t5 + 3436180*q);
} else if (q < 2) {
return (2.0/33.0)*shambase::constants::pi<Tscal>*(504*t2 - 8316*t1 + 55440*t11 - 180180*t10 + 249480*t9 - 27720*t8 + 232848*t7 - 1580040*t6 + 41580*t5 + 5147340*t4 - 8766564)/q - (2.0/33.0)*shambase::constants::pi<Tscal>*(42*t3 - 756*t2 + 5544*t1 - 20020*t11 + 31185*t10 - 3960*t9 + 38808*t8 - 316008*t7 + 10395*t6 + 1715780*t5 - 8766564*q - 21)/t4;
} else if (q < 3) {
return (2.0/33.0)*shambase::constants::pi<Tscal>*(-216*t2 + 5940*t1 - 71280*t11 + 485100*t10 - 2031480*t9 + 5294520*t8 - 8282736*t7 + 7543800*t6 - 6040980*t5 + 7174860*t4 - 8914020)/q - (2.0/33.0)*shambase::constants::pi<Tscal>*(-18*t3 + 540*t2 - 7128*t1 + 53900*t11 - 253935*t10 + 756360*t9 - 1380456*t8 + 1508760*t7 - 1510245*t6 + 2391620*t5 - 8914020*q + 49131)/t4;
} else if (q < 4) {
return (1.0/33.0)*shambase::constants::pi<Tscal>*(108*t2 - 4158*t1 + 71280*t11 - 713790*t10 + 4597560*t9 - 19722780*t8 + 56182896*t7 - 101829420*t6 + 104835060*t5 - 44108790*t4 - 8262102)/q - (1.0/33.0)*shambase::constants::pi<Tscal>*(9*t3 - 378*t2 + 7128*t1 - 79310*t11 + 574695*t10 - 2817540*t9 + 9363816*t8 - 20365884*t7 + 26208765*t6 - 14702930*t5 - 8262102*q - 4684707)/t4;
} else if (q < 5) {
return (1.0/33.0)*shambase::constants::pi<Tscal>*(-12*t2 + 594*t1 - 13200*t11 + 173250*t10 - 1485000*t9 + 8662500*t8 - 34650000*t7 + 92812500*t6 - 154687500*t5 + 128906250*t4 - 58593750)/q - (1.0/33.0)*shambase::constants::pi<Tscal>*(-t3 + 54*t2 - 1320*t1 + 19250*t11 - 185625*t10 + 1237500*t9 - 5775000*t8 + 18562500*t7 - 38671875*t6 + 42968750*t5 - 58593750*q + 28869725)/t4;
}
else
return 604800*shambase::constants::pi<Tscal>/t4;
}
};
------------------------------------------
Test the kernel
896 test_result_m10 = test_kernel(ret, tolerance=1e-8)
------------------------------------------
Testing kernel M10 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 2.7557319223985893e-06, expected=0.00000275573192239859, delta=0
Shamrock norm_2d = 1.1981975231318881e-06, expected=0.00000119819752313189, delta=0
Shamrock norm_3d = 5.263060287430401e-07, expected=5.26306028743040E-7, delta=1.05879118406788E-22
Testing kernel radius:
Shamrock Rkern = 5.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M10_f of PyCapsule object at 0x7fb9ee6b7630>
Shamrock df(q) = <built-in method M10_df of PyCapsule object at 0x7fb9ee6b76c0>
Shamrock ddf(q) = <built-in method M10_ddf of PyCapsule object at 0x7fb9ee8e7720>
Shamrock phi_tilde_3d(q) = <built-in method M10_phi_tilde_3d of PyCapsule object at 0x7fb9ee6b7780>
Shamrock phi_tilde_3d_prime(q) = <built-in method M10_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6b77b0>
Absolute error f(q) = 1.3969838619232178e-09
Absolute error df(q) = 2.928572939708829e-09
Absolute error ddf(q) = 3.725290298461914e-09
Absolute error phi_tilde_3d(q) = 0.00000059496615953398715061353189510027
Absolute error phi_tilde_3d_prime(q) = 0.00000075282548568718158583115382224247
------------------------------------------
C2 Kernel#
Generate c++ code for the kernel
905 ret = kernel_to_shamrock(c2)
908 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel C2
------------------------------------------
f(q) = Piecewise(((1 - q/2)**4*(2*q + 1), q < 2), (0, True))
df(q) = Piecewise((2*(1 - q/2)**4 - 2*(1 - q/2)**3*(2*q + 1), q < 2), (0, True))
ddf(q) = Piecewise((-8*(1 - q/2)**3 + 3*(1 - q/2)**2*(2*q + 1), q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**2*(3*q**5 - 30*q**4 + 112*q**3 - 168*q**2 + 224) - 384)/336, q < 2), (-16*pi/(21*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**2*(15*q**4 - 120*q**3 + 336*q**2 - 336*q) + 2*q*(3*q**5 - 30*q**4 + 112*q**3 - 168*q**2 + 224))/336, q < 2), (16*pi/(21*q**2), True))
------------------------------------------
911 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel C2
------------------------------------------
template<class Tscal>
class KernelDefC2 {
public:
inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 3.0/4.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (7.0/4.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (21.0/16.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<4>(1 - (1.0/2.0)*q);
if (q < 2) {
return t1*(2*q + 1);
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<3>(1 - (1.0/2.0)*q);
Tscal t2 = sham::pow_constexpr<4>(1 - (1.0/2.0)*q);
if (q < 2) {
return 2*t2 - 2*t1*(2*q + 1);
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(1 - (1.0/2.0)*q);
Tscal t2 = sham::pow_constexpr<3>(1 - (1.0/2.0)*q);
if (q < 2) {
return -8*t2 + 3*t1*(2*q + 1);
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
if (q < 2) {
return (1.0/336.0)*shambase::constants::pi<Tscal>*(t1*(3*t4 - 30*t3 + 112*t2 - 168*t1 + 224) - 384);
}
else
return -(16.0/21.0)*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
if (q < 2) {
return (1.0/336.0)*shambase::constants::pi<Tscal>*(t1*(15*t3 - 120*t2 + 336*t1 - 336*q) + 2*q*(3*t4 - 30*t3 + 112*t2 - 168*t1 + 224));
}
else
return (16.0/21.0)*shambase::constants::pi<Tscal>/t1;
}
};
------------------------------------------
Test the kernel
915 test_result_c2 = test_kernel(ret)
------------------------------------------
Testing kernel C2 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.75, expected=0.750000000000000, delta=0
Shamrock norm_2d = 0.5570423008216338, expected=0.557042300821634, delta=1.11022302462516E-16
Shamrock norm_3d = 0.4177817256162253, expected=0.417781725616225, delta=0
Testing kernel radius:
Shamrock Rkern = 2.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method C2_f of PyCapsule object at 0x7fb9ee6b7b10>
Shamrock df(q) = <built-in method C2_df of PyCapsule object at 0x7fb9ee6b7ba0>
Shamrock ddf(q) = <built-in method C2_ddf of PyCapsule object at 0x7fb9ee6b7c30>
Shamrock phi_tilde_3d(q) = <built-in method C2_phi_tilde_3d of PyCapsule object at 0x7fb9ee6b7c90>
Shamrock phi_tilde_3d_prime(q) = <built-in method C2_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6b7b70>
Absolute error f(q) = 4.211561406673559480149857592443e-16
Absolute error df(q) = 4.0113140431073352052618380525483e-16
Absolute error ddf(q) = 1.371769129953306271714661459359e-15
Absolute error phi_tilde_3d(q) = 4.3963002435041423632058897180998e-15
Absolute error phi_tilde_3d_prime(q) = 3.8772262773114693044442611002654e-15
------------------------------------------
C4 Kernel#
Generate c++ code for the kernel
924 ret = kernel_to_shamrock(c4)
927 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel C4
------------------------------------------
f(q) = Piecewise(((1 - q/2)**6*(35*q**2/12 + 3*q + 1), q < 2), (0, True))
df(q) = Piecewise(((1 - q/2)**6*(35*q/6 + 3) - 3*(1 - q/2)**5*(35*q**2/12 + 3*q + 1), q < 2), (0, True))
ddf(q) = Piecewise((35*(1 - q/2)**6/6 - 6*(1 - q/2)**5*(35*q/6 + 3) + 15*(1 - q/2)**4*(35*q**2/12 + 3*q + 1)/2, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**2*(105*q**8 - 1408*q**7 + 7700*q**6 - 21120*q**5 + 26400*q**4 - 29568*q**2 + 42240) - 56320)/63360, q < 2), (-256*pi/(495*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**2*(840*q**7 - 9856*q**6 + 46200*q**5 - 105600*q**4 + 105600*q**3 - 59136*q) + 2*q*(105*q**8 - 1408*q**7 + 7700*q**6 - 21120*q**5 + 26400*q**4 - 29568*q**2 + 42240))/63360, q < 2), (256*pi/(495*q**2), True))
------------------------------------------
930 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel C4
------------------------------------------
template<class Tscal>
class KernelDefC4 {
public:
inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 27.0/32.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (9.0/4.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (495.0/256.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<6>(1 - (1.0/2.0)*q);
if (q < 2) {
return t2*((35.0/12.0)*t1 + 3*q + 1);
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<5>(1 - (1.0/2.0)*q);
Tscal t3 = sham::pow_constexpr<6>(1 - (1.0/2.0)*q);
if (q < 2) {
return t3*((35.0/6.0)*q + 3) - 3*t2*((35.0/12.0)*t1 + 3*q + 1);
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<4>(1 - (1.0/2.0)*q);
Tscal t3 = sham::pow_constexpr<5>(1 - (1.0/2.0)*q);
Tscal t4 = sham::pow_constexpr<6>(1 - (1.0/2.0)*q);
if (q < 2) {
return (35.0/6.0)*t4 - 6*t3*((35.0/6.0)*q + 3) + (15.0/2.0)*t2*((35.0/12.0)*t1 + 3*q + 1);
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<4>(q);
Tscal t3 = sham::pow_constexpr<5>(q);
Tscal t4 = sham::pow_constexpr<6>(q);
Tscal t5 = sham::pow_constexpr<7>(q);
Tscal t6 = sham::pow_constexpr<8>(q);
if (q < 2) {
return (1.0/63360.0)*shambase::constants::pi<Tscal>*(t1*(105*t6 - 1408*t5 + 7700*t4 - 21120*t3 + 26400*t2 - 29568*t1 + 42240) - 56320);
}
else
return -(256.0/495.0)*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
Tscal t6 = sham::pow_constexpr<7>(q);
Tscal t7 = sham::pow_constexpr<8>(q);
if (q < 2) {
return (1.0/63360.0)*shambase::constants::pi<Tscal>*(t1*(840*t6 - 9856*t5 + 46200*t4 - 105600*t3 + 105600*t2 - 59136*q) + 2*q*(105*t7 - 1408*t6 + 7700*t5 - 21120*t4 + 26400*t3 - 29568*t1 + 42240));
}
else
return (256.0/495.0)*shambase::constants::pi<Tscal>/t1;
}
};
------------------------------------------
Test the kernel
934 test_result_c4 = test_kernel(ret)
------------------------------------------
Testing kernel C4 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.84375, expected=0.843750000000000, delta=0
Shamrock norm_2d = 0.716197243913529, expected=0.716197243913529, delta=0
Shamrock norm_3d = 0.6154820064881891, expected=0.615482006488189, delta=1.11022302462516E-16
Testing kernel radius:
Shamrock Rkern = 2.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method C4_f of PyCapsule object at 0x7fb9ee6d0090>
Shamrock df(q) = <built-in method C4_df of PyCapsule object at 0x7fb9ee6d0120>
Shamrock ddf(q) = <built-in method C4_ddf of PyCapsule object at 0x7fb9ee6d01b0>
Shamrock phi_tilde_3d(q) = <built-in method C4_phi_tilde_3d of PyCapsule object at 0x7fb9ee6d0210>
Shamrock phi_tilde_3d_prime(q) = <built-in method C4_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6d00f0>
Absolute error f(q) = 6.77057790166676168827606387144e-16
Absolute error df(q) = 7.3058119683808204154226284388856e-16
Absolute error ddf(q) = 4.5862635469006640758973444472878e-15
Absolute error phi_tilde_3d(q) = 2.9691588069996695982673972652105e-14
Absolute error phi_tilde_3d_prime(q) = 6.4917526524589430038090522815696e-14
------------------------------------------
C6 Kernel#
Generate c++ code for the kernel
943 ret = kernel_to_shamrock(c6)
946 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel C6
------------------------------------------
f(q) = Piecewise(((1 - q/2)**8*(4*q**3 + 25*q**2/4 + 4*q + 1), q < 2), (0, True))
df(q) = Piecewise(((1 - q/2)**8*(12*q**2 + 25*q/2 + 4) - 4*(1 - q/2)**7*(4*q**3 + 25*q**2/4 + 4*q + 1), q < 2), (0, True))
ddf(q) = Piecewise(((1 - q/2)**8*(24*q + 25/2) - 8*(1 - q/2)**7*(12*q**2 + 25*q/2 + 4) + 14*(1 - q/2)**6*(4*q**3 + 25*q**2/4 + 4*q + 1), q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**2*(480*q**11 - 8085*q**10 + 58240*q**9 - 229320*q**8 + 512512*q**7 - 560560*q**6 + 549120*q**4 - 768768*q**2 + 931840) - 1003520)/1397760, q < 2), (-512*pi/(1365*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**2*(5280*q**10 - 80850*q**9 + 524160*q**8 - 1834560*q**7 + 3587584*q**6 - 3363360*q**5 + 2196480*q**3 - 1537536*q) + 2*q*(480*q**11 - 8085*q**10 + 58240*q**9 - 229320*q**8 + 512512*q**7 - 560560*q**6 + 549120*q**4 - 768768*q**2 + 931840))/1397760, q < 2), (512*pi/(1365*q**2), True))
------------------------------------------
949 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel C6
------------------------------------------
template<class Tscal>
class KernelDefC6 {
public:
inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 15.0/16.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (39.0/14.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (1365.0/512.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<8>(1 - (1.0/2.0)*q);
if (q < 2) {
return t3*(4*t2 + (25.0/4.0)*t1 + 4*q + 1);
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<7>(1 - (1.0/2.0)*q);
Tscal t4 = sham::pow_constexpr<8>(1 - (1.0/2.0)*q);
if (q < 2) {
return t4*(12*t1 + (25.0/2.0)*q + 4) - 4*t3*(4*t2 + (25.0/4.0)*t1 + 4*q + 1);
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<6>(1 - (1.0/2.0)*q);
Tscal t4 = sham::pow_constexpr<7>(1 - (1.0/2.0)*q);
Tscal t5 = sham::pow_constexpr<8>(1 - (1.0/2.0)*q);
if (q < 2) {
return t5*(24*q + 25.0/2.0) - 8*t4*(12*t1 + (25.0/2.0)*q + 4) + 14*t3*(4*t2 + (25.0/4.0)*t1 + 4*q + 1);
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<10>(q);
Tscal t2 = sham::pow_constexpr<11>(q);
Tscal t3 = sham::pow_constexpr<2>(q);
Tscal t4 = sham::pow_constexpr<4>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
Tscal t6 = sham::pow_constexpr<7>(q);
Tscal t7 = sham::pow_constexpr<8>(q);
Tscal t8 = sham::pow_constexpr<9>(q);
if (q < 2) {
return (1.0/1397760.0)*shambase::constants::pi<Tscal>*(t3*(480*t2 - 8085*t1 + 58240*t8 - 229320*t7 + 512512*t6 - 560560*t5 + 549120*t4 - 768768*t3 + 931840) - 1003520);
}
else
return -(512.0/1365.0)*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<10>(q);
Tscal t2 = sham::pow_constexpr<11>(q);
Tscal t3 = sham::pow_constexpr<2>(q);
Tscal t4 = sham::pow_constexpr<3>(q);
Tscal t5 = sham::pow_constexpr<4>(q);
Tscal t6 = sham::pow_constexpr<5>(q);
Tscal t7 = sham::pow_constexpr<6>(q);
Tscal t8 = sham::pow_constexpr<7>(q);
Tscal t9 = sham::pow_constexpr<8>(q);
Tscal t10 = sham::pow_constexpr<9>(q);
if (q < 2) {
return (1.0/1397760.0)*shambase::constants::pi<Tscal>*(t3*(5280*t1 - 80850*t10 + 524160*t9 - 1834560*t8 + 3587584*t7 - 3363360*t6 + 2196480*t4 - 1537536*q) + 2*q*(480*t2 - 8085*t1 + 58240*t10 - 229320*t9 + 512512*t8 - 560560*t7 + 549120*t5 - 768768*t3 + 931840));
}
else
return (512.0/1365.0)*shambase::constants::pi<Tscal>/t3;
}
};
------------------------------------------
Test the kernel
953 test_result_c6 = test_kernel(ret)
------------------------------------------
Testing kernel C6 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.9375, expected=0.937500000000000, delta=0
Shamrock norm_2d = 0.886720397226274, expected=0.886720397226274, delta=0
Shamrock norm_3d = 0.8486191301579575, expected=0.848619130157958, delta=0
Testing kernel radius:
Shamrock Rkern = 2.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method C6_f of PyCapsule object at 0x7fb9ee6d05d0>
Shamrock df(q) = <built-in method C6_df of PyCapsule object at 0x7fb9ee6d0660>
Shamrock ddf(q) = <built-in method C6_ddf of PyCapsule object at 0x7fb9ee6d06f0>
Shamrock phi_tilde_3d(q) = <built-in method C6_phi_tilde_3d of PyCapsule object at 0x7fb9ee6d0750>
Shamrock phi_tilde_3d_prime(q) = <built-in method C6_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6d0780>
Absolute error f(q) = 8.0028218067564718410675074309679e-16
Absolute error df(q) = 1.2403417665354433235366686954591e-15
Absolute error ddf(q) = 8.6752247783722774983167486346621e-15
Absolute error phi_tilde_3d(q) = 1.5252718357359013011480267646023e-13
Absolute error phi_tilde_3d_prime(q) = 3.809628034813140135293395651618e-13
------------------------------------------
M4DH Kernel#
Generate c++ code for the kernel
963 ret = kernel_to_shamrock(m4dh)
966 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M4DH
------------------------------------------
f(q) = Piecewise((q**2*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (q**2*(2 - q)**3/4, q < 2), (0, True))
df(q) = Piecewise((q**2*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 2*q*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**2*(2 - q)**2/4 + q*(2 - q)**3/2, q < 2), (0, True))
ddf(q) = Piecewise((q**2*(9*q/2 - 3) + 4*q*(3*(1 - q)**2 - 3*(2 - q)**2/4) - 2*(1 - q)**3 + (2 - q)**3/2, q < 1), (-3*q**2*(2*q - 4)/4 - 3*q*(2 - q)**2 + (2 - q)**3/2, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**4*(15*q**3 - 40*q**2 + 56) - 248)/280, q < 1), (pi*(-5*q**8 + 40*q**7 - 112*q**6 + 112*q**5 - 256*q + 4)/(280*q), q < 2), (-9*pi/(10*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**4*(45*q**2 - 80*q) + 4*q**3*(15*q**3 - 40*q**2 + 56))/280, q < 1), (pi*(-40*q**7 + 280*q**6 - 672*q**5 + 560*q**4 - 256)/(280*q) - pi*(-5*q**8 + 40*q**7 - 112*q**6 + 112*q**5 - 256*q + 4)/(280*q**2), q < 2), (9*pi/(10*q**2), True))
------------------------------------------
969 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M4DH
------------------------------------------
template<class Tscal>
class KernelDefM4DH {
public:
inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 2;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (70.0/31.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (10.0/9.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(1 - q);
Tscal t3 = sham::pow_constexpr<3>(2 - q);
if (q < 1) {
return t1*(-t2 + (1.0/4.0)*t3);
} else if (q < 2) {
return (1.0/4.0)*t1*t3;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(1 - q);
Tscal t2 = sham::pow_constexpr<2>(2 - q);
Tscal t3 = sham::pow_constexpr<2>(q);
Tscal t4 = sham::pow_constexpr<3>(1 - q);
Tscal t5 = sham::pow_constexpr<3>(2 - q);
if (q < 1) {
return t3*(3*t1 - (3.0/4.0)*t2) + 2*q*(-t4 + (1.0/4.0)*t5);
} else if (q < 2) {
return -(3.0/4.0)*t3*t2 + (1.0/2.0)*q*t5;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(1 - q);
Tscal t2 = sham::pow_constexpr<2>(2 - q);
Tscal t3 = sham::pow_constexpr<2>(q);
Tscal t4 = sham::pow_constexpr<3>(1 - q);
Tscal t5 = sham::pow_constexpr<3>(2 - q);
if (q < 1) {
return t3*((9.0/2.0)*q - 3) + 4*q*(3*t1 - (3.0/4.0)*t2) - 2*t4 + (1.0/2.0)*t5;
} else if (q < 2) {
return -(3.0/4.0)*t3*(2*q - 4) - 3*q*t2 + (1.0/2.0)*t5;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
Tscal t6 = sham::pow_constexpr<7>(q);
Tscal t7 = sham::pow_constexpr<8>(q);
if (q < 1) {
return (1.0/280.0)*shambase::constants::pi<Tscal>*(t3*(15*t2 - 40*t1 + 56) - 248);
} else if (q < 2) {
return (1.0/280.0)*shambase::constants::pi<Tscal>*(-5*t7 + 40*t6 - 112*t5 + 112*t4 - 256*q + 4)/q;
}
else
return -(9.0/10.0)*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
Tscal t6 = sham::pow_constexpr<7>(q);
Tscal t7 = sham::pow_constexpr<8>(q);
if (q < 1) {
return (1.0/280.0)*shambase::constants::pi<Tscal>*(t3*(45*t1 - 80*q) + 4*t2*(15*t2 - 40*t1 + 56));
} else if (q < 2) {
return (1.0/280.0)*shambase::constants::pi<Tscal>*(-40*t6 + 280*t5 - 672*t4 + 560*t3 - 256)/q - (1.0/280.0)*shambase::constants::pi<Tscal>*(-5*t7 + 40*t6 - 112*t5 + 112*t4 - 256*q + 4)/t1;
}
else
return (9.0/10.0)*shambase::constants::pi<Tscal>/t1;
}
};
------------------------------------------
Test the kernel
973 test_result_m4dh = test_kernel(ret)
------------------------------------------
Testing kernel M4DH matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 2.0, expected=2.00000000000000, delta=0
Shamrock norm_2d = 0.7187642591246886, expected=0.718764259124689, delta=1.11022302462516E-16
Shamrock norm_3d = 0.35367765131532297, expected=0.353677651315323, delta=0
Testing kernel radius:
Shamrock Rkern = 2.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M4DH_f of PyCapsule object at 0x7fb9ee6d0ae0>
Shamrock df(q) = <built-in method M4DH_df of PyCapsule object at 0x7fb9ee6d0b10>
Shamrock ddf(q) = <built-in method M4DH_ddf of PyCapsule object at 0x7fb9ee6d0b40>
Shamrock phi_tilde_3d(q) = <built-in method M4DH_phi_tilde_3d of PyCapsule object at 0x7fb9ee6d0b70>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4DH_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6d0ba0>
Absolute error f(q) = 1.2402696483592525566706857128555e-16
Absolute error df(q) = 1.1457310014385301187257368207776e-15
Absolute error ddf(q) = 9.0591089023422617007004687927574e-16
Absolute error phi_tilde_3d(q) = 8.3097552699191310126976481785064e-15
Absolute error phi_tilde_3d_prime(q) = 3.3207246830111135096394910341455e-14
------------------------------------------
M4DH3 Kernel#
Generate c++ code for the kernel
983 ret = kernel_to_shamrock(m4dh3)
986 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M4DH3
------------------------------------------
f(q) = Piecewise((q**3*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (q**3*(2 - q)**3/4, q < 2), (0, True))
df(q) = Piecewise((q**3*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 3*q**2*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**3*(2 - q)**2/4 + 3*q**2*(2 - q)**3/4, q < 2), (0, True))
ddf(q) = Piecewise((q**3*(9*q/2 - 3) + 6*q**2*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 6*q*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**3*(2*q - 4)/4 - 9*q**2*(2 - q)**2/2 + 3*q*(2 - q)**3/2, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**5*(35*q**3 - 90*q**2 + 112) - 756)/840, q < 1), (pi*(-35*q**9 + 270*q**8 - 720*q**7 + 672*q**6 - 2304*q + 20)/(2520*q), q < 2), (-127*pi/(126*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**5*(105*q**2 - 180*q) + 5*q**4*(35*q**3 - 90*q**2 + 112))/840, q < 1), (pi*(-315*q**8 + 2160*q**7 - 5040*q**6 + 4032*q**5 - 2304)/(2520*q) - pi*(-35*q**9 + 270*q**8 - 720*q**7 + 672*q**6 - 2304*q + 20)/(2520*q**2), q < 2), (127*pi/(126*q**2), True))
------------------------------------------
989 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M4DH3
------------------------------------------
template<class Tscal>
class KernelDefM4DH3 {
public:
inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 70.0/31.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (20.0/9.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (126.0/127.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<3>(1 - q);
Tscal t2 = sham::pow_constexpr<3>(2 - q);
Tscal t3 = sham::pow_constexpr<3>(q);
if (q < 1) {
return t3*(-t1 + (1.0/4.0)*t2);
} else if (q < 2) {
return (1.0/4.0)*t3*t2;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(1 - q);
Tscal t2 = sham::pow_constexpr<2>(2 - q);
Tscal t3 = sham::pow_constexpr<2>(q);
Tscal t4 = sham::pow_constexpr<3>(1 - q);
Tscal t5 = sham::pow_constexpr<3>(2 - q);
Tscal t6 = sham::pow_constexpr<3>(q);
if (q < 1) {
return t6*(3*t1 - (3.0/4.0)*t2) + 3*t3*(-t4 + (1.0/4.0)*t5);
} else if (q < 2) {
return -(3.0/4.0)*t6*t2 + (3.0/4.0)*t3*t5;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(1 - q);
Tscal t2 = sham::pow_constexpr<2>(2 - q);
Tscal t3 = sham::pow_constexpr<2>(q);
Tscal t4 = sham::pow_constexpr<3>(1 - q);
Tscal t5 = sham::pow_constexpr<3>(2 - q);
Tscal t6 = sham::pow_constexpr<3>(q);
if (q < 1) {
return t6*((9.0/2.0)*q - 3) + 6*t3*(3*t1 - (3.0/4.0)*t2) + 6*q*(-t4 + (1.0/4.0)*t5);
} else if (q < 2) {
return -(3.0/4.0)*t6*(2*q - 4) - (9.0/2.0)*t3*t2 + (3.0/2.0)*q*t5;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<5>(q);
Tscal t4 = sham::pow_constexpr<6>(q);
Tscal t5 = sham::pow_constexpr<7>(q);
Tscal t6 = sham::pow_constexpr<8>(q);
Tscal t7 = sham::pow_constexpr<9>(q);
if (q < 1) {
return (1.0/840.0)*shambase::constants::pi<Tscal>*(t3*(35*t2 - 90*t1 + 112) - 756);
} else if (q < 2) {
return (1.0/2520.0)*shambase::constants::pi<Tscal>*(-35*t7 + 270*t6 - 720*t5 + 672*t4 - 2304*q + 20)/q;
}
else
return -(127.0/126.0)*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
Tscal t6 = sham::pow_constexpr<7>(q);
Tscal t7 = sham::pow_constexpr<8>(q);
Tscal t8 = sham::pow_constexpr<9>(q);
if (q < 1) {
return (1.0/840.0)*shambase::constants::pi<Tscal>*(t4*(105*t1 - 180*q) + 5*t3*(35*t2 - 90*t1 + 112));
} else if (q < 2) {
return (1.0/2520.0)*shambase::constants::pi<Tscal>*(-315*t7 + 2160*t6 - 5040*t5 + 4032*t4 - 2304)/q - (1.0/2520.0)*shambase::constants::pi<Tscal>*(-35*t8 + 270*t7 - 720*t6 + 672*t5 - 2304*q + 20)/t1;
}
else
return (127.0/126.0)*shambase::constants::pi<Tscal>/t1;
}
};
------------------------------------------
Test the kernel
993 test_result_m4dh3 = test_kernel(ret)
------------------------------------------
Testing kernel M4DH3 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 2.2580645161290325, expected=2.25806451612903, delta=4.44089209850063E-16
Shamrock norm_2d = 0.7073553026306459, expected=0.707355302630646, delta=0
Shamrock norm_3d = 0.31580350912722543, expected=0.315803509127225, delta=5.55111512312578E-17
Testing kernel radius:
Shamrock Rkern = 2.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M4DH3_f of PyCapsule object at 0x7fb9ee6d0db0>
Shamrock df(q) = <built-in method M4DH3_df of PyCapsule object at 0x7fb9ee6d0e10>
Shamrock ddf(q) = <built-in method M4DH3_ddf of PyCapsule object at 0x7fb9ee6d0e40>
Shamrock phi_tilde_3d(q) = <built-in method M4DH3_phi_tilde_3d of PyCapsule object at 0x7fb9ee6d0e70>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4DH3_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6d0ea0>
Absolute error f(q) = 1.2218718934561889186295748586557e-16
Absolute error df(q) = 2.2764428614529119995212313147607e-15
Absolute error ddf(q) = 1.0320055293124922608435806326406e-15
Absolute error phi_tilde_3d(q) = 1.8637731717387722473783472141369e-14
Absolute error phi_tilde_3d_prime(q) = 4.2032553794902827882340545209492e-14
------------------------------------------
M4DH5 Kernel#
Generate c++ code for the kernel
1003 ret = kernel_to_shamrock(m4dh5)
1006 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M4DH5
------------------------------------------
f(q) = Piecewise((q**5*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (q**5*(2 - q)**3/4, q < 2), (0, True))
df(q) = Piecewise((q**5*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 5*q**4*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**5*(2 - q)**2/4 + 5*q**4*(2 - q)**3/4, q < 2), (0, True))
ddf(q) = Piecewise((q**5*(9*q/2 - 3) + 10*q**4*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 20*q**3*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**5*(2*q - 4)/4 - 15*q**4*(2 - q)**2/2 + 5*q**3*(2 - q)**3, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**7*(63*q**3 - 154*q**2 + 165) - 2805)/2310, q < 1), (pi*(-21*q**11 + 154*q**10 - 385*q**9 + 330*q**8 - 2816*q + 7)/(2310*q), q < 2), (-511*pi/(330*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**7*(189*q**2 - 308*q) + 7*q**6*(63*q**3 - 154*q**2 + 165))/2310, q < 1), (pi*(-231*q**10 + 1540*q**9 - 3465*q**8 + 2640*q**7 - 2816)/(2310*q) - pi*(-21*q**11 + 154*q**10 - 385*q**9 + 330*q**8 - 2816*q + 7)/(2310*q**2), q < 2), (511*pi/(330*q**2), True))
------------------------------------------
1009 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M4DH5
------------------------------------------
template<class Tscal>
class KernelDefM4DH5 {
public:
inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 252.0/127.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (28.0/17.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (330.0/511.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<3>(1 - q);
Tscal t2 = sham::pow_constexpr<3>(2 - q);
Tscal t3 = sham::pow_constexpr<5>(q);
if (q < 1) {
return t3*(-t1 + (1.0/4.0)*t2);
} else if (q < 2) {
return (1.0/4.0)*t3*t2;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(1 - q);
Tscal t2 = sham::pow_constexpr<2>(2 - q);
Tscal t3 = sham::pow_constexpr<3>(1 - q);
Tscal t4 = sham::pow_constexpr<3>(2 - q);
Tscal t5 = sham::pow_constexpr<4>(q);
Tscal t6 = sham::pow_constexpr<5>(q);
if (q < 1) {
return t6*(3*t1 - (3.0/4.0)*t2) + 5*t5*(-t3 + (1.0/4.0)*t4);
} else if (q < 2) {
return -(3.0/4.0)*t6*t2 + (5.0/4.0)*t5*t4;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(1 - q);
Tscal t2 = sham::pow_constexpr<2>(2 - q);
Tscal t3 = sham::pow_constexpr<3>(1 - q);
Tscal t4 = sham::pow_constexpr<3>(2 - q);
Tscal t5 = sham::pow_constexpr<3>(q);
Tscal t6 = sham::pow_constexpr<4>(q);
Tscal t7 = sham::pow_constexpr<5>(q);
if (q < 1) {
return t7*((9.0/2.0)*q - 3) + 10*t6*(3*t1 - (3.0/4.0)*t2) + 20*t5*(-t3 + (1.0/4.0)*t4);
} else if (q < 2) {
return -(3.0/4.0)*t7*(2*q - 4) - (15.0/2.0)*t6*t2 + 5*t5*t4;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<10>(q);
Tscal t2 = sham::pow_constexpr<11>(q);
Tscal t3 = sham::pow_constexpr<2>(q);
Tscal t4 = sham::pow_constexpr<3>(q);
Tscal t5 = sham::pow_constexpr<7>(q);
Tscal t6 = sham::pow_constexpr<8>(q);
Tscal t7 = sham::pow_constexpr<9>(q);
if (q < 1) {
return (1.0/2310.0)*shambase::constants::pi<Tscal>*(t5*(63*t4 - 154*t3 + 165) - 2805);
} else if (q < 2) {
return (1.0/2310.0)*shambase::constants::pi<Tscal>*(-21*t2 + 154*t1 - 385*t7 + 330*t6 - 2816*q + 7)/q;
}
else
return -(511.0/330.0)*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<10>(q);
Tscal t2 = sham::pow_constexpr<11>(q);
Tscal t3 = sham::pow_constexpr<2>(q);
Tscal t4 = sham::pow_constexpr<3>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
Tscal t6 = sham::pow_constexpr<7>(q);
Tscal t7 = sham::pow_constexpr<8>(q);
Tscal t8 = sham::pow_constexpr<9>(q);
if (q < 1) {
return (1.0/2310.0)*shambase::constants::pi<Tscal>*(t6*(189*t3 - 308*q) + 7*t5*(63*t4 - 154*t3 + 165));
} else if (q < 2) {
return (1.0/2310.0)*shambase::constants::pi<Tscal>*(-231*t1 + 1540*t8 - 3465*t7 + 2640*t6 - 2816)/q - (1.0/2310.0)*shambase::constants::pi<Tscal>*(-21*t2 + 154*t1 - 385*t8 + 330*t7 - 2816*q + 7)/t3;
}
else
return (511.0/330.0)*shambase::constants::pi<Tscal>/t3;
}
};
------------------------------------------
Test the kernel
1013 test_result_m4dh5 = test_kernel(ret)
------------------------------------------
Testing kernel M4DH5 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 1.984251968503937, expected=1.98425196850394, delta=0
Shamrock norm_2d = 0.5242751066556552, expected=0.524275106655655, delta=0
Shamrock norm_3d = 0.20556215741810357, expected=0.205562157418104, delta=0
Testing kernel radius:
Shamrock Rkern = 2.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M4DH5_f of PyCapsule object at 0x7fb9ee6d10e0>
Shamrock df(q) = <built-in method M4DH5_df of PyCapsule object at 0x7fb9ee6d1110>
Shamrock ddf(q) = <built-in method M4DH5_ddf of PyCapsule object at 0x7fb9ee6d1140>
Shamrock phi_tilde_3d(q) = <built-in method M4DH5_phi_tilde_3d of PyCapsule object at 0x7fb9ee6d1170>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4DH5_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6d11a0>
Absolute error f(q) = 1.1325334002180934793019450988886e-16
Absolute error df(q) = 8.9872253031621204204154022587469e-15
Absolute error ddf(q) = 3.0518560430000101442400695924271e-15
Absolute error phi_tilde_3d(q) = 4.139306340840759038659513118751e-14
Absolute error phi_tilde_3d_prime(q) = 2.1996238639070152921228687986712e-13
------------------------------------------
M4DH7 Kernel#
Generate c++ code for the kernel
1023 ret = kernel_to_shamrock(m4dh7)
1026 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M4DH7
------------------------------------------
f(q) = Piecewise((q**7*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (q**7*(2 - q)**3/4, q < 2), (0, True))
df(q) = Piecewise((q**7*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 7*q**6*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**7*(2 - q)**2/4 + 7*q**6*(2 - q)**3/4, q < 2), (0, True))
ddf(q) = Piecewise((q**7*(9*q/2 - 3) + 14*q**6*(3*(1 - q)**2 - 3*(2 - q)**2/4) + 42*q**5*(-(1 - q)**3 + (2 - q)**3/4), q < 1), (-3*q**7*(2*q - 4)/4 - 21*q**6*(2 - q)**2/2 + 21*q**5*(2 - q)**3/2, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(q**9*(495*q**3 - 1170*q**2 + 1144) - 53196)/25740, q < 1), (pi*(-165*q**13 + 1170*q**12 - 2808*q**11 + 2288*q**10 - 53248*q + 36)/(25740*q), q < 2), (-2047*pi/(715*q), True))
phi_tilde_3d_prime(q) = Piecewise((pi*(q**9*(1485*q**2 - 2340*q) + 9*q**8*(495*q**3 - 1170*q**2 + 1144))/25740, q < 1), (pi*(-2145*q**12 + 14040*q**11 - 30888*q**10 + 22880*q**9 - 53248)/(25740*q) - pi*(-165*q**13 + 1170*q**12 - 2808*q**11 + 2288*q**10 - 53248*q + 36)/(25740*q**2), q < 2), (2047*pi/(715*q**2), True))
------------------------------------------
1029 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M4DH7
------------------------------------------
template<class Tscal>
class KernelDefM4DH7 {
public:
inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 660.0/511.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (30.0/31.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (715.0/2047.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<3>(1 - q);
Tscal t2 = sham::pow_constexpr<3>(2 - q);
Tscal t3 = sham::pow_constexpr<7>(q);
if (q < 1) {
return t3*(-t1 + (1.0/4.0)*t2);
} else if (q < 2) {
return (1.0/4.0)*t3*t2;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(1 - q);
Tscal t2 = sham::pow_constexpr<2>(2 - q);
Tscal t3 = sham::pow_constexpr<3>(1 - q);
Tscal t4 = sham::pow_constexpr<3>(2 - q);
Tscal t5 = sham::pow_constexpr<6>(q);
Tscal t6 = sham::pow_constexpr<7>(q);
if (q < 1) {
return t6*(3*t1 - (3.0/4.0)*t2) + 7*t5*(-t3 + (1.0/4.0)*t4);
} else if (q < 2) {
return -(3.0/4.0)*t6*t2 + (7.0/4.0)*t5*t4;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(1 - q);
Tscal t2 = sham::pow_constexpr<2>(2 - q);
Tscal t3 = sham::pow_constexpr<3>(1 - q);
Tscal t4 = sham::pow_constexpr<3>(2 - q);
Tscal t5 = sham::pow_constexpr<5>(q);
Tscal t6 = sham::pow_constexpr<6>(q);
Tscal t7 = sham::pow_constexpr<7>(q);
if (q < 1) {
return t7*((9.0/2.0)*q - 3) + 14*t6*(3*t1 - (3.0/4.0)*t2) + 42*t5*(-t3 + (1.0/4.0)*t4);
} else if (q < 2) {
return -(3.0/4.0)*t7*(2*q - 4) - (21.0/2.0)*t6*t2 + (21.0/2.0)*t5*t4;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<10>(q);
Tscal t2 = sham::pow_constexpr<11>(q);
Tscal t3 = sham::pow_constexpr<12>(q);
Tscal t4 = sham::pow_constexpr<13>(q);
Tscal t5 = sham::pow_constexpr<2>(q);
Tscal t6 = sham::pow_constexpr<3>(q);
Tscal t7 = sham::pow_constexpr<9>(q);
if (q < 1) {
return (1.0/25740.0)*shambase::constants::pi<Tscal>*(t7*(495*t6 - 1170*t5 + 1144) - 53196);
} else if (q < 2) {
return (1.0/25740.0)*shambase::constants::pi<Tscal>*(-165*t4 + 1170*t3 - 2808*t2 + 2288*t1 - 53248*q + 36)/q;
}
else
return -(2047.0/715.0)*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<10>(q);
Tscal t2 = sham::pow_constexpr<11>(q);
Tscal t3 = sham::pow_constexpr<12>(q);
Tscal t4 = sham::pow_constexpr<13>(q);
Tscal t5 = sham::pow_constexpr<2>(q);
Tscal t6 = sham::pow_constexpr<3>(q);
Tscal t7 = sham::pow_constexpr<8>(q);
Tscal t8 = sham::pow_constexpr<9>(q);
if (q < 1) {
return (1.0/25740.0)*shambase::constants::pi<Tscal>*(t8*(1485*t5 - 2340*q) + 9*t7*(495*t6 - 1170*t5 + 1144));
} else if (q < 2) {
return (1.0/25740.0)*shambase::constants::pi<Tscal>*(-2145*t3 + 14040*t2 - 30888*t1 + 22880*t8 - 53248)/q - (1.0/25740.0)*shambase::constants::pi<Tscal>*(-165*t4 + 1170*t3 - 2808*t2 + 2288*t1 - 53248*q + 36)/t5;
}
else
return (2047.0/715.0)*shambase::constants::pi<Tscal>/t5;
}
};
------------------------------------------
Test the kernel
1033 test_result_m4dh7 = test_kernel(ret)
------------------------------------------
Testing kernel M4DH7 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 1.2915851272015655, expected=1.29158512720157, delta=0
Shamrock norm_2d = 0.3080418253391523, expected=0.308041825339152, delta=5.55111512312578E-17
Shamrock norm_3d = 0.11118298418241834, expected=0.111182984182418, delta=1.38777878078145E-17
Testing kernel radius:
Shamrock Rkern = 2.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M4DH7_f of PyCapsule object at 0x7fb9ee6d1410>
Shamrock df(q) = <built-in method M4DH7_df of PyCapsule object at 0x7fb9ee6d1470>
Shamrock ddf(q) = <built-in method M4DH7_ddf of PyCapsule object at 0x7fb9ee6d14a0>
Shamrock phi_tilde_3d(q) = <built-in method M4DH7_phi_tilde_3d of PyCapsule object at 0x7fb9ee6d14d0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4DH7_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6d1500>
Absolute error f(q) = 2.7745488273644668238686811902901e-16
Absolute error df(q) = 3.5482969028749994264124678817261e-14
Absolute error ddf(q) = 8.9764994763351031453229854274657e-15
Absolute error phi_tilde_3d(q) = 1.1410837367852384316621894214667e-13
Absolute error phi_tilde_3d_prime(q) = 4.2978977992024483366058950511763e-13
------------------------------------------
M4Shift2 Kernel#
Generate c++ code for the kernel
1043 ret = kernel_to_shamrock(m4shift2)
1046 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M4Shift2
------------------------------------------
f(q) = Piecewise((1, q < 1), (-(3 - 2*q)**3 + (4 - 2*q)**3/4, q < 3/2), ((4 - 2*q)**3/4, q < 2), (0, True))
df(q) = Piecewise((0, q < 1), (6*(3 - 2*q)**2 - 3*(4 - 2*q)**2/2, q < 3/2), (-3*(4 - 2*q)**2/2, q < 2), (0, True))
ddf(q) = Piecewise((0, q < 1), (36*q - 48, q < 3/2), (24 - 12*q, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(40*q**2 - 231)/60, q < 1), (pi*(48*q**6 - 288*q**5 + 600*q**4 - 440*q**3 - 39*q - 72)/(60*q), q < 3/2), (pi*(-32*q**6 + 288*q**5 - 960*q**4 + 1280*q**3 - 1536*q + 585)/(120*q), q < 2), (-439*pi/(120*q), True))
phi_tilde_3d_prime(q) = Piecewise((4*pi*q/3, q < 1), (pi*(288*q**5 - 1440*q**4 + 2400*q**3 - 1320*q**2 - 39)/(60*q) - pi*(48*q**6 - 288*q**5 + 600*q**4 - 440*q**3 - 39*q - 72)/(60*q**2), q < 3/2), (pi*(-192*q**5 + 1440*q**4 - 3840*q**3 + 3840*q**2 - 1536)/(120*q) - pi*(-32*q**6 + 288*q**5 - 960*q**4 + 1280*q**3 - 1536*q + 585)/(120*q**2), q < 2), (439*pi/(120*q**2), True))
------------------------------------------
1049 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M4Shift2
------------------------------------------
template<class Tscal>
class KernelDefM4Shift2 {
public:
inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 4.0/11.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (40.0/77.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (120.0/439.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<3>(3 - 2*q);
Tscal t2 = sham::pow_constexpr<3>(4 - 2*q);
if (q < 1) {
return 1;
} else if (q < 3.0/2.0) {
return -t1 + (1.0/4.0)*t2;
} else if (q < 2) {
return (1.0/4.0)*t2;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(3 - 2*q);
Tscal t2 = sham::pow_constexpr<2>(4 - 2*q);
if (q < 1) {
return 0;
} else if (q < 3.0/2.0) {
return 6*t1 - (3.0/2.0)*t2;
} else if (q < 2) {
return -(3.0/2.0)*t2;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
if (q < 1) {
return 0;
} else if (q < 3.0/2.0) {
return 36*q - 48;
} else if (q < 2) {
return 24 - 12*q;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
if (q < 1) {
return (1.0/60.0)*shambase::constants::pi<Tscal>*(40*t1 - 231);
} else if (q < 3.0/2.0) {
return (1.0/60.0)*shambase::constants::pi<Tscal>*(48*t5 - 288*t4 + 600*t3 - 440*t2 - 39*q - 72)/q;
} else if (q < 2) {
return (1.0/120.0)*shambase::constants::pi<Tscal>*(-32*t5 + 288*t4 - 960*t3 + 1280*t2 - 1536*q + 585)/q;
}
else
return -(439.0/120.0)*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
if (q < 1) {
return (4.0/3.0)*shambase::constants::pi<Tscal>*q;
} else if (q < 3.0/2.0) {
return (1.0/60.0)*shambase::constants::pi<Tscal>*(288*t4 - 1440*t3 + 2400*t2 - 1320*t1 - 39)/q - (1.0/60.0)*shambase::constants::pi<Tscal>*(48*t5 - 288*t4 + 600*t3 - 440*t2 - 39*q - 72)/t1;
} else if (q < 2) {
return (1.0/120.0)*shambase::constants::pi<Tscal>*(-192*t4 + 1440*t3 - 3840*t2 + 3840*t1 - 1536)/q - (1.0/120.0)*shambase::constants::pi<Tscal>*(-32*t5 + 288*t4 - 960*t3 + 1280*t2 - 1536*q + 585)/t1;
}
else
return (439.0/120.0)*shambase::constants::pi<Tscal>/t1;
}
};
------------------------------------------
Test the kernel
1053 test_result_m4shift2 = test_kernel(ret)
------------------------------------------
Testing kernel M4Shift2 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.36363636363636365, expected=0.363636363636364, delta=0
Shamrock norm_2d = 0.1653557850305406, expected=0.165355785030541, delta=0
Shamrock norm_3d = 0.08700953608668538, expected=0.0870095360866854, delta=0
Testing kernel radius:
Shamrock Rkern = 2.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M4Shift2_f of PyCapsule object at 0x7fb9ee6d1740>
Shamrock df(q) = <built-in method M4Shift2_df of PyCapsule object at 0x7fb9ee6d1770>
Shamrock ddf(q) = <built-in method M4Shift2_ddf of PyCapsule object at 0x7fb9ee6d17a0>
Shamrock phi_tilde_3d(q) = <built-in method M4Shift2_phi_tilde_3d of PyCapsule object at 0x7fb9ee6d17d0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4Shift2_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6d1800>
Absolute error f(q) = 3.3306690738754696212708950042725e-16
Absolute error df(q) = 9.7144514654701197287067770957947e-16
Absolute error ddf(q) = 1.7763568394002505e-15
Absolute error phi_tilde_3d(q) = 4.8057488080825312374425800768037e-14
Absolute error phi_tilde_3d_prime(q) = 8.7604806542352871750438831278774e-14
------------------------------------------
M4Shift4 Kernel#
Generate c++ code for the kernel
1063 ret = kernel_to_shamrock(m4shift4)
1066 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M4Shift4
------------------------------------------
f(q) = Piecewise((1, q < 3/2), (-(7 - 4*q)**3 + (8 - 4*q)**3/4, q < 7/4), ((8 - 4*q)**3/4, q < 2), (0, True))
df(q) = Piecewise((0, q < 3/2), (12*(7 - 4*q)**2 - 3*(8 - 4*q)**2, q < 7/4), (-3*(8 - 4*q)**2, q < 2), (0, True))
ddf(q) = Piecewise((0, q < 3/2), (288*q - 480, q < 7/4), (192 - 96*q, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(160*q**2 - 1371)/240, q < 3/2), (pi*(1536*q**6 - 11520*q**5 + 31680*q**4 - 34400*q**3 + 25845*q - 14580)/(240*q), q < 7/4), (pi*(-2048*q**6 + 18432*q**5 - 61440*q**4 + 81920*q**3 - 98304*q + 59329)/(960*q), q < 2), (-2069*pi/(320*q), True))
phi_tilde_3d_prime(q) = Piecewise((4*pi*q/3, q < 3/2), (pi*(9216*q**5 - 57600*q**4 + 126720*q**3 - 103200*q**2 + 25845)/(240*q) - pi*(1536*q**6 - 11520*q**5 + 31680*q**4 - 34400*q**3 + 25845*q - 14580)/(240*q**2), q < 7/4), (pi*(-12288*q**5 + 92160*q**4 - 245760*q**3 + 245760*q**2 - 98304)/(960*q) - pi*(-2048*q**6 + 18432*q**5 - 61440*q**4 + 81920*q**3 - 98304*q + 59329)/(960*q**2), q < 2), (2069*pi/(320*q**2), True))
------------------------------------------
1069 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M4Shift4
------------------------------------------
template<class Tscal>
class KernelDefM4Shift4 {
public:
inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 8.0/27.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (160.0/457.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (320.0/2069.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<3>(7 - 4*q);
Tscal t2 = sham::pow_constexpr<3>(8 - 4*q);
if (q < 3.0/2.0) {
return 1;
} else if (q < 7.0/4.0) {
return -t1 + (1.0/4.0)*t2;
} else if (q < 2) {
return (1.0/4.0)*t2;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(7 - 4*q);
Tscal t2 = sham::pow_constexpr<2>(8 - 4*q);
if (q < 3.0/2.0) {
return 0;
} else if (q < 7.0/4.0) {
return 12*t1 - 3*t2;
} else if (q < 2) {
return -3*t2;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
if (q < 3.0/2.0) {
return 0;
} else if (q < 7.0/4.0) {
return 288*q - 480;
} else if (q < 2) {
return 192 - 96*q;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
if (q < 3.0/2.0) {
return (1.0/240.0)*shambase::constants::pi<Tscal>*(160*t1 - 1371);
} else if (q < 7.0/4.0) {
return (1.0/240.0)*shambase::constants::pi<Tscal>*(1536*t5 - 11520*t4 + 31680*t3 - 34400*t2 + 25845*q - 14580)/q;
} else if (q < 2) {
return (1.0/960.0)*shambase::constants::pi<Tscal>*(-2048*t5 + 18432*t4 - 61440*t3 + 81920*t2 - 98304*q + 59329)/q;
}
else
return -(2069.0/320.0)*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
if (q < 3.0/2.0) {
return (4.0/3.0)*shambase::constants::pi<Tscal>*q;
} else if (q < 7.0/4.0) {
return (1.0/240.0)*shambase::constants::pi<Tscal>*(9216*t4 - 57600*t3 + 126720*t2 - 103200*t1 + 25845)/q - (1.0/240.0)*shambase::constants::pi<Tscal>*(1536*t5 - 11520*t4 + 31680*t3 - 34400*t2 + 25845*q - 14580)/t1;
} else if (q < 2) {
return (1.0/960.0)*shambase::constants::pi<Tscal>*(-12288*t4 + 92160*t3 - 245760*t2 + 245760*t1 - 98304)/q - (1.0/960.0)*shambase::constants::pi<Tscal>*(-2048*t5 + 18432*t4 - 61440*t3 + 81920*t2 - 98304*q + 59329)/t1;
}
else
return (2069.0/320.0)*shambase::constants::pi<Tscal>/t1;
}
};
------------------------------------------
Test the kernel
1073 test_result_m4shift4 = test_kernel(ret)
------------------------------------------
Testing kernel M4Shift4 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.2962962962962963, expected=0.296296296296296, delta=0
Shamrock norm_2d = 0.11144328619126151, expected=0.111443286191262, delta=0
Shamrock norm_3d = 0.049231108544617215, expected=0.0492311085446172, delta=0
Testing kernel radius:
Shamrock Rkern = 2.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M4Shift4_f of PyCapsule object at 0x7fb9ee6d1a40>
Shamrock df(q) = <built-in method M4Shift4_df of PyCapsule object at 0x7fb9ee6d1a70>
Shamrock ddf(q) = <built-in method M4Shift4_ddf of PyCapsule object at 0x7fb9ee6d1aa0>
Shamrock phi_tilde_3d(q) = <built-in method M4Shift4_phi_tilde_3d of PyCapsule object at 0x7fb9ee6d1ad0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4Shift4_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6d1b00>
Absolute error f(q) = 2.7755575615628913510590791702271e-16
Absolute error df(q) = 2.220446049250313e-15
Absolute error ddf(q) = 2.842170943040401e-14
Absolute error phi_tilde_3d(q) = 4.4156542197657458622817453458733e-13
Absolute error phi_tilde_3d_prime(q) = 1.0539828786480343310056803196174e-12
------------------------------------------
M4Shift8 Kernel#
Generate c++ code for the kernel
1083 ret = kernel_to_shamrock(m4shift8)
1086 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M4Shift8
------------------------------------------
f(q) = Piecewise((1, q < 7/4), (-(15 - 8*q)**3 + (16 - 8*q)**3/4, q < 15/8), ((16 - 8*q)**3/4, q < 2), (0, True))
df(q) = Piecewise((0, q < 7/4), (24*(15 - 8*q)**2 - 6*(16 - 8*q)**2, q < 15/8), (-6*(16 - 8*q)**2, q < 2), (0, True))
ddf(q) = Piecewise((0, q < 7/4), (2304*q - 4224, q < 15/8), (1536 - 768*q, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(640*q**2 - 6531)/960, q < 7/4), (pi*(49152*q**6 - 405504*q**5 + 1236480*q**4 - 1504640*q**3 + 1491693*q - 907578)/(960*q), q < 15/8), (pi*(-131072*q**6 + 1179648*q**5 - 3932160*q**4 + 5242880*q**3 - 6291456*q + 4130001)/(7680*q), q < 2), (-64303*pi/(7680*q), True))
phi_tilde_3d_prime(q) = Piecewise((4*pi*q/3, q < 7/4), (pi*(294912*q**5 - 2027520*q**4 + 4945920*q**3 - 4513920*q**2 + 1491693)/(960*q) - pi*(49152*q**6 - 405504*q**5 + 1236480*q**4 - 1504640*q**3 + 1491693*q - 907578)/(960*q**2), q < 15/8), (pi*(-786432*q**5 + 5898240*q**4 - 15728640*q**3 + 15728640*q**2 - 6291456)/(7680*q) - pi*(-131072*q**6 + 1179648*q**5 - 3932160*q**4 + 5242880*q**3 - 6291456*q + 4130001)/(7680*q**2), q < 2), (64303*pi/(7680*q**2), True))
------------------------------------------
1089 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M4Shift8
------------------------------------------
template<class Tscal>
class KernelDefM4Shift8 {
public:
inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 16.0/59.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (640.0/2177.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (7680.0/64303.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<3>(15 - 8*q);
Tscal t2 = sham::pow_constexpr<3>(16 - 8*q);
if (q < 7.0/4.0) {
return 1;
} else if (q < 15.0/8.0) {
return -t1 + (1.0/4.0)*t2;
} else if (q < 2) {
return (1.0/4.0)*t2;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(15 - 8*q);
Tscal t2 = sham::pow_constexpr<2>(16 - 8*q);
if (q < 7.0/4.0) {
return 0;
} else if (q < 15.0/8.0) {
return 24*t1 - 6*t2;
} else if (q < 2) {
return -6*t2;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
if (q < 7.0/4.0) {
return 0;
} else if (q < 15.0/8.0) {
return 2304*q - 4224;
} else if (q < 2) {
return 1536 - 768*q;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
if (q < 7.0/4.0) {
return (1.0/960.0)*shambase::constants::pi<Tscal>*(640*t1 - 6531);
} else if (q < 15.0/8.0) {
return (1.0/960.0)*shambase::constants::pi<Tscal>*(49152*t5 - 405504*t4 + 1236480*t3 - 1504640*t2 + 1491693*q - 907578)/q;
} else if (q < 2) {
return (1.0/7680.0)*shambase::constants::pi<Tscal>*(-131072*t5 + 1179648*t4 - 3932160*t3 + 5242880*t2 - 6291456*q + 4130001)/q;
}
else
return -(64303.0/7680.0)*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
if (q < 7.0/4.0) {
return (4.0/3.0)*shambase::constants::pi<Tscal>*q;
} else if (q < 15.0/8.0) {
return (1.0/960.0)*shambase::constants::pi<Tscal>*(294912*t4 - 2027520*t3 + 4945920*t2 - 4513920*t1 + 1491693)/q - (1.0/960.0)*shambase::constants::pi<Tscal>*(49152*t5 - 405504*t4 + 1236480*t3 - 1504640*t2 + 1491693*q - 907578)/t1;
} else if (q < 2) {
return (1.0/7680.0)*shambase::constants::pi<Tscal>*(-786432*t4 + 5898240*t3 - 15728640*t2 + 15728640*t1 - 6291456)/q - (1.0/7680.0)*shambase::constants::pi<Tscal>*(-131072*t5 + 1179648*t4 - 3932160*t3 + 5242880*t2 - 6291456*q + 4130001)/t1;
}
else
return (64303.0/7680.0)*shambase::constants::pi<Tscal>/t1;
}
};
------------------------------------------
Test the kernel
1093 test_result_m4shift8 = test_kernel(ret)
------------------------------------------
Testing kernel M4Shift8 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.2711864406779661, expected=0.271186440677966, delta=0
Shamrock norm_2d = 0.09357755037098117, expected=0.0935775503709812, delta=1.38777878078145E-17
Shamrock norm_3d = 0.03801719866711526, expected=0.0380171986671153, delta=6.93889390390723E-18
Testing kernel radius:
Shamrock Rkern = 2.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M4Shift8_f of PyCapsule object at 0x7fb9ee6d1d40>
Shamrock df(q) = <built-in method M4Shift8_df of PyCapsule object at 0x7fb9ee6d1d70>
Shamrock ddf(q) = <built-in method M4Shift8_ddf of PyCapsule object at 0x7fb9ee6d1da0>
Shamrock phi_tilde_3d(q) = <built-in method M4Shift8_phi_tilde_3d of PyCapsule object at 0x7fb9ee6d1dd0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4Shift8_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6d1e00>
Absolute error f(q) = 3.3306690738754696212708950042725e-16
Absolute error df(q) = 5.773159728050814e-15
Absolute error ddf(q) = 2.8421709430404007e-13
Absolute error phi_tilde_3d(q) = 4.8963625190753593596772964897216e-12
Absolute error phi_tilde_3d_prime(q) = 8.8836344431557297504000270626971e-12
------------------------------------------
M4Shift16 Kernel#
Generate c++ code for the kernel
1103 ret = kernel_to_shamrock(m4shift16)
1106 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel M4Shift16
------------------------------------------
f(q) = Piecewise((1, q < 15/8), (-(31 - 16*q)**3 + (32 - 16*q)**3/4, q < 31/16), ((32 - 16*q)**3/4, q < 2), (0, True))
df(q) = Piecewise((0, q < 15/8), (48*(31 - 16*q)**2 - 12*(32 - 16*q)**2, q < 31/16), (-12*(32 - 16*q)**2, q < 2), (0, True))
ddf(q) = Piecewise((0, q < 15/8), (18432*q - 35328, q < 31/16), (12288 - 6144*q, q < 2), (0, True))
phi_tilde_3d(q) = Piecewise((pi*(2560*q**2 - 28371)/3840, q < 15/8), (pi*(1572864*q**6 - 13565952*q**5 + 43315200*q**4 - 55293440*q**3 + 60721629*q - 38728125)/(3840*q), q < 31/16), (pi*(-8388608*q**6 + 75497472*q**5 - 251658240*q**4 + 335544320*q**3 - 402653184*q + 267853681)/(61440*q), q < 2), (-38785*pi/(4096*q), True))
phi_tilde_3d_prime(q) = Piecewise((4*pi*q/3, q < 15/8), (pi*(9437184*q**5 - 67829760*q**4 + 173260800*q**3 - 165880320*q**2 + 60721629)/(3840*q) - pi*(1572864*q**6 - 13565952*q**5 + 43315200*q**4 - 55293440*q**3 + 60721629*q - 38728125)/(3840*q**2), q < 31/16), (pi*(-50331648*q**5 + 377487360*q**4 - 1006632960*q**3 + 1006632960*q**2 - 402653184)/(61440*q) - pi*(-8388608*q**6 + 75497472*q**5 - 251658240*q**4 + 335544320*q**3 - 402653184*q + 267853681)/(61440*q**2), q < 2), (38785*pi/(4096*q**2), True))
------------------------------------------
1109 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel M4Shift16
------------------------------------------
template<class Tscal>
class KernelDefM4Shift16 {
public:
inline static constexpr Tscal Rkern = 2; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 32.0/123.0;
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = (2560.0/9457.0)/shambase::constants::pi<Tscal>;
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (4096.0/38785.0)/shambase::constants::pi<Tscal>;
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<3>(31 - 16*q);
Tscal t2 = sham::pow_constexpr<3>(32 - 16*q);
if (q < 15.0/8.0) {
return 1;
} else if (q < 31.0/16.0) {
return -t1 + (1.0/4.0)*t2;
} else if (q < 2) {
return (1.0/4.0)*t2;
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(31 - 16*q);
Tscal t2 = sham::pow_constexpr<2>(32 - 16*q);
if (q < 15.0/8.0) {
return 0;
} else if (q < 31.0/16.0) {
return 48*t1 - 12*t2;
} else if (q < 2) {
return -12*t2;
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
if (q < 15.0/8.0) {
return 0;
} else if (q < 31.0/16.0) {
return 18432*q - 35328;
} else if (q < 2) {
return 12288 - 6144*q;
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
if (q < 15.0/8.0) {
return (1.0/3840.0)*shambase::constants::pi<Tscal>*(2560*t1 - 28371);
} else if (q < 31.0/16.0) {
return (1.0/3840.0)*shambase::constants::pi<Tscal>*(1572864*t5 - 13565952*t4 + 43315200*t3 - 55293440*t2 + 60721629*q - 38728125)/q;
} else if (q < 2) {
return (1.0/61440.0)*shambase::constants::pi<Tscal>*(-8388608*t5 + 75497472*t4 - 251658240*t3 + 335544320*t2 - 402653184*q + 267853681)/q;
}
else
return -(38785.0/4096.0)*shambase::constants::pi<Tscal>/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
Tscal t2 = sham::pow_constexpr<3>(q);
Tscal t3 = sham::pow_constexpr<4>(q);
Tscal t4 = sham::pow_constexpr<5>(q);
Tscal t5 = sham::pow_constexpr<6>(q);
if (q < 15.0/8.0) {
return (4.0/3.0)*shambase::constants::pi<Tscal>*q;
} else if (q < 31.0/16.0) {
return (1.0/3840.0)*shambase::constants::pi<Tscal>*(9437184*t4 - 67829760*t3 + 173260800*t2 - 165880320*t1 + 60721629)/q - (1.0/3840.0)*shambase::constants::pi<Tscal>*(1572864*t5 - 13565952*t4 + 43315200*t3 - 55293440*t2 + 60721629*q - 38728125)/t1;
} else if (q < 2) {
return (1.0/61440.0)*shambase::constants::pi<Tscal>*(-50331648*t4 + 377487360*t3 - 1006632960*t2 + 1006632960*t1 - 402653184)/q - (1.0/61440.0)*shambase::constants::pi<Tscal>*(-8388608*t5 + 75497472*t4 - 251658240*t3 + 335544320*t2 - 402653184*q + 267853681)/t1;
}
else
return (38785.0/4096.0)*shambase::constants::pi<Tscal>/t1;
}
};
------------------------------------------
Test the kernel
1113 test_result_m4shift16 = test_kernel(ret)
------------------------------------------
Testing kernel M4Shift16 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.2601626016260163, expected=0.260162601626016, delta=5.55111512312578E-17
Shamrock norm_2d = 0.08616615296928246, expected=0.0861661529692825, delta=0
Shamrock norm_3d = 0.03361601891991251, expected=0.0336160189199125, delta=6.93889390390723E-18
Testing kernel radius:
Shamrock Rkern = 2.0, delta=0
Testing kernel functions:
Shamrock f(q) = <built-in method M4Shift16_f of PyCapsule object at 0x7fb9ee6d2040>
Shamrock df(q) = <built-in method M4Shift16_df of PyCapsule object at 0x7fb9ee6d2070>
Shamrock ddf(q) = <built-in method M4Shift16_ddf of PyCapsule object at 0x7fb9ee6d20a0>
Shamrock phi_tilde_3d(q) = <built-in method M4Shift16_phi_tilde_3d of PyCapsule object at 0x7fb9ee6d20d0>
Shamrock phi_tilde_3d_prime(q) = <built-in method M4Shift16_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6d2100>
Absolute error f(q) = 3.3306690738754696212708950042725e-16
Absolute error df(q) = 3.552713678800501e-15
Absolute error ddf(q) = 2.2737367544323206e-12
Absolute error phi_tilde_3d(q) = 5.1146640037939727802974118892356e-11
Absolute error phi_tilde_3d_prime(q) = 7.6346989922726302352262500496407e-11
------------------------------------------
TGauss3 Kernel#
Generate c++ code for the kernel
1122 ret = kernel_to_shamrock(tgauss3)
1124 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel TGauss3
------------------------------------------
f(q) = Piecewise((exp(-q**2), q < 3), (0, True))
df(q) = Piecewise((-2*q*exp(-q**2), q < 3), (0, True))
ddf(q) = Piecewise((4*q**2*exp(-q**2) - 2*exp(-q**2), q < 3), (0, True))
phi_tilde_3d(q) = Piecewise((2*pi*exp(-9) - pi**(3/2)*erf(q)/q, q < 3), ((-pi**(3/2)*erf(3) + 6*pi*exp(-9))/q, True))
phi_tilde_3d_prime(q) = Piecewise((-2*pi*exp(-q**2)/q + pi**(3/2)*erf(q)/q**2, q < 3), (-(-pi**(3/2)*erf(3) + 6*pi*exp(-9))/q**2, True))
------------------------------------------
1126 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel TGauss3
------------------------------------------
template<class Tscal>
class KernelDefTGauss3 {
public:
inline static constexpr Tscal Rkern = 3; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 1/(sqrt(shambase::constants::pi<Tscal>)*erf(3));
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = 1.0/(-shambase::constants::pi<Tscal>*exp(-9) + shambase::constants::pi<Tscal>);
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (1.0/4.0)/(shambase::constants::pi<Tscal>*(-3.0/2.0*exp(-9) + (1.0/4.0)*sqrt(shambase::constants::pi<Tscal>)*erf(3)));
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
if (q < 3) {
return exp(-t1);
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
if (q < 3) {
return -2*q*exp(-t1);
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
if (q < 3) {
return 4*t1*exp(-t1) - 2*exp(-t1);
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
if (q < 3) {
return 2*shambase::constants::pi<Tscal>*exp(-9) - pow(shambase::constants::pi<Tscal>, 3.0/2.0)*erf(q)/q;
}
else
return (-pow(shambase::constants::pi<Tscal>, 3.0/2.0)*erf(3) + 6*shambase::constants::pi<Tscal>*exp(-9))/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
if (q < 3) {
return -2*shambase::constants::pi<Tscal>*exp(-t1)/q + pow(shambase::constants::pi<Tscal>, 3.0/2.0)*erf(q)/t1;
}
else
return -(-pow(shambase::constants::pi<Tscal>, 3.0/2.0)*erf(3) + 6*shambase::constants::pi<Tscal>*exp(-9))/t1;
}
};
------------------------------------------
Test the kernel
1129 test_result_tgauss3 = test_kernel(ret)
------------------------------------------
Testing kernel TGauss3 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.564202047051383, expected=0.564202047051383, delta=3.33066907387547E-16
Shamrock norm_2d = 0.318349173592935, expected=0.318349173592935, delta=1.11022302462516E-16
Shamrock norm_3d = 0.179666148218087, expected=0.179666148218087, delta=1.94289029309402E-16
Testing kernel radius:
Shamrock Rkern = 3.0, delta=0.0
Testing kernel functions:
Shamrock f(q) = <built-in method TGauss3_f of PyCapsule object at 0x7fb9ee6d2310>
Shamrock df(q) = <built-in method TGauss3_df of PyCapsule object at 0x7fb9ee6d2340>
Shamrock ddf(q) = <built-in method TGauss3_ddf of PyCapsule object at 0x7fb9ee6d2370>
Shamrock phi_tilde_3d(q) = <built-in method TGauss3_phi_tilde_3d of PyCapsule object at 0x7fb9ee6d23a0>
Shamrock phi_tilde_3d_prime(q) = <built-in method TGauss3_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6d23d0>
ZeroDivisionError: using limits towards 0
limit_f = 1.0 | shamrock_f[0] = 1.0
limit_df = 0.0 | shamrock_df[0] = -0.0
limit_ddf = -2.0 | shamrock_ddf[0] = -2.0
limit_phi_tilde_3d = -6.282409900511787 | shamrock_phi_tilde_3d[0] = -6.282409900511787
limit_phi_tilde_3d_prime = 0.0 | shamrock_phi_tilde_3d_prime[0] = 0.0
Absolute error f(q) = 5.5328869015178000863705218429141e-17
Absolute error df(q) = 1.0852585178826740881498128270435e-16
Absolute error ddf(q) = 2.6110952805685855689621535940408e-16
Absolute error phi_tilde_3d(q) = 1.5032018857438210402578376798422e-15
Absolute error phi_tilde_3d_prime(q) = 4.3266324410581327910712008064036e-14
------------------------------------------
TGauss5 Kernel#
Generate c++ code for the kernel
1138 ret = kernel_to_shamrock(tgauss5)
1140 print_kernel_info(ret)
------------------------------------------
Sympy expression for kernel TGauss5
------------------------------------------
f(q) = Piecewise((exp(-q**2), q < 5), (0, True))
df(q) = Piecewise((-2*q*exp(-q**2), q < 5), (0, True))
ddf(q) = Piecewise((4*q**2*exp(-q**2) - 2*exp(-q**2), q < 5), (0, True))
phi_tilde_3d(q) = Piecewise((2*pi*exp(-25) - pi**(3/2)*erf(q)/q, q < 5), ((-pi**(3/2)*erf(5) + 10*pi*exp(-25))/q, True))
phi_tilde_3d_prime(q) = Piecewise((-2*pi*exp(-q**2)/q + pi**(3/2)*erf(q)/q**2, q < 5), (-(-pi**(3/2)*erf(5) + 10*pi*exp(-25))/q**2, True))
------------------------------------------
1142 print_kernel_cpp_code(ret)
------------------------------------------
C++ generated code for kernel TGauss5
------------------------------------------
template<class Tscal>
class KernelDefTGauss5 {
public:
inline static constexpr Tscal Rkern = 5; ///< Compact support radius of the kernel
/// default hfact to be used for this kernel
inline static constexpr Tscal hfactd = 1.2;
/// 1D norm of the kernel
inline static constexpr Tscal norm_1d = 1/(sqrt(shambase::constants::pi<Tscal>)*erf(5));
/// 2D norm of the kernel
inline static constexpr Tscal norm_2d = 1.0/(-shambase::constants::pi<Tscal>*exp(-25) + shambase::constants::pi<Tscal>);
/// 3D norm of the kernel
inline static constexpr Tscal norm_3d = (1.0/4.0)/(shambase::constants::pi<Tscal>*(-5.0/2.0*exp(-25) + (1.0/4.0)*sqrt(shambase::constants::pi<Tscal>)*erf(5)));
inline static Tscal f(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
if (q < 5) {
return exp(-t1);
}
else
return 0;
}
inline static Tscal df(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
if (q < 5) {
return -2*q*exp(-t1);
}
else
return 0;
}
inline static Tscal ddf(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
if (q < 5) {
return 4*t1*exp(-t1) - 2*exp(-t1);
}
else
return 0;
}
inline static Tscal phi_tilde_3d(Tscal q) {
if (q < 5) {
return 2*shambase::constants::pi<Tscal>*exp(-25) - pow(shambase::constants::pi<Tscal>, 3.0/2.0)*erf(q)/q;
}
else
return (-pow(shambase::constants::pi<Tscal>, 3.0/2.0)*erf(5) + 10*shambase::constants::pi<Tscal>*exp(-25))/q;
}
inline static Tscal phi_tilde_3d_prime(Tscal q) {
Tscal t1 = sham::pow_constexpr<2>(q);
if (q < 5) {
return -2*shambase::constants::pi<Tscal>*exp(-t1)/q + pow(shambase::constants::pi<Tscal>, 3.0/2.0)*erf(q)/t1;
}
else
return -(-pow(shambase::constants::pi<Tscal>, 3.0/2.0)*erf(5) + 10*shambase::constants::pi<Tscal>*exp(-25))/t1;
}
};
------------------------------------------
Test the kernel
1146 test_result_tgauss5 = test_kernel(ret)
------------------------------------------
Testing kernel TGauss5 matching with Shamrock code
------------------------------------------
Testing norms:
Shamrock norm_1d = 0.564189583548624, expected=0.564189583548624, delta=3.33066907387547E-16
Shamrock norm_2d = 0.318309886188211, expected=0.318309886188211, delta=3.33066907387547E-16
Shamrock norm_3d = 0.179587122139514, expected=0.179587122139514, delta=1.11022302462516E-16
Testing kernel radius:
Shamrock Rkern = 5.0, delta=0.0
Testing kernel functions:
Shamrock f(q) = <built-in method TGauss5_f of PyCapsule object at 0x7fb9ee6d25e0>
Shamrock df(q) = <built-in method TGauss5_df of PyCapsule object at 0x7fb9ee6d2610>
Shamrock ddf(q) = <built-in method TGauss5_ddf of PyCapsule object at 0x7fb9ee6d2640>
Shamrock phi_tilde_3d(q) = <built-in method TGauss5_phi_tilde_3d of PyCapsule object at 0x7fb9ee6d2670>
Shamrock phi_tilde_3d_prime(q) = <built-in method TGauss5_phi_tilde_3d_prime of PyCapsule object at 0x7fb9ee6d26a0>
ZeroDivisionError: using limits towards 0
limit_f = 1.0 | shamrock_f[0] = 1.0
limit_df = 0.0 | shamrock_df[0] = -0.0
limit_ddf = -2.0 | shamrock_ddf[0] = -2.0
limit_phi_tilde_3d = -6.283185307092326 | shamrock_phi_tilde_3d[0] = -6.283185307092326
limit_phi_tilde_3d_prime = 0.0 | shamrock_phi_tilde_3d_prime[0] = 0.0
Absolute error f(q) = 5.5661915480185534057588247389932e-17
Absolute error df(q) = 1.0170157165450439613793689371898e-16
Absolute error ddf(q) = 2.6643412974770109075058790857962e-16
Absolute error phi_tilde_3d(q) = 1.6838964663548441663865583230091e-15
Absolute error phi_tilde_3d_prime(q) = 6.3308030348904237984452863555703e-14
------------------------------------------
Plot the kernels#
1153 fig_sz = (6.4, 12)
Plotting helper functions
1158 def create_kernel_plot_figure(fig_sz):
1159 """Create a figure with 5 subplots for kernel plotting"""
1160 plt.figure(figsize=fig_sz)
1161 ax_f = plt.subplot(5, 1, 1)
1162 ax_df = plt.subplot(5, 1, 2)
1163 ax_ddf = plt.subplot(5, 1, 3)
1164 ax_phi_tilde_3d = plt.subplot(5, 1, 4)
1165 ax_phi_tilde_3d_prime = plt.subplot(5, 1, 5)
1166 return ax_f, ax_df, ax_ddf, ax_phi_tilde_3d, ax_phi_tilde_3d_prime
1167
1168
1169 def plot_kernel_result(axes, test_result, kernel_label):
1170 """Plot a single kernel result on the given axes"""
1171 ax_f, ax_df, ax_ddf, ax_phi_tilde_3d, ax_phi_tilde_3d_prime = axes
1172 q_arr = test_result["q_arr"]
1173
1174 ax_f.plot(q_arr, test_result["shamrock_Cf"], label=f"C_3d f_{kernel_label}(q)")
1175 ax_df.plot(q_arr, test_result["shamrock_Cdf"], label=f"C_3d df_{kernel_label}(q)")
1176 ax_ddf.plot(q_arr, test_result["shamrock_Cddf"], label=f"C_3d ddf_{kernel_label}(q)")
1177 ax_phi_tilde_3d.plot(
1178 q_arr, test_result["shamrock_Cphi_tilde_3d"], label=f"C_3d phi_tilde_3d_{kernel_label}(q)"
1179 )
1180 ax_phi_tilde_3d_prime.plot(
1181 q_arr,
1182 test_result["shamrock_Cphi_tilde_3d_prime"],
1183 label=f"C_3d phi_tilde_3d_prime_{kernel_label}(q)",
1184 )
1185
1186
1187 def finalize_kernel_plot(axes):
1188 """Add titles, labels, and legends to the kernel plot"""
1189 ax_f, ax_df, ax_ddf, ax_phi_tilde_3d, ax_phi_tilde_3d_prime = axes
1190
1191 # Get current axis limits before adding reference line
1192 xlim = ax_phi_tilde_3d.get_xlim()
1193 ylim = ax_phi_tilde_3d.get_ylim()
1194
1195 ylim = (1.2 * ylim[0], ylim[1])
1196
1197 # Add -1/r reference line (only within current x range)
1198 q = np.linspace(max(1e-6, xlim[0]), xlim[1], 1000)
1199 one_over_r = -1 / q
1200 ax_phi_tilde_3d.plot(q, one_over_r, "--", color="grey", label="-1/r")
1201
1202 # Restore original limits (ignore reference line for autoscaling)
1203 ax_phi_tilde_3d.set_xlim(xlim)
1204 ax_phi_tilde_3d.set_ylim(ylim)
1205
1206 # Get current axis limits before adding reference line
1207 xlim = ax_phi_tilde_3d_prime.get_xlim()
1208 ylim = ax_phi_tilde_3d_prime.get_ylim()
1209
1210 ylim = (ylim[0], 1.5 * ylim[1])
1211
1212 # Add 1/r^2 reference line (only within current x range)
1213 one_over_r_squared = 1 / q**2
1214 ax_phi_tilde_3d_prime.plot(q, one_over_r_squared, "--", color="grey", label="1/r^2")
1215
1216 # Restore original limits (ignore reference line for autoscaling)
1217 ax_phi_tilde_3d_prime.set_xlim(xlim)
1218 ax_phi_tilde_3d_prime.set_ylim(ylim)
1219
1220 ax_f.set_title("C_3d f(q)", fontsize=10)
1221 ax_df.set_title("C_3d df(q)", fontsize=10)
1222 ax_ddf.set_title("C_3d ddf(q)", fontsize=10)
1223 ax_phi_tilde_3d.set_title("C_3d phi_tilde_3d(q)", fontsize=10)
1224 ax_phi_tilde_3d_prime.set_title("C_3d phi_tilde_3d_prime(q)", fontsize=10)
1225
1226 ax_f.set_xlabel("q")
1227 ax_df.set_xlabel("q")
1228 ax_ddf.set_xlabel("q")
1229 ax_phi_tilde_3d.set_xlabel("q")
1230 ax_phi_tilde_3d_prime.set_xlabel("q")
1231
1232 ax_f.legend(loc="right", fontsize=8)
1233 ax_df.legend(loc="right", fontsize=8)
1234 ax_ddf.legend(loc="right", fontsize=8)
1235 ax_phi_tilde_3d.legend(loc="right", fontsize=8)
1236 ax_phi_tilde_3d_prime.legend(loc="right", fontsize=8)
1237
1238 plt.tight_layout()
1239 plt.show()
Cubic splines kernels (M-series)
1245 axes = create_kernel_plot_figure(fig_sz)
1246 plot_kernel_result(axes, test_result_m4, "m4")
1247 plot_kernel_result(axes, test_result_m5, "m5")
1248 plot_kernel_result(axes, test_result_m6, "m6")
1249 plot_kernel_result(axes, test_result_m7, "m7")
1250 plot_kernel_result(axes, test_result_m8, "m8")
1251 plot_kernel_result(axes, test_result_m9, "m9")
1252 plot_kernel_result(axes, test_result_m10, "m10")
1253 finalize_kernel_plot(axes)
1254 plt.show()
Wendland kernels (C-series)
1259 axes = create_kernel_plot_figure(fig_sz)
1260 plot_kernel_result(axes, test_result_c2, "c2")
1261 plot_kernel_result(axes, test_result_c4, "c4")
1262 plot_kernel_result(axes, test_result_c6, "c6")
1263 finalize_kernel_plot(axes)
1264 plt.show()
Double hump kernels
1269 axes = create_kernel_plot_figure(fig_sz)
1270 plot_kernel_result(axes, test_result_m4dh, "m4dh")
1271 plot_kernel_result(axes, test_result_m4dh3, "m4dh3")
1272 plot_kernel_result(axes, test_result_m4dh5, "m4dh5")
1273 plot_kernel_result(axes, test_result_m4dh7, "m4dh7")
1274 finalize_kernel_plot(axes)
1275 plt.show()
M4Shift kernels
1281 axes = create_kernel_plot_figure(fig_sz)
1282 plot_kernel_result(axes, test_result_m4shift2, "m4shift2")
1283 plot_kernel_result(axes, test_result_m4shift4, "m4shift4")
1284 plot_kernel_result(axes, test_result_m4shift8, "m4shift8")
1285 plot_kernel_result(axes, test_result_m4shift16, "m4shift16")
1286 finalize_kernel_plot(axes)
1287 plt.show()
TGauss kernels
1292 axes = create_kernel_plot_figure(fig_sz)
1293 plot_kernel_result(axes, test_result_tgauss3, "tgauss3")
1294 plot_kernel_result(axes, test_result_tgauss5, "tgauss5")
1295 finalize_kernel_plot(axes)
1296 plt.show()
Total running time of the script: (1 minutes 43.851 seconds)
Estimated memory usage: 264 MB