diff --git a/.github/workflows/test-struphy.yml b/.github/workflows/test-struphy.yml index 807fd0220..3fef4e200 100644 --- a/.github/workflows/test-struphy.yml +++ b/.github/workflows/test-struphy.yml @@ -58,6 +58,7 @@ jobs: echo "Psydac location for this branch" pip show psydac pip uninstall psydac -y + git checkout 108-psydac-change-renaming-of-globalpojector python -m pip install ".[phys]" --no-cache-dir echo "Psydac location after installing struphy" pip show psydac diff --git a/docs/source/modules/feec.multipatch.rst b/docs/source/modules/feec.multipatch.rst index 64bf49a12..85bf08af1 100644 --- a/docs/source/modules/feec.multipatch.rst +++ b/docs/source/modules/feec.multipatch.rst @@ -7,11 +7,6 @@ feec.multipatch :toctree: STUBDIR :template: autosummary/module.rst - multipatch.api - multipatch.fem_linear_operators multipatch.multipatch_domain_utilities - multipatch.non_matching_operators - multipatch.operators - multipatch.plotting_utilities multipatch.utilities multipatch.utils_conga_2d diff --git a/docs/source/modules/feec.rst b/docs/source/modules/feec.rst index 5d424f396..f75dfb8cc 100644 --- a/docs/source/modules/feec.rst +++ b/docs/source/modules/feec.rst @@ -7,8 +7,10 @@ feec :toctree: STUBDIR :template: autosummary/module.rst + feec.conforming_projectors feec.derivatives - feec.global_projectors + feec.global_geometric_projectors + feec.hodge feec.pull_push feec.pushforward diff --git a/docs/source/modules/linalg.rst b/docs/source/modules/linalg.rst index 7ccbbf491..a625296ef 100644 --- a/docs/source/modules/linalg.rst +++ b/docs/source/modules/linalg.rst @@ -14,6 +14,7 @@ linalg linalg.kernels linalg.kron linalg.solvers + linalg.sparse linalg.stencil linalg.topetsc linalg.utilities diff --git a/performance/compare_3d_matrix_assembly_speed.py b/performance/compare_3d_matrix_assembly_speed.py new file mode 100644 index 000000000..0bc8678a5 --- /dev/null +++ b/performance/compare_3d_matrix_assembly_speed.py @@ -0,0 +1,643 @@ +import os +import shutil +import time +from pathlib import Path +from datetime import datetime + +from mpi4py import MPI +import numpy as np +#import matplotlib.pyplot as plt +from sympy import sin + +from sympde.calculus import dot, cross, grad, curl +from sympde.expr import BilinearForm, integral +from sympde.topology import element_of, elements_of, Cube, Mapping, ScalarFunctionSpace, Domain, Derham + +from psydac.api.discretization import discretize +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL +from psydac.cad.geometry import Geometry +from psydac.fem.basic import FemField +from psydac.mapping.discrete import SplineMapping + +datetime_md = datetime.today().strftime('%Y-%m-%d %H:%M:%S') +datetime_file = datetime.today().strftime('%Y-%m-%d_%H:%M:%S') + +""" +Executing this file will add new tables to the matrix_assembly_speed_log.md file in which we keep track of assembly and discretization speed. +Estimated runtime: ~4 min. + +""" + +comm = MPI.COMM_WORLD +backend = PSYDAC_BACKEND_GPYCCEL +mpi_rank = comm.rank + +if mpi_rank == 0: + print('Expected runtime: 4 min.') + print() + +class SquareTorus(Mapping): + + _expressions = {'x': 'x1 * cos(x2)', + 'y': 'x1 * sin(x2)', + 'z': 'x3'} + + _ldim = 3 + _pdim = 3 + +def make_square_torus_geometry_3d(ncells, degree, comm=None): + + if comm is not None: + mpi_rank = comm.Get_rank() + else: + mpi_rank = 0 + + if (ncells[0] == ncells[1]) and (ncells[0] == ncells[2]): + nc = f'{ncells[0]}' + else: + nc = f'{ncells[0]}_{ncells[1]}_{ncells[2]}' + + if (degree[0] == degree[1]) and (degree[0] == degree[2]): + de = f'{degree[0]}' + else: + de = f'{degree[0]}_{degree[1]}_{degree[2]}' + + name = f'st_3d_nc_{nc}_d_{de}' + + r = 0.5 + R = 1. + logical_domain = Cube('C', bounds1=(r, R), bounds2=(0, 2*np.pi), bounds3=(0, 1)) + domain_h = discretize(logical_domain, ncells=ncells, comm=comm) + + V = ScalarFunctionSpace('V', logical_domain) + V_h = discretize(V, domain_h, degree=degree) + + mapping = SquareTorus('S') + map_discrete = SplineMapping.from_mapping(V_h, mapping.get_callable_mapping()) + + geometry = Geometry.from_discrete_mapping(map_discrete, comm=comm) + + if mpi_rank == 0: + if not os.path.isdir('geometry'): + os.makedirs('geometry') + os.makedirs('geometry/files') + else: + if not os.path.isdir('geometry/files'): + os.makedirs('geometry/files') + + geometry.export(f'geometry/files/{name}.h5') + + return f'geometry/files/{name}.h5' + +# ---------- 1 ---------- +#t0_1_glob = time.time() +# +#ncells = [32, 32, 32] +#degree_list = [[2, 2, 2], [3, 3, 3], [4, 4, 4], [5, 5, 5]] +#degree_list2 = [[2, 2, 2], [3, 3, 3], [4, 4, 4]] +#periodic = [False, False, False] +# +##mapping = HalfHollowTorusMapping3D('M', R=2, r=1) +#mapping = SquareTorus('S') +##logical_domain = Cube('C', bounds1=(0,1), bounds2=(0,1), bounds3=(0,1)) +#r = 0.5 +#R = 1. +#logical_domain = Cube('C', bounds1=(r, R), bounds2=(0, 2*np.pi), bounds3=(0, 1)) +# +#domain = mapping(logical_domain) +#derham = Derham(domain) +#plot_domain(domain, draw=True, isolines=True) +# +#ass_time_H1_old = [[] for _ in degree_list] +#ass_time_H1_new = [[] for _ in degree_list] +# +#ass_time_Hcurl_old = [[] for _ in degree_list2] +#ass_time_Hcurl_new = [[] for _ in degree_list2] +# +#for i, degree in enumerate(degree_list): +# +# domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm) +# derham_h = discretize(derham, domain_h, degree=degree) +# +# V0 = derham.V0 +# V0h = derham_h.V0 +# +# u, v = elements_of(V0, names='u, v') +# +# a = BilinearForm((u, v), integral(domain, u*v)) +# +# a_h = discretize(a, domain_h, (V0h, V0h), backend=backend, sum_factorization=False) +# +# t0_ao = time.time() +# M_o = a_h.assemble() +# t1_ao = time.time() +# +# a_h = discretize(a, domain_h, (V0h, V0h), backend=backend) +# +# t0_an = time.time() +# M_n = a_h.assemble() +# t1_an = time.time() +# +# ass_time_H1_old[i].append(t1_ao - t0_ao) +# ass_time_H1_new[i].append(t1_an - t0_an) +# +# if degree != [5, 5, 5]: +# +# V1 = derham.V1 +# V1h = derham_h.V1 +# +# u, v = elements_of(V1, names='u, v') +# +# a = BilinearForm((u, v), integral(domain, dot(u, v))) +# +# a_h = discretize(a, domain_h, (V1h, V1h), backend=backend, sum_factorization=False) +# +# t0_ao = time.time() +# M_o = a_h.assemble() +# t1_ao = time.time() +# +# a_h = discretize(a, domain_h, (V1h, V1h), backend=backend) +# +# t0_an = time.time() +# M_n = a_h.assemble() +# t1_an = time.time() +# +# ass_time_Hcurl_old[i].append(t1_ao - t0_ao) +# ass_time_Hcurl_new[i].append(t1_an - t0_an) +# +#d = [degree[0] for degree in degree_list] +#d2 = [degree[0] for degree in degree_list2] +# +#ass_time_H1_old_d = [ass_time_H1_old[i][0] for i, _ in enumerate(degree_list)] +#ass_time_H1_new_d = [ass_time_H1_new[i][0] for i, _ in enumerate(degree_list)] +# +#ass_time_Hcurl_old_d = [ass_time_Hcurl_old[i][0] for i, _ in enumerate(degree_list2)] +#ass_time_Hcurl_new_d = [ass_time_Hcurl_new[i][0] for i, _ in enumerate(degree_list2)] +# +#if mpi_rank == 0: +# plt.plot(d, ass_time_H1_old_d, '--.', label=f'old Algorithm') +# plt.plot(d, ass_time_H1_new_d, '--.', label=f'new Algorithm') +# plt.title(r'Assembly times for the $H^1(\Omega)$ mass matrix') +# plt.legend() +# plt.yscale('log') +# plt.ylabel('Wallclock Time [s]') +# plt.xlabel('d - [d, d, d] Bspline degrees') +# plt.xticks(d) +# #plt.savefig(f'figures/H1_{datetime_file}.png') +# plt.show() +# plt.clf() +# +# plt.plot(d2, ass_time_Hcurl_old_d, '--.', label=f'old Algorithm') +# plt.plot(d2, ass_time_Hcurl_new_d, '--.', label=f'new Algorithm') +# plt.title(r'Assembly times for the $H(curl;\Omega)$ mass matrix') +# plt.legend() +# plt.yscale('log') +# plt.ylabel('Wallclock Time [s]') +# plt.xlabel('d - [d, d, d] Bspline degrees') +# plt.xticks(d2) +# #plt.savefig(f'figures/Hcurl_{datetime_file}.png') +# plt.show() +# +#t1_1_glob = time.time() +#if mpi_rank == 0: +# print(f'Part 1 out of 3 done after {(t1_1_glob-t0_1_glob)/60:.2g}min') +# ----------------------- + +# ---------- 2 ---------- +t0_2_glob = time.time() + +ncells = [32, 32, 32] +degree = [3, 3, 3] +periodic = [False, False, False] + +r = 0.5 +R = 1. +logical_domain = Cube('C', bounds1=(r, R), bounds2=(0, 2*np.pi), bounds3=(0, 1)) +logical_derham = Derham(logical_domain) + +logical_domain_h = discretize(logical_domain, ncells=ncells, periodic=periodic, comm=comm) +logical_derham_h = discretize(logical_derham, logical_domain_h, degree=degree) + +filename = make_square_torus_geometry_3d(ncells, degree, comm=comm) + +bspline_domain = Domain.from_file(filename) +bspline_derham = Derham(bspline_domain) + +bspline_domain_h = discretize(bspline_domain, filename=filename, comm=comm) +bspline_derham_h = discretize(bspline_derham, bspline_domain_h, degree=bspline_domain.mapping.get_callable_mapping().space.degree) + +mapping = SquareTorus('S') + +analytical_domain = mapping(logical_domain) +analytical_derham = Derham(analytical_domain) + +analytical_domain_h = discretize(analytical_domain, ncells=ncells, periodic=periodic, comm=comm) +analytical_derham_h = discretize(analytical_derham, analytical_domain_h, degree=degree) + +ax, ay, az = analytical_domain.coordinates +agamma = ax*ay*az + sin(ax*ay+az)**2 + +# 2.1 +dom = logical_domain +domh = logical_domain_h +V = logical_derham.V1 +Vh = logical_derham_h.V1 + +F = element_of(V, name='F') +f = Vh.coeff_space.zeros() +f[0]._data = np.ones(f[0]._data.shape) +f[1]._data = np.ones(f[1]._data.shape) +f[2]._data = np.ones(f[2]._data.shape) +f_field = FemField(Vh, f) + +u, v = elements_of(V, names='u, v') +a = BilinearForm((u, v), integral(dom, dot(cross(F, u), cross(F, v)))) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend, sum_factorization=False) +t1 = time.time() +old_disc_21 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble(F=f_field) +t1 = time.time() +old_ass_21 = round(t1-t0, 3) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend) +t1 = time.time() +new_disc_21 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble(F=f_field) +t1 = time.time() +new_ass_21 = round(t1-t0, 3) + +# 2.2 +dom = analytical_domain +domh = analytical_domain_h +V = analytical_derham.V2 +Vh = analytical_derham_h.V2 + +u, v = elements_of(V, names='u, v') +a = BilinearForm((u, v), integral(dom, dot(u, v)*agamma)) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend, sum_factorization=False) +t1 = time.time() +old_disc_22 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +old_ass_22 = round(t1-t0, 3) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend) +t1 = time.time() +new_disc_22 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +new_ass_22 = round(t1-t0, 3) + +# 2.3 +dom = bspline_domain +domh = bspline_domain_h +V = bspline_derham.V1 +Vh = bspline_derham_h.V1 + +u, v = elements_of(V, names='u, v') +a = BilinearForm((u, v), integral(dom, dot(curl(u), curl(v)))) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend, sum_factorization=False) +t1 = time.time() +old_disc_23 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +old_ass_23 = round(t1-t0, 3) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend) +t1 = time.time() +new_disc_23 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +new_ass_23 = round(t1-t0, 3) + +t1_2_glob = time.time() +if mpi_rank == 0: + print(f'Part 2 out of 3 done after {(t1_2_glob-t0_2_glob)/60:.2g}min') +# ----------------------- + +# ---------- 3 ---------- +t0_3_glob = time.time() + +ncells = [16, 8, 32] +degree = [2, 4, 3] +periodic = [False, True, False] + +r = 0.5 +R = 1. +logical_domain = Cube('C', bounds1=(r, R), bounds2=(0, 2*np.pi), bounds3=(0, 1)) +logical_derham = Derham(logical_domain) + +logical_domain_h = discretize(logical_domain, ncells=ncells, periodic=periodic, comm=comm) +logical_derham_h = discretize(logical_derham, logical_domain_h, degree=degree) + +filename = make_square_torus_geometry_3d(ncells, degree, comm=comm) + +bspline_domain = Domain.from_file(filename) +bspline_derham = Derham(bspline_domain) + +bspline_domain_h = discretize(bspline_domain, filename=filename, comm=comm) +bspline_derham_h = discretize(bspline_derham, bspline_domain_h, degree=bspline_domain.mapping.get_callable_mapping().space.degree) + +mapping = SquareTorus('S') + +analytical_domain = mapping(logical_domain) +analytical_derham = Derham(analytical_domain) + +analytical_domain_h = discretize(analytical_domain, ncells=ncells, periodic=periodic, comm=comm) +analytical_derham_h = discretize(analytical_derham, analytical_domain_h, degree=degree) + +ax, ay, az = analytical_domain.coordinates +agamma = ax*ay*az + sin(ax*ay+az)**2 + +# 3.1.1 +dom = logical_domain +domh = logical_domain_h +V = logical_derham.V0 +Vh = logical_derham_h.V0 + +u, v = elements_of(V, names='u, v') +a = BilinearForm((u, v), integral(dom, dot(grad(u), grad(v)))) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend, sum_factorization=False) +t1 = time.time() +old_disc_311 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +old_ass_311 = round(t1-t0, 3) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend) +t1 = time.time() +new_disc_311 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +new_ass_311 = round(t1-t0, 3) + +# 3.1.2 +dom = analytical_domain +domh = analytical_domain_h +V = analytical_derham.V0 +Vh = analytical_derham_h.V0 + +u, v = elements_of(V, names='u, v') +a = BilinearForm((u, v), integral(dom, dot(grad(u), grad(v)))) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend, sum_factorization=False) +t1 = time.time() +old_disc_312 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +old_ass_312 = round(t1-t0, 3) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend) +t1 = time.time() +new_disc_312 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +new_ass_312 = round(t1-t0, 3) + +# 3.1.3 +dom = bspline_domain +domh = bspline_domain_h +V = bspline_derham.V0 +Vh = bspline_derham_h.V0 + +u, v = elements_of(V, names='u, v') +a = BilinearForm((u, v), integral(dom, dot(grad(u), grad(v)))) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend, sum_factorization=False) +t1 = time.time() +old_disc_313 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +old_ass_313 = round(t1-t0, 3) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend) +t1 = time.time() +new_disc_313 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +new_ass_313 = round(t1-t0, 3) + +# 3.2 +dom = analytical_domain +domh = analytical_domain_h +V = ScalarFunctionSpace('V', analytical_domain) +Vh = discretize(V, analytical_domain_h, degree=degree) + +u, v = elements_of(V, names='u, v') +a = BilinearForm((u, v), integral(dom, dot(grad(u), grad(v)))) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend, sum_factorization=False) +t1 = time.time() +old_disc_32 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +old_ass_32 = round(t1-t0, 3) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend) +t1 = time.time() +new_disc_32 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +new_ass_32 = round(t1-t0, 3) + +# 3.3 +dom = analytical_domain +domh = analytical_domain_h +V = ScalarFunctionSpace('V', analytical_domain) +Vh = discretize(V, analytical_domain_h, degree=degree) +W = analytical_derham.V0 +Wh = analytical_derham_h.V0 + +u = element_of(V, name='u') +v = element_of(W, name='v') +a = BilinearForm((u, v), integral(dom, dot(grad(u), grad(v)))) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend, sum_factorization=False) +t1 = time.time() +old_disc_33 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +old_ass_33 = round(t1-t0, 3) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend) +t1 = time.time() +new_disc_33 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +new_ass_33 = round(t1-t0, 3) + +# ------------------ +dom = analytical_domain +domh = analytical_domain_h +V = analytical_derham.V0 +Vh = analytical_derham_h.V0 + +u, v = elements_of(V, names='u, v') +# ------------------ + +# 3.4 +a = BilinearForm((u, v), integral(dom, dot(grad(u), grad(v)) * agamma)) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend, sum_factorization=False) +t1 = time.time() +old_disc_34 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +old_ass_34 = round(t1-t0, 3) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend) +t1 = time.time() +new_disc_34 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble() +t1 = time.time() +new_ass_34 = round(t1-t0, 3) + +# 3.5 +F = element_of(V, name='F') +f = Vh.coeff_space.zeros() +f._data = np.ones(f._data.shape) +f_field = FemField(Vh, f) + +a = BilinearForm((u, v), integral(dom, dot(grad(u), grad(v)) * F)) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend, sum_factorization=False) +t1 = time.time() +old_disc_35 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble(F=f_field) +t1 = time.time() +old_ass_35 = round(t1-t0, 3) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend) +t1 = time.time() +new_disc_35 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble(F=f_field) +t1 = time.time() +new_ass_35 = round(t1-t0, 3) + +# 3.6 +mult = [1, 3, 2] +analytical_derham_h = discretize(analytical_derham, analytical_domain_h, degree=degree, multiplicity=mult) +Vh = analytical_derham_h.V0 + +a = BilinearForm((u, v), integral(dom, dot(grad(u), grad(v)) * F)) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend, sum_factorization=False) +t1 = time.time() +old_disc_36 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble(F=f_field) +t1 = time.time() +old_ass_36 = round(t1-t0, 3) + +t0 = time.time() +ah = discretize(a, domh, (Vh, Vh), backend=backend) +t1 = time.time() +new_disc_36 = round(t1-t0, 3) +t0 = time.time() +M = ah.assemble(F=f_field) +t1 = time.time() +new_ass_36 = round(t1-t0, 3) + +t1_3_glob = time.time() +if mpi_rank == 0: + print(f'Part 3 out of 3 done after {(t1_3_glob-t0_3_glob)/60:.2g}min') +# ----------------------- + +template = '''| Test case | old assembly | new assembly | old discretization | new discretization | +| --- | --- | --- | --- | --- | +| 2.1 | {old_ass_21} | {new_ass_21} | {old_disc_21} | {new_disc_21} | +| 2.2 | {old_ass_22} | {new_ass_22} | {old_disc_22} | {new_disc_22} | +| 2.3 | {old_ass_23} | {new_ass_23} | {old_disc_23} | {new_disc_23} | +| 3.1.1 | {old_ass_311} | {new_ass_311} | {old_disc_311} | {new_disc_311} | +| 3.1.2 | {old_ass_312} | {new_ass_312} | {old_disc_312} | {new_disc_312} | +| 3.1.3 | {old_ass_313} | {new_ass_313} | {old_disc_313} | {new_disc_313} | +| 3.2 | {old_ass_32} | {new_ass_32} | {old_disc_32} | {new_disc_32} | +| 3.3 | {old_ass_33} | {new_ass_33} | {old_disc_33} | {new_disc_33} | +| 3.4 | {old_ass_34} | {new_ass_34} | {old_disc_34} | {new_disc_34} | +| 3.5 | {old_ass_35} | {new_ass_35} | {old_disc_35} | {new_disc_35} | +| 3.6 | {old_ass_36} | {new_ass_36} | {old_disc_36} | {new_disc_36} |''' + +txt = '' +txt += f'{datetime_md}\n' +txt += f'----------\n\n' +#txt += f'![](tests/figures/H1_{datetime_file}.png)\n' +#txt += f'![](tests/figures/Hcurl_{datetime_file}.png)\n\n' +txt += template.format(old_ass_21=old_ass_21, new_ass_21=new_ass_21, old_disc_21=old_disc_21, new_disc_21=new_disc_21, + old_ass_22=old_ass_22, new_ass_22=new_ass_22, old_disc_22=old_disc_22, new_disc_22=new_disc_22, + old_ass_23=old_ass_23, new_ass_23=new_ass_23, old_disc_23=old_disc_23, new_disc_23=new_disc_23, + old_ass_311=old_ass_311, new_ass_311=new_ass_311, old_disc_311=old_disc_311, new_disc_311=new_disc_311, + old_ass_312=old_ass_312, new_ass_312=new_ass_312, old_disc_312=old_disc_312, new_disc_312=new_disc_312, + old_ass_313=old_ass_313, new_ass_313=new_ass_313, old_disc_313=old_disc_313, new_disc_313=new_disc_313, + old_ass_32=old_ass_32, new_ass_32=new_ass_32, old_disc_32=old_disc_32, new_disc_32=new_disc_32, + old_ass_33=old_ass_33, new_ass_33=new_ass_33, old_disc_33=old_disc_33, new_disc_33=new_disc_33, + old_ass_34=old_ass_34, new_ass_34=new_ass_34, old_disc_34=old_disc_34, new_disc_34=new_disc_34, + old_ass_35=old_ass_35, new_ass_35=new_ass_35, old_disc_35=old_disc_35, new_disc_35=new_disc_35, + old_ass_36=old_ass_36, new_ass_36=new_ass_36, old_disc_36=old_disc_36, new_disc_36=new_disc_36) +txt += '\n\n' + +if mpi_rank == 0: + # Write performance table to MarkDown file + with open('matrix_assembly_speed_log.md', 'a') as f: + f.write(txt) + + # Remove temporary folders + base_dir = Path(__file__).parent + dirs_to_remove = [ + "geometry", + "__psydac__", + "__pycache__", + "__epyccel__", + "__gpyccel__" + ] + for d in dirs_to_remove: + path = base_dir / d + if path.exists(): + try: + shutil.rmtree(path) + except Exception as e: + print(f"Failed to remove {path}: {e}") diff --git a/performance/matrix_assembly_speed_log.md b/performance/matrix_assembly_speed_log.md new file mode 100644 index 000000000..996975d26 --- /dev/null +++ b/performance/matrix_assembly_speed_log.md @@ -0,0 +1,86 @@ +New Matrix Assembly for Psydac +------------------------------ + +Here we keep track on the performance of the new assembly algorithm (sum factorization). +We measure both the discretization time of `BilinearForms` as well as the matrix assembly time of a `DiscreteBilinearForm` +and compare it to the old algorithm. Executing `compare_3d_matrix_assembly_speed.py` will add new data to this file. +This allows us to detect whether any future changes have a positive or negative impact. + +(Of course, runtime depends on the machine used to execute this file. Hence, an even decrease in runtime +is no reason to celebrate, and an even increase in runtime no reason to worry.) + +Test cases +---------- + +1. Scaling Analysis \ (**REMOVED**) + We assemble the H1 and Hcurl mass matrices ($\Omega$ is a half hollow torus (analytical mapping), ncells = [32, 32, 32], periodic=[False, False, False]) for varying degrees + and report the assembly times in Figures. + +2. Three specific test cases \ + We report assembly and discretization time for the following test cases \ + 2.1 "Q" + - $(u, v; F)\mapsto\int_{\Omega}(F\times u)\circ(F\times v)$ + - $u, v, F\in H(curl;\Omega)$, $F$ a "free field" + - no mapping ($\Omega=(0,1)^3$) + - ncells = [32, 32, 32], degree = [3, 3, 3], periodic = [False, False, False] + + 2.2 "weighted $Hdiv$ mass matrix" + - $(u, v; \gamma)\mapsto\int_{\Omega}u\circ v\ \gamma$ + - $u, v\in H(div;\Omega)$, $\gamma$ an analytic scalar weight function + - analytical mapping ($\Omega$ a half hollow torus) + - ncells = [32, 32, 32], degree = [3, 3, 3], periodic = [False, False, False] + + 2.3 "curl curl" + - $(u, v)\mapsto\int_{\Omega}(\nabla\times u)\circ(\nabla\times v)$ + - $u, v\in H(curl;\Omega)$ + - Bspline mapping ($\Omega$ a half hollow torus) + - ncells = [32, 32, 32], degree = [3, 3, 3], periodic = [False, False, False] + +3. More variations! + - $(u, v)\mapsto\int_{\Omega}\nabla u\circ\nabla v$ + - ncells = [16, 8, 32], degree = [2, 4, 3], periodic = [False, True, False] + - 3.1.1 + - `u, v = elements_of(derham.V0, names='u, v')` + - no mapping + - 3.1.2 + - `u, v = elements_of(derham.V0, names='u, v')` + - analytical mapping + - 3.1.3 + - `u, v = elements_of(derham.V0, names='u, v')` + - Bspline mapping + - 3.2 + - `u, v = elements_of(ScalarFunctionSpace('Vs', domain), names='u, v')` + - analytical mapping + - 3.3 + - `u = element_of(derham.V0, name='u')` + - `v = element_of(ScalarFunctionSpace('Vs', domain), name='v')` + - analytical mapping + - `u, v = elements_of(derham.V0, names='u, v')` + - analytical mapping + - 3.4 + - additional analytical weight function $\gamma$ + - 3.5 + - additional scalar free field $F\in H^1(\Omega)$ + - 3.6 + - multiplicity vector [1, 3, 2] + +Data +---- + +2025-09-09 17:28:19 (added by Julian O. - ThinkPad T14 on performance mode) +---------- + +| Test case | old assembly | new assembly | old discretization | new discretization | +| --- | --- | --- | --- | --- | +| 2.1 | 21.358 | 1.676 | 16.568 | 16.589 | +| 2.2 | 2.348 | 0.322 | 4.533 | 3.38 | +| 2.3 | 36.858 | 3.171 | 16.749 | 19.413 | +| 3.1.1 | 0.236 | 0.044 | 1.462 | 1.482 | +| 3.1.2 | 0.54 | 0.051 | 1.722 | 1.52 | +| 3.1.3 | 0.802 | 0.169 | 3.579 | 2.414 | +| 3.2 | 0.517 | 0.045 | 1.027 | 1.425 | +| 3.3 | 0.525 | 0.06 | 1.023 | 1.527 | +| 3.4 | 0.524 | 0.044 | 1.085 | 1.547 | +| 3.5 | 0.577 | 0.057 | 1.51 | 1.784 | +| 3.6 | 0.598 | 0.209 | 2.411 | 1.848 | + diff --git a/psydac/api/fem_bilinear_form.py b/psydac/api/fem_bilinear_form.py new file mode 100644 index 000000000..8af05df83 --- /dev/null +++ b/psydac/api/fem_bilinear_form.py @@ -0,0 +1,2226 @@ +import sys +import os +import importlib + +import numpy as np + +from sympy import ImmutableDenseMatrix, Matrix, Symbol, sympify +from sympy.tensor.indexed import Indexed, IndexedBase +from sympy.simplify import cse_main + +from pyccel import epyccel + +from sympde.topology.basic import Boundary, Interface +from sympde.topology.mapping import Mapping, SymbolicExpr +from sympde.topology.space import ScalarFunction, VectorFunction, IndexedVectorFunction +from sympde.topology.derivatives import get_atom_logical_derivatives +from sympde.topology.derivatives import _logical_partial_derivatives +from sympde.topology.derivatives import get_index_logical_derivatives +from sympde.topology.derivatives import get_max_logical_partial_derivatives # NOTE [YG 31.07.2025]: Maybe use the one in ast.utilities +from sympde.expr.expr import BilinearForm +from sympde.expr.evaluation import KernelExpression, TerminalExpr +from sympde.calculus.core import PlusInterfaceOperator + +from psydac.cad.geometry import Geometry +from psydac.mapping.discrete import SplineMapping, NurbsMapping +from psydac.fem.basic import FemSpace, FemField +from psydac.fem.vector import VectorFemSpace +from psydac.linalg.stencil import StencilMatrix +from psydac.linalg.block import BlockVectorSpace, BlockLinearOperator +from psydac.api.grid import QuadratureGrid, BasisValues +from psydac.api.settings import PSYDAC_BACKENDS +from psydac.api.utilities import flatten, random_string +from psydac.api.fem_common import ( + compute_imports, + compute_max_nderiv, + compute_free_arguments, + construct_test_space_arguments, + construct_trial_space_arguments, + construct_quad_grids_arguments, + reset_arrays, + do_nothing, + extract_stencil_mats, +) + +# TODO [YG 01.08.2025]: Avoid importing anything from psydac.pyccel +from psydac.pyccel.ast.core import _atomic, Assign + +__all__ = ('DiscreteBilinearForm',) + +NoneType = type(None) + +#============================================================================== +class DiscreteBilinearForm: + """ + Discrete bilinear form ready to be assembled into a matrix. + + This class represents the concept of a discrete bilinear form in Psydac. + Instances of this class generate an appropriate matrix assembly kernel, + allocate the matrix if not provided, and prepare a list of arguments for + the kernel. + + An implementation of the sum factorization algorithm is used to assemble + the matrix. + + Parameters + ---------- + + expr : sympde.expr.expr.BilinearForm + The symbolic bilinear form. + + kernel_expr : list or tuple of sympde.expr.evaluation.KernelExpression + The atomic representation of the bilinear form. + + domain_h : psydac.cad.geometry.Geometry + The discretized domain. + + spaces : list of psydac.fem.basic.FemSpace + The discrete trial and test spaces. + + nquads : list or tuple of int + The number of quadrature points used in the assembly kernel along each + direction. + + matrix : psydac.linalg.stencil.StencilMatrix or psydac.linalg.block.BlockLinearOperator, optional + The matrix that we assemble into. If not provided, a new matrix is + created with the appropriate domain and codomain (default: None). + + update_ghost_regions : bool, default=True + Accumulate the contributions of the neighbouring processes. + + backend : dict, optional + The backend used to accelerate the computing kernels. + The backend dictionaries are defined in the file psydac/api/settings.py + + assembly_backend : dict, optional + The backend used to accelerate the assembly kernel. + The backend dictionaries are defined in the file psydac/api/settings.py + + linalg_backend : dict, optional + The backend used to accelerate the computing kernels of the linear operator. + The backend dictionaries are defined in the file psydac/api/settings.py + + symbolic_mapping : sympde.topology.mapping.Mapping, optional + The symbolic mapping which defines the physical domain of the bilinear form. + + See Also + -------- + DiscreteLinearForm + DiscreteFunctional + DiscreteSumForm + + """ + def __init__(self, expr, kernel_expr, domain_h, spaces, *, nquads, + matrix=None, update_ghost_regions=True, backend=None, + linalg_backend=None, assembly_backend=None, + symbolic_mapping=None): + + #... Sanity checks + assert isinstance(expr, BilinearForm) + assert isinstance(domain_h, Geometry) + for space in spaces: + assert isinstance(space, FemSpace) + for nquad in nquads: + assert isinstance(nquad, int) + assert nquad > 0 + assert isinstance(matrix, (NoneType, StencilMatrix, BlockLinearOperator)) + assert isinstance(update_ghost_regions, bool) + assert isinstance( backend, (NoneType, dict)) + assert isinstance( linalg_backend, (NoneType, dict)) + assert isinstance(assembly_backend, (NoneType, dict)) + assert isinstance(symbolic_mapping, (NoneType, Mapping)) + #... + + if isinstance(kernel_expr, (tuple, list)): + if len(kernel_expr) == 1: + kernel_expr = kernel_expr[0] + else: + raise ValueError('> Expecting only one kernel') + assert isinstance(kernel_expr, KernelExpression) + + self._kernel_expr = kernel_expr + self._expr = expr + self._target = kernel_expr.target + self._domain = domain_h.domain + self._spaces = spaces + self._matrix = matrix + + domain = self.domain + target = self.target + + # ... + if len(domain) > 1: + i, j = self.get_space_indices_from_target(domain, target) + test_space = self.spaces[1].spaces[i] + trial_space = self.spaces[0].spaces[j] + if isinstance(target, Interface): + m,_ = self.get_space_indices_from_target(domain, target.minus) + p,_ = self.get_space_indices_from_target(domain, target.plus) + mapping_m = list(domain_h.mappings.values())[m] + mapping_p = list(domain_h.mappings.values())[p] + mapping = (mapping_m, mapping_p) if mapping_m else None + else: + mapping = list(domain_h.mappings.values())[i] + else: + trial_space = self.spaces[0] + test_space = self.spaces[1] + mapping = list(domain_h.mappings.values())[0] + + self._mapping = mapping + + is_rational_mapping = False + mapping_space = None + if (mapping is not None) and not isinstance(target, Interface): + is_rational_mapping = isinstance(mapping, NurbsMapping) + mapping_space = mapping.space + elif (mapping is not None) and isinstance(target, Interface): + is_rational_mapping = (isinstance(mapping[0], NurbsMapping), isinstance(mapping[1], NurbsMapping)) + mapping_space = (mapping[0].space, mapping[1].space) + + self._is_rational_mapping = is_rational_mapping + # ... + + if isinstance(test_space.coeff_space, BlockVectorSpace): + coeff_space = test_space.coeff_space.spaces[0] + else: + coeff_space = test_space.coeff_space + + self._coeff_space = coeff_space + self._num_threads = 1 + if coeff_space.parallel and coeff_space.cart.num_threads > 1: + self._num_threads = coeff_space.cart.num_threads + + self._update_ghost_regions = update_ghost_regions + + # In case of multiple patches, if the communicator is MPI_COMM_NULL, we do not generate the assembly code + # because the patch is not owned by the MPI rank. + if coeff_space.parallel and coeff_space.cart.is_comm_null: + self._free_args = () + self._func = do_nothing + self._args = () + self._threads_args = () + self._global_matrices = () + self._update_ghost_regions = False + return + + # ... + test_ext = None + trial_ext = None + if isinstance(target, Boundary): + axis = target.axis + test_ext = target.ext + trial_ext = target.ext + elif isinstance(target, Interface): + # this part treats the cases of: + # integral(v_minus * u_plus) + # integral(v_plus * u_minus) + # the other cases, integral(v_minus * u_minus) and integral(v_plus * u_plus) + # are converted to boundary integrals by Sympde + axis = target.axis + test = self.kernel_expr.test + trial = self.kernel_expr.trial + test_target = target.plus if isinstance( test, PlusInterfaceOperator) else target.minus + trial_target = target.plus if isinstance(trial, PlusInterfaceOperator) else target.minus + test_ext = test_target.ext + trial_ext = trial_target.ext + ncells = tuple(max(i, j) for i, j in zip(test_space.ncells, trial_space.ncells)) + if isinstance(trial_space, VectorFemSpace): + spaces = [] + for sp in trial_space.spaces: + if (trial_target.axis, trial_target.ext) in sp.interfaces: + spaces.append(sp.get_refined_space(ncells).interfaces[trial_target.axis, trial_target.ext]) + + if len(spaces) == len(trial_space.spaces): + sym_space = trial_space.symbolic_space + trial_space = VectorFemSpace(*spaces) + trial_space.symbolic_space = sym_space + + elif (trial_target.axis, trial_target.ext) in trial_space.interfaces: + sym_space = trial_space.symbolic_space + trial_space = trial_space.get_refined_space(ncells).interfaces[trial_target.axis, trial_target.ext] + trial_space.symbolic_space = sym_space + + test_space = test_space.get_refined_space(ncells) + self._test_ext = test_target.ext + self._trial_ext = trial_target.ext + + #... + + # Assuming that all vector spaces (and their Cartesian decomposition, + # if any) are compatible with each other, extract the first available + # vector space from which (starts, ends, npts) will be read: + starts = coeff_space.starts + ends = coeff_space.ends + npts = coeff_space.npts + + # MPI communicator + comm = coeff_space.cart.comm if coeff_space.parallel else None + + # Store the MPI communicator (or None) + self._comm = comm + + #... + # Get default backend from environment, or use 'python'. + default_backend = PSYDAC_BACKENDS.get(os.environ.get('PSYDAC_BACKEND'))\ + or PSYDAC_BACKENDS['python'] + + # Backends for code generation + assembly_backend = backend or assembly_backend + linalg_backend = backend or linalg_backend + + # Store backend dictionary + self._backend = assembly_backend or default_backend + #... + + # TODO: remove + # BasicDiscrete generates the assembly code and sets the following attributes that are used afterwards: + # self._func, self._free_args, self._max_nderiv and self._backend +# BasicDiscrete.__init__(self, expr, kernel_expr, comm=comm, root=0, discrete_space=discrete_space, +# nquads=nquads, is_rational_mapping=is_rational_mapping, mapping=symbolic_mapping, +# mapping_space=mapping_space, num_threads=self._num_threads, backend=assembly_backend) + + + #... Compute the string with all the imports + texpr = kernel_expr + sym_expr = SymbolicExpr(texpr.expr) + imports = compute_imports(sym_expr, spaces=(trial_space, test_space), openmp=False) + indent = 4 + glue = '\n' + ' '* indent + imports_str = glue.join([f"from {m} import {', '.join(vars)}" + for m, vars in imports.items()]) + + # Broadcast the import information (sqrt, sin, pi, ...) to all processes + if (comm is not None) and (comm.size > 1): + imports_str = comm.bcast(imports_str, root=0) + + # Store the imports string as it will be used by make_file() + self._imports_string = imports_str + #... + + # Compute the highest order of derivation in the kernel expression + self._max_nderiv = compute_max_nderiv(kernel_expr) + + # TODO [YG 31.07.2025]: Implement this + self._free_args = compute_free_arguments(expr, kernel_expr) + + #... Handle the special case where the current MPI process does not need to do anything + if isinstance(target, (Boundary, Interface)): + + # If process does not own the boundary or interface, do not assemble anything + if test_ext == -1: + if starts[axis] != 0: + self._func = do_nothing + + elif test_ext == 1: + if ends[axis] != npts[axis]-1: + self._func = do_nothing + + # In case of target==Interface, we only use the MPI ranks that are on the interface to assemble the BilinearForm + if self._func == do_nothing and isinstance(target, Interface): + self._free_args = () + self._args = () + self._global_matrices = () + self._threads_args = () + return + #... + + #... Build the quadrature grids + if isinstance(target, Boundary): + test_grid = QuadratureGrid( test_space, axis=axis, ext= test_ext, nquads=nquads) + trial_grid = QuadratureGrid(trial_space, axis=axis, ext=trial_ext, nquads=nquads) + self._grid = (test_grid,) + elif isinstance(target, Interface): + # this part treats the cases of: + # integral(v_minus * u_plus) + # integral(v_plus * u_minus) + # the other cases, integral(v_minus * u_minus) and integral(v_plus * u_plus) + # are converted to boundary integrals by Sympde + test_grid = QuadratureGrid( test_space, axis=axis, ext= test_ext, nquads=nquads) + trial_grid = QuadratureGrid(trial_space, axis=axis, ext=trial_ext, nquads=nquads) + self._grid = (test_grid, trial_grid) if test_target == target.minus else (trial_grid, test_grid) + self._test_ext = test_target.ext + self._trial_ext = trial_target.ext + else: + test_grid = QuadratureGrid( test_space, nquads=nquads) + trial_grid = QuadratureGrid(trial_space, nquads=nquads) + self._grid = (test_grid,) + #... + + # Extract the basis function values on the quadrature grids + self._test_basis = BasisValues( + test_space, + nderiv = self.max_nderiv, + nquads = nquads, + trial = False, + grid = test_grid + ) + self._trial_basis = BasisValues( + trial_space, + nderiv = self.max_nderiv, + nquads = nquads, + trial = True , + grid = trial_grid + ) + + # Allocate the output matrix, if needed + self.allocate_matrices(linalg_backend) + + # Determine whether OpenMP instructions were generated + self._with_openmp = (assembly_backend['name'] == 'pyccel' and assembly_backend['openmp']) if assembly_backend else False + + # Construct the arguments to be passed to the assemble() function, which is stored in self._func + # First we generate the assembly file + + # pyccelize process of computing the test_trial arrays + # currently set to False, as a Python 3.9 test fails, and due to the "speed up" not being significant + self._pyccelize_test_trial_computation = False + + # no openmp support yet: with_openmp is not passed + self._args, self._threads_args = self.construct_arguments_generate_assembly_file() + + #-------------------------------------------------------------------------- + @property + def comm(self): + return self._comm + + @property + def expr(self): + return self._expr + + @property + def kernel_expr(self): + return self._kernel_expr + + @property + def domain(self): + return self._domain + + @property + def mapping(self): + return self._mapping + + @property + def is_rational_mapping(self): + return self._is_rational_mapping + + @property + def target(self): + return self._target + + @property + def spaces(self): + return self._spaces + + @property + def test_basis(self): + return self._test_basis + + @property + def trial_basis(self): + return self._trial_basis + + @property + def grid(self): + return self._grid + + @property + def nquads(self): + return self._grid[0].nquads + + @property + def free_args(self): + return self._free_args + + @property + def max_nderiv(self): + # TODO: compute with read_BilinearForm and store + return self._max_nderiv + + @property + def backend(self): + return self._backend + + @property + def args(self): + return self._args + + @property + def global_matrices(self): + return self._global_matrices + + #-------------------------------------------------------------------------- + def allocate_matrices(self, backend=None): + """ + Allocate the global matrices used in the assembly method. + In this method we allocate only the matrices that are computed in the self._target domain, + we also avoid double allocation if we have many DiscreteLinearForm that are defined on the same self._target domain. + + Parameters + ---------- + backend : dict + The backend used to accelerate the computing kernels. + + """ + global_mats = {} + + expr = self.kernel_expr.expr + target = self.kernel_expr.target + test_degree = np.array(self.test_basis.space.degree) + trial_degree = np.array(self.trial_basis.space.degree) + test_space = self.spaces[1].coeff_space + trial_space = self.spaces[0].coeff_space + test_fem_space = self.spaces[1] + trial_fem_space = self.spaces[0] + domain = self.domain + is_broken = len(domain) > 1 + is_conformal = True + + if isinstance(expr, (ImmutableDenseMatrix, Matrix)): + if not isinstance(test_degree[0],(list, tuple, np.ndarray)): + test_degree = [test_degree] + + if not isinstance(trial_degree[0],(list, tuple, np.ndarray)): + trial_degree = [trial_degree] + + pads = np.empty((len(test_degree),len(trial_degree),len(test_degree[0])), dtype=int) + for i in range(len(test_degree)): + for j in range(len(trial_degree)): + td = test_degree[i] + trd = trial_degree[j] + pads[i,j][:] = np.array([td, trd]).max(axis=0) + else: + pads = np.maximum(test_degree, trial_degree) + + if self._matrix is None and (is_broken or isinstance(expr, (ImmutableDenseMatrix, Matrix))): + self._matrix = BlockLinearOperator(trial_space, test_space) + + if is_broken: + i, j = self.get_space_indices_from_target(domain, target) + test_fem_space = self.spaces[1].spaces[i] + trial_fem_space = self.spaces[0].spaces[j] + test_space = test_space.spaces[i] + trial_space = trial_space.spaces[j] + ncells = tuple(max(i,j) for i,j in zip(test_fem_space.ncells, trial_fem_space.ncells)) + is_conformal = tuple(test_fem_space.ncells) == ncells and tuple(trial_fem_space.ncells) == ncells + if is_broken and not is_conformal and not i==j: + use_restriction = all(trn>=tn for trn,tn in zip(trial_fem_space.ncells, test_fem_space.ncells)) + use_prolongation = not use_restriction + + else: + ncells = tuple(max(i,j) for i,j in zip(test_fem_space.ncells, trial_fem_space.ncells)) + i=0 + j=0 + #else so initialisation causing bug on line 682 + + if isinstance(expr, (ImmutableDenseMatrix, Matrix)): # case of system of equations + + if is_broken: #multi patch + if not self._matrix[i,j]: + mat = BlockLinearOperator(trial_fem_space.get_refined_space(ncells).coeff_space, test_fem_space.get_refined_space(ncells).coeff_space) + if not is_conformal and not i==j: + if use_restriction: + Ps = [knot_insertion_projection_operator(ts.get_refined_space(ncells), ts) for ts in test_fem_space.spaces] + P = BlockLinearOperator(test_fem_space.get_refined_space(ncells).coeff_space, test_fem_space.coeff_space) + for ni,Pi in enumerate(Ps): + P[ni,ni] = Pi + + mat = ComposedLinearOperator(trial_space, test_space, P, mat) + + elif use_prolongation: + Ps = [knot_insertion_projection_operator(trs, trs.get_refined_space(ncells)) for trs in trial_fem_space.spaces] + P = BlockLinearOperator(trial_fem_space.coeff_space, trial_fem_space.get_refined_space(ncells).coeff_space) + for ni,Pi in enumerate(Ps): + P[ni,ni] = Pi + + mat = ComposedLinearOperator(trial_space, test_space, mat, P) + + self._matrix[i,j] = mat + + matrix = self._matrix[i,j] + else: # single patch + matrix = self._matrix + + shape = expr.shape + for k1 in range(shape[0]): + for k2 in range(shape[1]): + if expr[k1,k2].is_zero: + continue + + if isinstance(test_fem_space, VectorFemSpace): + ts_space = test_fem_space.get_refined_space(ncells).coeff_space.spaces[k1] + else: + ts_space = test_fem_space.get_refined_space(ncells).coeff_space + + if isinstance(trial_fem_space, VectorFemSpace): + tr_space = trial_fem_space.get_refined_space(ncells).coeff_space.spaces[k2] + else: + tr_space = trial_fem_space.get_refined_space(ncells).coeff_space + + if is_conformal and matrix[k1, k2]: + global_mats[k1, k2] = matrix[k1, k2] + elif not i == j: # assembling in an interface (type(target) == Interface) + axis = target.axis + ext_d = self._trial_ext + ext_c = self._test_ext + test_n = self. test_basis.space.spaces[k1].spaces[axis].nbasis + test_s = self. test_basis.space.spaces[k1].coeff_space.starts[axis] + trial_n = self.trial_basis.space.spaces[k2].spaces[axis].nbasis + cart = self.trial_basis.space.spaces[k2].coeff_space.cart + trial_s = cart.global_starts[axis][cart._coords[axis]] + + s_d = trial_n - trial_s - trial_degree[k2][axis] - 1 if ext_d == 1 else 0 + s_c = test_n - trial_s - test_degree[k1][axis] - 1 if ext_c == 1 else 0 + + # We only handle the case where direction = 1 + direction = target.ornt + if domain.dim == 2: + assert direction == 1 + elif domain.dim == 3: + assert all(d==1 for d in direction) + + direction = 1 + flip = [direction]*domain.dim + flip[axis] = 1 + if self._func != do_nothing: + global_mats[k1, k2] = StencilInterfaceMatrix(tr_space, ts_space, + s_d, s_c, + axis, axis, + ext_d, ext_c, + pads=tuple(pads[k1, k2]), + flip=flip) + else: + global_mats[k1, k2] = StencilMatrix(tr_space, ts_space, pads = tuple(pads[k1, k2])) + + if is_conformal: + matrix[k1, k2] = global_mats[k1, k2] + elif use_restriction: + matrix.multiplicants[-1][k1, k2] = global_mats[k1, k2] + elif use_prolongation: + matrix.multiplicants[0][k1, k2] = global_mats[k1, k2] + + else: # case of scalar equation + if is_broken: # multi-patch + if self._matrix[i, j]: + global_mats[i, j] = self._matrix[i, j] + + elif not i == j: # assembling in an interface (type(target) == Interface) + axis = target.axis + ext_d = self._trial_ext + ext_c = self._test_ext + test_n = self.test_basis.space.spaces[axis].nbasis + test_s = self.test_basis.space.coeff_space.starts[axis] + trial_n = self.trial_basis.space.spaces[axis].nbasis + cart = self.trial_basis.space.coeff_space.cart + trial_s = cart.global_starts[axis][cart._coords[axis]] + + s_d = trial_n - trial_s - trial_degree[axis] - 1 if ext_d == 1 else 0 + s_c = test_n - trial_s - test_degree[axis] - 1 if ext_c == 1 else 0 + + # We only handle the case where direction = 1 + direction = target.ornt + if domain.dim == 2: + assert direction == 1 + elif domain.dim == 3: + assert all(d==1 for d in direction) + + direction = 1 + flip = [direction]*domain.dim + flip[axis] = 1 + + if self._func != do_nothing: + mat = StencilInterfaceMatrix(trial_fem_space.get_refined_space(ncells).coeff_space, + test_fem_space.get_refined_space(ncells).coeff_space, + s_d, s_c, + axis, axis, + ext_d, ext_c, + flip=flip) + if not is_conformal: + if use_restriction: + P = knot_insertion_projection_operator(test_fem_space.get_refined_space(ncells), test_fem_space) + mat = ComposedLinearOperator(trial_space, test_space, P, mat) + elif use_prolongation: + P = knot_insertion_projection_operator(trial_fem_space, trial_fem_space.get_refined_space(ncells)) + mat = ComposedLinearOperator(trial_space, test_space, mat, P) + + global_mats[i, j] = mat + + # define part of the global matrix as a StencilMatrix + else: + global_mats[i, j] = StencilMatrix(trial_space, test_space, pads=tuple(pads)) + + if (i, j) in global_mats: + self._matrix[i, j] = global_mats[i, j] + + + # in single patch case, we define the matrices needed for the patch + else: + if self._matrix: + global_mats[0, 0] = self._matrix + else: + global_mats[0, 0] = StencilMatrix(trial_space, test_space, pads=tuple(pads)) + + self._matrix = global_mats[0, 0] + + # Set the backend of our matrices if given + if backend is not None and is_broken: + for mat in global_mats.values(): + mat.set_backend(backend) + elif backend is not None: + self._matrix.set_backend(backend) + + self._global_matrices = [M._data for M in extract_stencil_mats(global_mats.values())] + + #-------------------------------------------------------------------------- + def assemble(self, *, reset=True, **kwargs): + """ + This method assembles the left hand side Matrix by calling the private method `self._func` with proper arguments. + + In the complex case, this function returns the matrix conjugate. This comes from the fact that the + problem `a(u,v)=b(v)` is discretized as `A @ conj(U) = B` due to the antilinearity of `a` in the first variable. + Thus, to obtain `U`, the assemble function returns `conj(A)`. + + TODO: remove these lines when the dot product is changed for complex. + For now, since the dot product does not compute the conjugate in the complex case. We do not use the conjugate in the assemble function. + It should work if the complex only comes from the `rhs` in the linear form. + """ + + if self._free_args: + basis = [] + spans = [] + degrees = [] + pads = [] + coeffs = [] + consts = [] + + for key in self._free_args: + v = kwargs[key] + + if len(self.domain) > 1 and isinstance(v, FemField) and (v.space.is_multipatch or v.space.is_vector_valued): + assert v.space.is_multipatch ## [MCP 27.03.2025] should hold since len(domain) > 1. If Ok we can simplify above if + i, j = self.get_space_indices_from_target(self.domain, self.target) + assert i == j + v = v[i] + if isinstance(v, FemField): + assert len(self.grid) == 1 + if not v.coeffs.ghost_regions_in_sync: + v.coeffs.update_ghost_regions() + basis_v = BasisValues( + v.space, + nderiv = self.max_nderiv, + nquads = self.nquads, + trial = True, + grid = self.grid[0] + ) + bs, d, s, p, mult = construct_test_space_arguments(basis_v) + basis += bs + spans += s + degrees += [np.int64(a) for a in d] + pads += [np.int64(a) for a in p] + if v.space.is_multipatch or v.space.is_vector_valued: + coeffs += (e._data for e in v.coeffs) + else: + coeffs += (v.coeffs._data, ) + else: + consts += (v, ) + + args = (*self.args, *basis, *spans, *degrees, *pads, *coeffs, *consts) + + else: + args = self._args + + if reset: + reset_arrays(*self.global_matrices) + + self._func(*args, *self._threads_args) + if self._matrix and self._update_ghost_regions: + self._matrix.exchange_assembly_data() + + # TODO : uncomment this line when the conjugate is applied on the dot product in the complex case + #self._matrix.conjugate(out=self._matrix) + + if self._matrix: + self._matrix.ghost_regions_in_sync = False + + return self._matrix + + #-------------------------------------------------------------------------- + @property + def _assembly_template_head(self): + """A template for the 'head' of the assembly function. Only used with the sum factorization algorithm.""" + code = '''def assemble_matrix_{FILE_ID}({MAPPING_PART_1} +{SPAN} {MAPPING_PART_2} + global_x1 : "float64[:,:]", global_x2 : "float64[:,:]", global_x3 : "float64[:,:]", + {MAPPING_PART_3} + n_element_1 : "int64", n_element_2 : "int64", n_element_3 : "int64", + nq1 : "int64", nq2 : "int64", nq3 : "int64", + pad1 : "int64", pad2 : "int64", pad3 : "int64", + {MAPPING_PART_4} +{G_MAT}{NEW_ARGS}{FIELD_ARGS}): + + from numpy import abs as Abs + {imports} +''' + return code + + #-------------------------------------------------------------------------- + @property + def _assembly_template_body_bspline(self): + """A template for the 'body' of the assembly function (when using a spline mapping). Only used with the sum factorization algorithm.""" + code = ''' + arr_coeffs_x = zeros((1 + test_mapping_p1, 1 + test_mapping_p2, 1 + test_mapping_p3), dtype='float64') + arr_coeffs_y = zeros((1 + test_mapping_p1, 1 + test_mapping_p2, 1 + test_mapping_p3), dtype='float64') + arr_coeffs_z = zeros((1 + test_mapping_p1, 1 + test_mapping_p2, 1 + test_mapping_p3), dtype='float64') + +{F_COEFFS_ZEROS} + +{KEYS} + for k_1 in range(n_element_1): + span_mapping_1 = global_span_mapping_1[k_1] +{LOCAL_SPAN}{F_SPAN_1}{A1} + for q_1 in range(nq1): + for k_2 in range(n_element_2): + span_mapping_2 = global_span_mapping_2[k_2] +{F_SPAN_2} + for q_2 in range(nq2): + for k_3 in range(n_element_3): + span_mapping_3 = global_span_mapping_3[k_3] +{F_SPAN_3}{F_COEFFS} + arr_coeffs_x[:,:,:] = global_arr_coeffs_x[test_mapping_p1 + span_mapping_1 - test_mapping_p1:test_mapping_p1 + 1 + span_mapping_1,test_mapping_p2 + span_mapping_2 - test_mapping_p2:test_mapping_p2 + 1 + span_mapping_2,test_mapping_p3 + span_mapping_3 - test_mapping_p3:test_mapping_p3 + 1 + span_mapping_3] + arr_coeffs_y[:,:,:] = global_arr_coeffs_y[test_mapping_p1 + span_mapping_1 - test_mapping_p1:test_mapping_p1 + 1 + span_mapping_1,test_mapping_p2 + span_mapping_2 - test_mapping_p2:test_mapping_p2 + 1 + span_mapping_2,test_mapping_p3 + span_mapping_3 - test_mapping_p3:test_mapping_p3 + 1 + span_mapping_3] + arr_coeffs_z[:,:,:] = global_arr_coeffs_z[test_mapping_p1 + span_mapping_1 - test_mapping_p1:test_mapping_p1 + 1 + span_mapping_1,test_mapping_p2 + span_mapping_2 - test_mapping_p2:test_mapping_p2 + 1 + span_mapping_2,test_mapping_p3 + span_mapping_3 - test_mapping_p3:test_mapping_p3 + 1 + span_mapping_3] + for q_3 in range(nq3): + x = 0.0 + y = 0.0 + z = 0.0 + + x_x1 = 0.0 + x_x2 = 0.0 + x_x3 = 0.0 + y_x1 = 0.0 + y_x2 = 0.0 + y_x3 = 0.0 + z_x1 = 0.0 + z_x2 = 0.0 + z_x3 = 0.0 +{D2_1} + +{F_INIT} + +{F_ASSIGN_LOOP} + + for i_1 in range(test_mapping_p1+1): + mapping_1 = global_basis_mapping_1[k_1, i_1, 0, q_1] + mapping_1_x1 = global_basis_mapping_1[k_1, i_1, 1, q_1] + {D2_2} + for i_2 in range(test_mapping_p2+1): + mapping_2 = global_basis_mapping_2[k_2, i_2, 0, q_2] + mapping_2_x2 = global_basis_mapping_2[k_2, i_2, 1, q_2] + {D2_3} + for i_3 in range(test_mapping_p3+1): + mapping_3 = global_basis_mapping_3[k_3, i_3, 0, q_3] + mapping_3_x3 = global_basis_mapping_3[k_3, i_3, 1, q_3] + {D2_4} + + coeff_x = arr_coeffs_x[i_1,i_2,i_3] + coeff_y = arr_coeffs_y[i_1,i_2,i_3] + coeff_z = arr_coeffs_z[i_1,i_2,i_3] + + mapping = mapping_1*mapping_2*mapping_3 + mapping_x1 = mapping_1_x1*mapping_2*mapping_3 + mapping_x2 = mapping_1*mapping_2_x2*mapping_3 + mapping_x3 = mapping_1*mapping_2*mapping_3_x3 + +{D2_5} + + x += mapping*coeff_x + y += mapping*coeff_y + z += mapping*coeff_z + + x_x1 += mapping_x1*coeff_x + x_x2 += mapping_x2*coeff_x + x_x3 += mapping_x3*coeff_x + y_x1 += mapping_x1*coeff_y + y_x2 += mapping_x2*coeff_y + y_x3 += mapping_x3*coeff_y + z_x1 += mapping_x1*coeff_z + z_x2 += mapping_x2*coeff_z + z_x3 += mapping_x3*coeff_z + +{D2_6} + +{TEMPS} +{COUPLING_TERMS} +''' + return code + + #-------------------------------------------------------------------------- + @property + def _assembly_template_body_analytic(self): + """A template for the 'body' of the assembly function (when using an analytic or no mapping). Only used with the sum factorization algorithm.""" + code = ''' + local_x1 = zeros_like(global_x1[0,:]) + local_x2 = zeros_like(global_x2[0,:]) + local_x3 = zeros_like(global_x3[0,:]) + +{F_COEFFS_ZEROS} + +{KEYS} + for k_1 in range(n_element_1): + local_x1[:] = global_x1[k_1,:] +{LOCAL_SPAN}{F_SPAN_1}{A1} + for q_1 in range(nq1): + x1 = local_x1[q_1] + for k_2 in range(n_element_2): + local_x2[:] = global_x2[k_2,:] +{F_SPAN_2} + for q_2 in range(nq2): + x2 = local_x2[q_2] + for k_3 in range(n_element_3): + local_x3[:] = global_x3[k_3,:] +{F_SPAN_3}{F_COEFFS} + for q_3 in range(nq3): + x3 = local_x3[q_3] + +{F_INIT} + +{F_ASSIGN_LOOP} + +{TEMPS} +{COUPLING_TERMS} +''' + return code + + #-------------------------------------------------------------------------- + @property + def _assembly_template_loop(self): + """A template for the 'loop' of the assembly function. Only used with the sum factorization algorithm.""" + code = ''' + {A2}[:] = 0.0 + for k_2 in range(n_element_2): + {SPAN_2} = {GLOBAL_SPAN_2}[k_2] + for q_2 in range(nq2): + {A3}[:] = 0.0 + for k_3 in range(n_element_3): + {SPAN_3} = {GLOBAL_SPAN_3}[k_3] + for q_3 in range(nq3): + a4 = {COUPLING_TERMS}[k_2, q_2, k_3, q_3, :] + for i_3 in range({TEST_V_P3} + 1): + for j_3 in range({TRIAL_U_P3} + 1): + for e in range({NEXPR}): + {A3}[e, {SPAN_3} - {TEST_V_P3} + i_3, {MAX_P3} - {I_3} + j_3] += {TEST_TRIAL_3}[k_3, q_3, i_3, j_3, {KEYS_3}[2*e], {KEYS_3}[2*e+1]] * a4[e] + for i_2 in range({TEST_V_P2} + 1): + for j_2 in range({TRIAL_U_P2} + 1): + for e in range({NEXPR}): + {A2}[e, {SPAN_2} - {TEST_V_P2} + i_2, :, {MAX_P2} - {I_2} + j_2, :] += {TEST_TRIAL_2}[k_2, q_2, i_2, j_2, {KEYS_2}[2*e], {KEYS_2}[2*e+1]] * {A3}[e,:,:] + for i_1 in range({TEST_V_P1} + 1): + for j_1 in range({TRIAL_U_P1} + 1): + {A1}[i_1, :, :, {MAX_P1} - {I_1} + j_1, :, :] += {A2_TEMP} +''' + return code + + #-------------------------------------------------------------------------- + def make_file(self, temps, ordered_stmts, field_derivatives, max_logical_derivative, test_mult, trial_mult, test_v_p, trial_u_p, keys_1, keys_2, keys_3, mapping_option): + """ + Part of the sum factorization algorithm implementation. + Generates the correct assembly file. + Used at the end of construct_arguments_generate_assembly_file, before eventually pyccelizing that file. + + Parameters + ---------- + temps : tuple + Tuple of Assign statements defining temporary values. + Arithmetic combinations of these make up the coupling terms. + + ordered_stmts : dict + Dictionary defining the coupling terms. Keys are combinations of + test and trial function components, values are Assign statements + in terms of temporaries appearing in temps. + + field_derivatives : dict + Dictionary containing information on the derivatives of free FemFields. + Keys are components of free FemFields. Values are dictionaries again. + Their keys are names, as appearing in the assembly file, of partial derivatives of the + corresponding FemField component, and their values are dictionaries again. + Example: {F1_0_x3 : {'x1': 0, 'x2': 0, 'x3': 1}, F1_0_x2 : ...} + Meaning: There exists a free FemField named F1. Among other, the partial derivative w.r.t. x3 + of its first component F1_0 appears. + + max_logical_derivative : int + The largest appearing derivative order. + + test_mult : list + List of length 3(scalar test function) or 9(vector test function) including multiplicity information. + + trial_mult : list + List of length 3(scalar trial function) or 9(vector trial function) including multiplicity information. + + test_v_p : dict + Dictionary of length 1(scalar test function) or length 3(vector test function). + Each key corresponds to a component of the funciton (space), and each corresponding value + is a list of Bspline degrees of length 3. Example: Discretizing a de de Rham sequence using + a degree vector [2, 3, 4] means that test_v_p for a test function belonging to H(curl) will be + {0: [1, 3, 4], 1: [2, 2, 4], 2: [2, 3, 3]} + + trial_u_p : dict + Dictionary of length 1(scalar trial function) or length 3(vector trial function). + Each key corresponds to a component of the funciton (space), and each corresponding value + is a list of Bspline degrees of length 3. Example: Discretizing a de de Rham sequence using + a degree vector [2, 3, 4] means that trial_u_p for a trial function belonging to H^1 will be + {0: [2, 3, 4]} + + keys_1 : dict + Dictionary relating subexpressions to x1-derivative combinations. + Keys are combinations of test and trial function components. + Values are lists, each entry corresponding to one appearing partial derivative + combination of these components. + Example: keys_1[(u[0], v[1])][3] = [1,0] means that the fourth ([3]) + sub-expression (partial derivative combination) corresponding to the trial-test-function-component-product + u[0] * v[1] involves a first derivative in x1 direction of the trial function + and no derivative in x1 direction of the test function. + Information on appearing partial derivatives in x2 and x3 direction is stored in keys_2 and keys_3. + + keys_2 : dict + See keys_1. + + keys_3 : dict + See keys_1. + + mapping_option : None | 'Bspline' + None in case of no mapping or an analytical mapping, 'Bspline' in case of a Bspline mapping. + + Returns + ------- + + file_id : str + random string of length 8, corresponding to the assembly file name located in __psydac__/ + + """ + + #------------------------- FILE_ID ------------------------- + comm = self.comm + + # Root process generates a random string to be used as file_id + if comm is None or comm.rank == 0: + file_id = random_string(size=8) + else: + file_id = None + + # Parallel case: root process broadcasts file_id to all processes + if comm is not None and comm.size > 1: + file_id = comm.bcast(file_id, root=0) + + # ----- free FemField related strings ----- + + # used as {FIELD_ARGS} in _assembly_template_head + # adding the right arguments for free FemFields to the assembly function header + basis_args_block = [f'global_test_basis_'+'{field}'+f'_{i+1} : "float64[:,:,:,:]"' for i in range(3)] + basis_args_block = ", ".join(basis_args_block) + "," + basis_args_block = [basis_args_block.format(field=field) for field in field_derivatives] + basis_args = " " + "\n ".join(basis_args_block) + "\n" + span_args_block = [f'global_span_'+'{field}'+f'_{i+1} : "int64[:]"' for i in range(3)] + span_args_block = ", ".join(span_args_block) + "," + span_args_block = [span_args_block.format(field=field) for field in field_derivatives] + span_args = " " + "\n ".join(span_args_block) + "\n" + degree_args_block = [f'test_'+'{field}'+f'_p{i+1} : "int64"' for i in range(3)] + degree_args_block = ", ".join(degree_args_block) + "," + degree_args_block = [degree_args_block.format(field=field) for field in field_derivatives] + degree_args = " " + "\n ".join(degree_args_block) + "\n" + pad_args_block = [f'pad_'+'{field}'+f'_{i+1} : "int64"' for i in range(3)] + pad_args_block = ", ".join(pad_args_block) + "," + pad_args_block = [pad_args_block.format(field=field) for field in field_derivatives] + pad_args = " " + "\n ".join(pad_args_block) + "\n" + coeff_args_block = [f'global_arr_coeffs_{field} : "float64[:,:,:]"' for field in field_derivatives] + coeff_args = " " + ", ".join(coeff_args_block) + FIELD_ARGS = basis_args+span_args+degree_args+pad_args+coeff_args + + # {F_COEFFS_ZEROS} in both _assembly_template_body_bspline & _analytic + F_COEFFS_ZEROS = "\n".join([f" arr_coeffs_{field} = zeros((1 + test_{field}_p1, 1 + test_{field}_p2, 1 + test_{field}_p3), dtype='float64')" for field in field_derivatives]) + + # {F_SPAN_1}, {F_SPAN_2}, {F_SPAN_3} in both _assembly_template_body_bspline & _analytic + F_SPAN_1 = "\n".join([f" span_{field}_1 = global_span_{field}_1[k_1]" for field in field_derivatives]) + "\n" + F_SPAN_2 = "\n".join([f" span_{field}_2 = global_span_{field}_2[k_2]" for field in field_derivatives]) + "\n" + F_SPAN_3 = "\n".join([f" span_{field}_3 = global_span_{field}_3[k_3]" for field in field_derivatives]) + "\n" + + # {F_COEFFS} in both _assembly_template_body_bspline & _analytic + coeff_ranges = ", ".join([f"pad_"+"{field}"+f"_{i+1} + span_"+"{field}"+f"_{i+1} - test_"+"{field}"+f"_p{i+1}:1 + pad_"+"{field}"+f"_{i+1} + span_"+"{field}"+f"_{i+1}" for i in range(3)]) + F_COEFFS = "\n".join([f" arr_coeffs_{field}[:,:,:] = global_arr_coeffs_{field}[{coeff_ranges.format(field=field)}]" for i, field in enumerate(field_derivatives)]) + + # {F_INIT} + F_INIT = "\n".join([f" {derivative} = 0.0" for field in field_derivatives for derivative in field_derivatives[field]]) + + # + # field_init assigns 0 to appearing free FemField derivatives (F_x1 = 0.0 \n F_x2 = 0.0 \n ...) + # In the following, we assemble loops that correctly compute those free FemField derivatives at + # a specific quadrature point (q_1, q_2, q_3). Those values will then be used in the computation + # of the temps or directly in the computation of the coupling terms + # + assign_loop_contents = {'1':{}, '2':{}, '3':{}} + multiplication_info = {} + + for field, derivatives in field_derivatives.items(): + multiplication_info[field] = {} + assign_statements = {'1':[], '2':[], '3':[]} + for derivative, dxs in derivatives.items(): + multiplication_info[field][derivative] = [] + dx1 = dxs['x1'] + dx2 = dxs['x2'] + dx3 = dxs['x3'] + for i, dx in enumerate([dx1, dx2, dx3]): + name = f"{field}_{i+1}" if dx == 0 else f"{field}_{i+1}_{dx*f'x{i+1}'}" + multiplication_info[field][derivative].append(name) + if dx == 0: + assign_statement = f"{name} = global_test_basis_{field}_{i+1}[k_{i+1}, i_{i+1}, 0, q_{i+1}]" + else: + assign_statement = f"{name} = global_test_basis_{field}_{i+1}[k_{i+1}, i_{i+1}, {dx}, q_{i+1}]" + if assign_statement not in assign_statements[f"{i+1}"]: + assign_statements[f"{i+1}"].append(assign_statement) + for i in range(3): + content = ("\n"+(8+i)*" ").join(assign_statements[f"{i+1}"]) + assign_loop_contents[f"{i+1}"][field] = content + tab = 7*" " + assign = [] + for field in field_derivatives: + txt = f"{tab}for i_1 in range(1 + test_{field}_p1):\n" + \ + f"{tab} {assign_loop_contents['1'][field]}\n" + \ + f"{tab} for i_2 in range(1 + test_{field}_p2):\n" + \ + f"{tab} {assign_loop_contents['2'][field]}\n" + \ + f"{tab} for i_3 in range(1 + test_{field}_p3):\n" + \ + f"{tab} {assign_loop_contents['3'][field]}\n" + \ + f"{tab} coeff_{field} = arr_coeffs_{field}[i_1, i_2, i_3]\n" + for derivative in multiplication_info[field]: + factors = " * ".join(multiplication_info[field][derivative]) + txt += f"{tab} {derivative} += {factors} * coeff_{field}\n" + txt += "\n" + assign.append(txt) + + # {F_ASSIGN_LOOP} in both _assembly_template_body_bspline & _analytic + F_ASSIGN_LOOP = "\n".join(assign) + + # ----------------------------------------- + + # ----- load the templates ----- + # + # head for the function header and imports + # body for the computation of coupling terms + # loop (part of function body): one loop per block ( e.g. (u[0], v[1]) ), each loop effectively + # assembles one StencilMatrix per sub expression ( e.g. (dx1(u[0]), dx3(v[1])) ) + code_head = self._assembly_template_head + code_loop = self._assembly_template_loop + if mapping_option == 'Bspline': + code_body = self._assembly_template_body_bspline + else: + code_body = self._assembly_template_body_analytic + # ------------------------------ + + # ---- obtain basic information not explicitely passed in the args ----- + blocks = ordered_stmts.keys() + block_list = list(blocks) + trial_components = [block[0] for block in block_list] + test_components = [block[1] for block in block_list] + nu = len(set(trial_components)) + nv = len(set(test_components)) + d = 3 + assert d == 3 + # ---------------------------------------------------------------------- + + # Prepare strings and string templates depending on whether the trial and test function are vector-valued or not (nu, nv > 1 or == 1) + + # ------------------------- STRINGS HEAD ------------------------- + + global_span_v_str = 'global_span_v_{v_j}_' if nv > 1 else 'global_span_v_' + + if mapping_option == 'Bspline': + MAPPING_PART_1 = 'global_basis_mapping_1 : "float64[:,:,:,:]", global_basis_mapping_2 : "float64[:,:,:,:]", global_basis_mapping_3 : "float64[:,:,:,:]", ' + MAPPING_PART_2 = 'global_span_mapping_1 : "int64[:]", global_span_mapping_2 : "int64[:]", global_span_mapping_3 : "int64[:]", ' + MAPPING_PART_3 = 'test_mapping_p1 : "int64", test_mapping_p2 : "int64", test_mapping_p3 : "int64", ' + MAPPING_PART_4 = 'global_arr_coeffs_x : "float64[:,:,:]", global_arr_coeffs_y : "float64[:,:,:]", global_arr_coeffs_z : "float64[:,:,:]", ' + else: + MAPPING_PART_1 = '' + MAPPING_PART_2 = '' + MAPPING_PART_3 = '' + MAPPING_PART_4 = '' + + if nv > 1: + tt1_str = 'test_trial_1_u_{u_i}_v_{v_j}' if nu > 1 else 'test_trial_1_u_v_{v_j}' + tt2_str = 'test_trial_2_u_{u_i}_v_{v_j}' if nu > 1 else 'test_trial_2_u_v_{v_j}' + tt3_str = 'test_trial_3_u_{u_i}_v_{v_j}' if nu > 1 else 'test_trial_3_u_v_{v_j}' + a3_str = 'a3_u_{u_i}_v_{v_j}' if nu > 1 else 'a3_u_v_{v_j}' + a2_str = 'a2_u_{u_i}_v_{v_j}' if nu > 1 else 'a2_u_v_{v_j}' + ct_str = 'coupling_terms_u_{u_i}_v_{v_j}' if nu > 1 else 'coupling_terms_u_v_{v_j}' + g_mat_str = 'g_mat_u_{u_i}_v_{v_j}' if nu > 1 else 'g_mat_u_v_{v_j}' + else: + tt1_str = 'test_trial_1_u_{u_i}_v' if nu > 1 else 'test_trial_1_u_v' + tt2_str = 'test_trial_2_u_{u_i}_v' if nu > 1 else 'test_trial_2_u_v' + tt3_str = 'test_trial_3_u_{u_i}_v' if nu > 1 else 'test_trial_3_u_v' + a3_str = 'a3_u_{u_i}_v' if nu > 1 else 'a3_u_v' + a2_str = 'a2_u_{u_i}_v' if nu > 1 else 'a2_u_v' + ct_str = 'coupling_terms_u_{u_i}_v' if nu > 1 else 'coupling_terms_u_v' + g_mat_str = 'g_mat_u_{u_i}_v' if nu > 1 else 'g_mat_u_v' + + # ---------------------------------------------------------------- + + # ------------------------- STRINGS BODY ------------------------- + + span_v_1_str = 'span_v_{v_j}_1' if nv > 1 else 'span_v_1' + test_v_p1_str = 'test_v_{v_j}_p1' if nv > 1 else 'test_v_p1' + + if nv > 1: + keys_2_str = 'keys_2_u_{u_i}_v_{v_j}' if nu > 1 else 'keys_2_u_v_{v_j}' + keys_3_str = 'keys_3_u_{u_i}_v_{v_j}' if nu > 1 else 'keys_3_u_v_{v_j}' + a1_str = 'a1_u_{u_i}_v_{v_j}' if nu > 1 else 'a1_u_v_{v_j}' + else: + keys_2_str = 'keys_2_u_{u_i}_v' if nu > 1 else 'keys_2_u_v' + keys_3_str = 'keys_3_u_{u_i}_v' if nu > 1 else 'keys_3_u_v' + a1_str = 'a1_u_{u_i}_v' if nu > 1 else 'a1_u_v' + + # ---------------------------------------------------------------- + + # ------------------------- STRINGS LOOP ------------------------- + + span_2_str = 'span_v_{v_j}_2' if nv > 1 else 'span_v_2' + span_3_str = 'span_v_{v_j}_3' if nv > 1 else 'span_v_3' + global_span_2_str = 'global_span_v_{v_j}_2' if nv > 1 else 'global_span_v_2' + global_span_3_str = 'global_span_v_{v_j}_3' if nv > 1 else 'global_span_v_3' + + # ---------------------------------------------------------------- + + #------------------------- MAKE HEAD ------------------------- + SPAN = '' + G_MAT = '' + + TT1 = ' ' + TT2 = ' ' + TT3 = ' ' + A3 = ' ' + A2 = ' ' + CT = ' ' + + for v_j in range(nv): + global_span_v = global_span_v_str.format(v_j=v_j) + SPAN += ' ' + for di in range(d): + SPAN += f'{global_span_v}{di+1} : "int64[:]", ' + SPAN = SPAN[:-1] + '\n' + + for block in blocks: + u_i = block[0].indices[0] if nu > 1 else 0 + v_j = block[1].indices[0] if nv > 1 else 0 + + # reverse order intended + if ((nu > 1) and (nv > 1)): + g_mat = g_mat_str.format(u_i=v_j, v_j=u_i) + else: + g_mat = g_mat_str.format(u_i=u_i, v_j=v_j) + G_MAT += f' {g_mat} : "float64[:,:,:,:,:,:]",\n' + + TT1 += tt1_str.format(u_i=u_i, v_j=v_j) + ' : "float64[:,:,:,:,:,:]", ' + TT2 += tt2_str.format(u_i=u_i, v_j=v_j) + ' : "float64[:,:,:,:,:,:]", ' + TT3 += tt3_str.format(u_i=u_i, v_j=v_j) + ' : "float64[:,:,:,:,:,:]", ' + A3 += a3_str.format(u_i=u_i, v_j=v_j) + ' : "float64[:,:,:]", ' + A2 += a2_str.format(u_i=u_i, v_j=v_j) + ' : "float64[:,:,:,:,:]", ' + CT += ct_str.format(u_i=u_i, v_j=v_j) + ' : "float64[:,:,:,:,:]", ' + + TT1 += '\n' + TT2 += '\n' + TT3 += '\n' + A3 += '\n' + A2 += '\n' + CT += '\n' + NEW_ARGS = TT1 + TT2 + TT3 + A3 + A2 + CT + IMPORTS = self._imports_string + + head = code_head.format(FILE_ID = file_id, + SPAN = SPAN, + G_MAT = G_MAT, + NEW_ARGS = NEW_ARGS, + MAPPING_PART_1 = MAPPING_PART_1, + MAPPING_PART_2 = MAPPING_PART_2, + MAPPING_PART_3 = MAPPING_PART_3, + MAPPING_PART_4 = MAPPING_PART_4, + FIELD_ARGS = FIELD_ARGS, + imports = IMPORTS) + + #------------------------- MAKE BODY ------------------------- + A1 = '' + KEYS_2 = '' + KEYS_3 = '' + LOCAL_SPAN = '' + TEMPS = '' + COUPLING_TERMS = '' + + for block in blocks: + u_i = block[0].indices[0] if nu > 1 else 0 + v_j = block[1].indices[0] if nv > 1 else 0 + + keys2 = keys_2[block].copy() + keys3 = keys_3[block].copy() + keys2 = ','.join(str(i) for i in keys2.flatten()) + keys3 = ','.join(str(i) for i in keys3.flatten()) + KEYS2 = keys_2_str.format(u_i=u_i, v_j=v_j) + KEYS3 = keys_3_str.format(u_i=u_i, v_j=v_j) + KEYS_2 += f' {KEYS2} = array([{keys2}])\n' + KEYS_3 += f' {KEYS3} = array([{keys3}])\n' + + test_v_p1, test_v_p2, test_v_p3 = test_v_p[v_j] + a1 = a1_str.format(u_i=u_i, v_j=v_j) + g_mat = g_mat_str.format(u_i=u_i, v_j=v_j) + TEST_V_P1 = test_v_p1_str.format(v_j=v_j) + SPAN_V_1 = span_v_1_str.format(v_j=v_j) + + A1_1 = f'{test_mult[0]}*pad1 + {SPAN_V_1} - {test_v_p1} : {test_mult[0]}*pad1 + {SPAN_V_1} + 1' if test_mult[0] > 1 else f'pad1 + {SPAN_V_1} - {test_v_p1} : pad1 + {SPAN_V_1} + 1' + A1_2 = f'{test_mult[1]}*pad2 : {test_mult[1]}*pad2 + n_element_2 + {test_v_p2} + ({test_mult[1]}-1)*(n_element_2-1)' if test_mult[1] > 1 else f'pad2 : pad2 + n_element_2 + {test_v_p2}' + A1_3 = f'{test_mult[2]}*pad3 : {test_mult[2]}*pad3 + n_element_3 + {test_v_p3} + ({test_mult[2]}-1)*(n_element_3-1)' if test_mult[2] > 1 else f'pad3 : pad3 + n_element_3 + {test_v_p3}' + A1 += f' {a1} = {g_mat}[{A1_1}, {A1_2}, {A1_3}, :, :, :]\n' + + for v_j in range(nv): + local_span_v_1 = span_v_1_str.format(v_j=v_j) + global_span_v = global_span_v_str.format(v_j=v_j) + LOCAL_SPAN += f' {local_span_v_1} = {global_span_v}1[k_1]\n' + + for temp in temps: + TEMPS += f' {temp.lhs} = {temp.rhs}\n' + for block in blocks: + for stmt in ordered_stmts[block]: + COUPLING_TERMS += f' {stmt.lhs} = {stmt.rhs}\n' + + KEYS = KEYS_2 + KEYS_3 + + # This part is interesting. Right now, below you find hardcoded rules regarding lines of code + # that need to be included when max_logical_derivative == 2 ( and mapping_option == 'Bspline'). + # E.g., the bilinear form corresponding to a bilaplacian problem ( laplace(laplace(u)) = f ) satisfies this assumption. + # This hardcoded set of rules could be generalized to n-th max derivatives - if needed! + # But for now, higher than second order derivatives on either trial or test function are not supported. + # + # Additional note: Given a Bspline mapping, the code computing the first order derivatives of mapping related terms is always required! + # Even in the case of a trivial bilinear form without derivatives. But only when there are second order partial derivatives involved + # do we need to compute second derivatives of the mapping (chain rule). + if (mapping_option == 'Bspline') and (max_logical_derivative == 2): + D2_1 = '\n' + spaces1 = ' ' + spaces2 = spaces1 + ' ' + for symbol in ('x', 'y', 'z'): + for d1 in range(1, 4): + for d2 in range(1, 4): + if d2 >= d1: + D2_1 += spaces1 + f'{symbol}_x{d1}x{d2} = 0.0\n' + D2_1 += '\n' + D2_2 = 'mapping_1_x1x1 = global_basis_mapping_1[k_1, i_1, 2, q_1]' + D2_3 = 'mapping_2_x2x2 = global_basis_mapping_2[k_2, i_2, 2, q_2]' + D2_4 = 'mapping_3_x3x3 = global_basis_mapping_3[k_3, i_3, 2, q_3]' + D2_5 = spaces2+'mapping_x1x1 = mapping_1_x1x1 * mapping_2 * mapping_3\n'+spaces2 + D2_5 += 'mapping_x1x2 = mapping_1_x1 * mapping_2_x2 * mapping_3\n'+spaces2 + D2_5 += 'mapping_x1x3 = mapping_1_x1 * mapping_2 * mapping_3_x3\n'+spaces2 + D2_5 += 'mapping_x2x2 = mapping_1 * mapping_2_x2x2 * mapping_3\n'+spaces2 + D2_5 += 'mapping_x2x3 = mapping_1 * mapping_2_x2 * mapping_3_x3\n'+spaces2 + D2_5 += 'mapping_x3x3 = mapping_1 * mapping_2 * mapping_3_x3x3\n' + D2_6 = '' + for symbol in ('x', 'y', 'z'): + for d1 in range(1, 4): + for d2 in range(1, 4): + if d2 >= d1: + D2_6 += f'{spaces2}{symbol}_x{d1}x{d2} += mapping_x{d1}x{d2} * coeff_{symbol}\n' + D2_6 += '\n' + else: + D2_1 = '' + D2_2 = '' + D2_3 = '' + D2_4 = '' + D2_5 = '' + D2_6 = '' + + body = code_body.format(LOCAL_SPAN = LOCAL_SPAN, + KEYS = KEYS, + A1 = A1, + TEMPS = TEMPS, + COUPLING_TERMS = COUPLING_TERMS, + F_COEFFS_ZEROS = F_COEFFS_ZEROS, + F_SPAN_1 = F_SPAN_1, + F_SPAN_2 = F_SPAN_2, + F_SPAN_3 = F_SPAN_3, + F_COEFFS = F_COEFFS, + F_INIT = F_INIT, + F_ASSIGN_LOOP = F_ASSIGN_LOOP, + D2_1 = D2_1, + D2_2 = D2_2, + D2_3 = D2_3, + D2_4 = D2_4, + D2_5 = D2_5, + D2_6 = D2_6) + + #------------------------- MAKE LOOP ------------------------- + assembly_code = head + body + loop_str = '' + + for block in blocks: + u_i = block[0].indices[0] if nu > 1 else 0 + v_j = block[1].indices[0] if nv > 1 else 0 + + A1 = a1_str.format(u_i=u_i, v_j=v_j) + A2 = a2_str.format(u_i=u_i, v_j=v_j) + A3 = a3_str.format(u_i=u_i, v_j=v_j) + TEST_TRIAL_2 = tt2_str.format(u_i=u_i, v_j=v_j) + TEST_TRIAL_3 = tt3_str.format(u_i=u_i, v_j=v_j) + SPAN_2 = span_2_str.format(u_i=u_i, v_j=v_j) + SPAN_3 = span_3_str.format(u_i=u_i, v_j=v_j) + GLOBAL_SPAN_2 = global_span_2_str.format(u_i=u_i, v_j=v_j) + GLOBAL_SPAN_3 = global_span_3_str.format(u_i=u_i, v_j=v_j) + KEYS_3 = keys_3_str.format(u_i=u_i, v_j=v_j) + KEYS_2 = keys_2_str.format(u_i=u_i, v_j=v_j) + COUPLING_TERMS = ct_str.format(u_i=u_i, v_j=v_j) + + TEST_V_P1, TEST_V_P2, TEST_V_P3 = test_v_p[v_j] + TRIAL_U_P1, TRIAL_U_P2, TRIAL_U_P3 = trial_u_p[u_i] + MAX_P1 = max(TEST_V_P1, TRIAL_U_P1) + MAX_P2 = max(TEST_V_P2, TRIAL_U_P2) + MAX_P3 = max(TEST_V_P3, TRIAL_U_P3) + NEXPR = len(ordered_stmts[block]) + + keys1 = keys_1[block] + TEST_TRIAL_1 = tt1_str.format(u_i=u_i, v_j=v_j) + A2_TEMP = " + ".join([f"{TEST_TRIAL_1}[k_1, q_1, i_1, j_1, {keys1[e][0]}, {keys1[e][1]}] * {A2}[{e},:,:,:,:]" for e in range(NEXPR)]) + + I_1 = f'int(floor(i_1/{test_mult[0]})*{trial_mult[0]})' if max(test_mult[0], trial_mult[0]) > 1 else 'i_1' + I_2 = f'int(floor(i_2/{test_mult[1]})*{trial_mult[1]})' if max(test_mult[1], trial_mult[1]) > 1 else 'i_2' + I_3 = f'int(floor(i_3/{test_mult[2]})*{trial_mult[2]})' if max(test_mult[2], trial_mult[2]) > 1 else 'i_3' + #MAX_P1 = max(int( ( MAX_P1 + np.floor(MAX_P1 / test_mult[0]) * trial_mult[0] ) / 2 ), MAX_P1) if max(test_mult[0], trial_mult[0]) > 1 else MAX_P1 + #MAX_P2 = max(int( ( MAX_P2 + np.floor(MAX_P2 / test_mult[1]) * trial_mult[1] ) / 2 ), MAX_P2) if max(test_mult[1], trial_mult[1]) > 1 else MAX_P2 + #MAX_P3 = max(int( ( MAX_P3 + np.floor(MAX_P3 / test_mult[2]) * trial_mult[2] ) / 2 ), MAX_P3) if max(test_mult[2], trial_mult[2]) > 1 else MAX_P3 + n_cols_x1 = max( int(MAX_P1 + 1 + np.floor(MAX_P1 / test_mult[0]) * trial_mult[0]), 2*MAX_P1+1 ) + n_cols_x2 = max( int(MAX_P2 + 1 + np.floor(MAX_P2 / test_mult[1]) * trial_mult[1]), 2*MAX_P2+1 ) + n_cols_x3 = max( int(MAX_P3 + 1 + np.floor(MAX_P3 / test_mult[2]) * trial_mult[2]), 2*MAX_P3+1 ) + MAX_P1 = n_cols_x1 - MAX_P1 - 1 + MAX_P2 = n_cols_x2 - MAX_P2 - 1 + MAX_P3 = n_cols_x3 - MAX_P3 - 1 + + loop = code_loop.format(A1 = A1, + A2 = A2, + A3 = A3, + TEST_TRIAL_2 = TEST_TRIAL_2, + TEST_TRIAL_3 = TEST_TRIAL_3, + SPAN_2 = SPAN_2, + SPAN_3 = SPAN_3, + GLOBAL_SPAN_2 = GLOBAL_SPAN_2, + GLOBAL_SPAN_3 = GLOBAL_SPAN_3, + KEYS_2 = KEYS_2, + KEYS_3 = KEYS_3, + COUPLING_TERMS = COUPLING_TERMS, + TEST_V_P1 = TEST_V_P1, + TEST_V_P2 = TEST_V_P2, + TEST_V_P3 = TEST_V_P3, + TRIAL_U_P1 = TRIAL_U_P1, + TRIAL_U_P2 = TRIAL_U_P2, + TRIAL_U_P3 = TRIAL_U_P3, + MAX_P1 = MAX_P1, + MAX_P2 = MAX_P2, + MAX_P3 = MAX_P3, + NEXPR = NEXPR, + A2_TEMP = A2_TEMP, + I_1 = I_1, + I_2 = I_2, + I_3 = I_3) + + loop_str += loop + + assembly_code += loop_str + assembly_code += '\n return\n' + + #------------------------- MAKE FILE ------------------------- + import os + if not os.path.isdir('__psydac__'): + os.makedirs('__psydac__') + + # Root process writes the assembly code to a file + if comm is None or comm.rank == 0: + filename = f'__psydac__/assemble_{file_id}.py' + f = open(filename, 'w') + f.writelines(assembly_code) + f.close() + + # Parallel case: wait for the file to be closed before proceeding + if comm is not None and comm.size > 1: + _ = comm.bcast(None, root=0) + + return file_id + + #-------------------------------------------------------------------------- + def read_BilinearForm(self): + """ + Part of the sum factorization algorithm implementation. + Used at the beginning of construct_arguments_generate_assembly_file(). + It's output determines both the design of the assembly function, and the arguments passed to it. + + Returns + ------- + + temps : tuple + tuple of Assign objects. Often times usable building blocks of complicated coupling terms. + + ordered_stmts : dict + assigns each block (trial&test component combination) a list of coupling term assignment + + ordered_sub_exprs_keys : dict + relates each coupling term assignment of ordered_stmts a partial derivative combination + + mapping_option : str | None + 'Bspline' if a spline mapping is involved, None if an analytical or no mapping is involved + + field_derivatives : dict + contains information regarding appearing free FemFields and appearing partial derivatives of those + + g_mat_information_false : list + possibly wrong list of non-zero blocks + + g_mat_information_true : list + correct list of non-zero blocks + + max_logical_derivative : int + maximum appearing partial derivative (in any fixed direction) + + """ + + a = self.expr + domain = a.domain + + # Because an analytical mapping only changes the expression, only the case of a Bspline mapping has to be treated + # entirely different + mapping_option = 'Bspline' if isinstance(self._mapping, SplineMapping) else None + + # The following are tuples consisting of test, trial and free FemField functions appearing, e.g. + # u, v, F1, F2 = elements_of(V, names='u, v, F1, F2) + # a = BilinearForm((u, v), integral(domain, dot(u, F1) * dot(v, F2))) + # tests = (v, ), trials = (u, ) fields = (F1, F2) - Note: The order of F1 & F2 is apparently random and changes from time to time! + # tuple entries are either sympde.topology.space.ScalarFunction or sympde.topology.space.VectorFunction objects + tests = a.test_functions + trials = a.trial_functions + fields = a.fields + + # A sympde.expr.evaluation.DomainExpression object + # TODO [YG 31.07.2025]: Why not using self.terminal_expr[0] instead? + texpr = TerminalExpr(a, domain)[0] + + # We extract all appearing components of test, trial and free FemFields, as well as appearing partial derivatives of these. + # e.g. atoms = [F1[1], F2[1], v[0], u[0], F1[2], F2[2], F1[0], v[1], v[2], F2[0], u[1], u[2]] + # for a bilinear form, without derivatives, involving two vector valued Fem fields F1 & F2 and vector valued test & trial functions v and u + atoms_types = (ScalarFunction, VectorFunction, IndexedVectorFunction) + atoms = _atomic(texpr, cls=atoms_types+_logical_partial_derivatives) + + # Preparing to sort all atoms into test_, trial_ and field_atoms + test_atoms = {} + for v in tests: + if isinstance(v, VectorFunction): + for i in range(domain.dim): + test_atoms[v[i]] = [] + else: + test_atoms[v] = [] + + trial_atoms = {} + for u in trials: + if isinstance(u, VectorFunction): + for i in range(domain.dim): + trial_atoms[u[i]] = [] + else: + trial_atoms[u] = [] + + field_atoms = {} + for f in fields: + if isinstance(f, VectorFunction): + for i in range(domain.dim): + field_atoms[f[i]] = [] + else: + field_atoms[f] = [] + + # atoms can consist of scalar functions (u, v), partial derivatives of scalar functions (dx1(u), dx3(v), ...), + # components of vector valued functions (u[0], v[1], ...), partial derivatives of components of vector valued functions + # (dx1(u[0]), dx3(v[1]), ...), and the same thing but for free FemFields. + # With + # get_atom_logical_derivatives(atom) + # we obtain the component without partial derivatives (u -> u ; dx1(u) -> u ; dx2(v[2]) -> v[2] ; ...) + # This way we can gather subexpressions belonging to the same block + for atom in atoms: + a = get_atom_logical_derivatives(atom) + # IF: NOT Indexed Mapping AND NOT VectorFunction + # I guess: <=> IF ScalarFunction + if not ((isinstance(a, Indexed) and isinstance(a.base, Mapping)) or (isinstance(a, IndexedVectorFunction))): + if a in tests: + # tests is a tuple, e.g. (v, ), hence tests[0] = v + test_atoms[tests[0]].append(atom) + elif a in trials: + trial_atoms[trials[0]].append(atom) + elif a in fields: + # while there can only be one trial and one test function, there can be multiple free FemFields. + for f in field_atoms: + if f == a: + field_atoms[f].append(atom) + else: + raise NotImplementedError(f"atoms of type {str(atom)} are not supported") + # IF VectorFunction + elif isinstance(a, IndexedVectorFunction): + # .base returns ... the base of a VectorFunction! E.g., u[2] -> u, v[0] -> v + if a.base in tests: + for vi in test_atoms: + if vi == a: + test_atoms[vi].append(atom) + break + elif a.base in trials: + for ui in trial_atoms: + if ui == a: + trial_atoms[ui].append(atom) + break + elif a.base in fields: + for fi in field_atoms: + if fi == a: + field_atoms[fi].append(atom) + break + else: + raise NotImplementedError(f"atoms of type {str(atom)} are not supported") + + # ----- Julian O. 11.06.25 ----- + # Regarding the code that follows: + # When dealing with a DiscreteBilinearForm depending on two or more free FemFields, + # the order of the dictionary `field_derivatives` must be the same as the order + # of the free FemFields in `self._free_args`. + # For some reason, the order of all appearing "atoms" in a BilinearForm (trial function, test function, free fields, .?.) + # as obtained in the __init__ of AST + # atoms = terminal_expr.expr.atoms(ScalarFunction, VectorFunction) + # is random and changes from code execution to code execution. + # This order of atoms however determines the order of the free FemFields appearing in `self._free_args`. + # In particular, this order only sometimes matches the order of `field_derivatives`, which results in wrong matrices. + # + # Below is the old version of the code that follows: + #field_derivatives = {} + #for key in field_atoms: + # sym_key = SymbolicExpr(key) + # field_derivatives[sym_key] = {} + # for f in field_atoms[key]: + # field_derivatives[sym_key][SymbolicExpr(f)] = get_index_logical_derivatives(f) + # ------------------------------ + + # For the computation of the coupling terms, among other we need to organize information + # related to free FemFields. For now, we have the dictionary field_atoms, whose keys are + # components of appearing fields, and whose values are appearing partial derivatives of these, e.g., + # field_atoms = {'F1[0]':[dx1(F1[0]), ], 'F1[1]':[dx2(F1[1]), ], 'F1[2]':[dx3(F1[2]), ], 'F2':[F2, ]} + # + # We now create the dictionary field_derivatives. + # It's keys are SymbolicExpr of the previous keys (F1[0] -> F1_0, F1[1] -> F1_1, F1[2] -> F1_2, F2 -> F2) + # and its values are again dictionaries, whose keys are symbolic expressions of the appearing partial derivatives, e.g. + # dx1(F1[0]) -> F1_0_x1, dx2(F1[1]) -> F1_1_x2, dx3(F1[2]) -> F1_2_x3, F2 -> F2, + # and whose values are dictionaries that store the respective derivative information. + # Consider for example the BilinearForm (u, v) \mapsto integral(domain, dot(u, grad(Fs)) * dot(v, grad(Fs2)): + # The corresponding field_derivatives dict will be + # {Fs: {Fs_x3: {'x1': 0, 'x2': 0, 'x3': 1}, Fs_x2: {'x1': 0, 'x2': 1, 'x3': 0}, Fs_x1: {'x1': 1, 'x2': 0, 'x3': 0}}, Fs2: {Fs2_x3: {'x1': 0, 'x2': 0, 'x3': 1}, Fs2_x2: {'x1': 0, 'x2': 1, 'x3': 0}, Fs2_x1: {'x1': 1, 'x2': 0, 'x3': 0}}} + + # Amount of free FemFields (NOT counting each component individually) + n_free_fields = len(self._free_args) + field_derivatives = {} + # The keys in field_derivatives will be in the same order as the fields appearing in self._free_args + for n in range(n_free_fields): + # The key might be F1[0], but we want to check whether F1 == self._free_args[0], and ... + for key in field_atoms: + # ... field_name does exactly that + field_name = str(key.base) if hasattr(key, 'base') else str(key) + if field_name == self._free_args[n]: + # SymbolicExpr transforms something like F1[0] into F1_0 (part of the name of a variable in the assembly code later) + sym_key = SymbolicExpr(key) + field_derivatives[sym_key] = {} + for f in field_atoms[key]: + # And similarly f, which might look like dx1(F1[0]), will be transformed to F1_0_x2 + # while get_index_logical_derivatives(dx1(F1[0])) = {'x1': 1, 'x2': 0, 'x3': 0} + field_derivatives[sym_key][SymbolicExpr(f)] = get_index_logical_derivatives(f) + + # This part was proposed by Said at some point + #syme = False + #if syme: + # from symengine import sympify as syme_sympify + # sym_test_atoms = {k:[syme_sympify(SymbolicExpr(ai)) for ai in a] for k,a in test_atoms.items()} + # sym_trial_atoms = {k:[syme_sympify(SymbolicExpr(ai)) for ai in a] for k,a in trial_atoms.items()} + # sym_expr = syme_sympify(SymbolicExpr(texpr.expr)) + #else: + # sym_test_atoms = {k:[SymbolicExpr(ai) for ai in a] for k,a in test_atoms.items()} + # sym_trial_atoms = {k:[SymbolicExpr(ai) for ai in a] for k,a in trial_atoms.items()} + # sym_expr = SymbolicExpr(texpr.expr) + + # test_atoms is a dict whose values are components of the test function and whose values + # are arrays with appearing partial derivatives of those components. + # sym_test_atoms has the same structure, but replaces the appearing partial derivatives with + # symbolic expressions of those partial derivatives. E.g., + # test_atoms: {v2[0]: [dx3(v2[0]), dx2(v2[0])], v2[1]: [dx1(v2[1]), dx3(v2[1])], v2[2]: [dx2(v2[2]), dx1(v2[2])]} + # sym_test_atoms: {v2[0]: [v2_0_x3, v2_0_x2], v2[1]: [v2_1_x1, v2_1_x3], v2[2]: [v2_2_x2, v2_2_x1]} + # In the following, we will gather all (coupling) terms of a specific combination of a sym_test_atom with a sym_trial_atom in sym_expr + sym_test_atoms = {k:[SymbolicExpr(ai) for ai in a] for k,a in test_atoms.items()} + sym_trial_atoms = {k:[SymbolicExpr(ai) for ai in a] for k,a in trial_atoms.items()} + sym_expr = SymbolicExpr(texpr.expr) + + # ----- temps, rhs ----- + + trials_subs = {ui:0 for u in sym_trial_atoms for ui in sym_trial_atoms[u]} + tests_subs = {vi:0 for v in sym_test_atoms for vi in sym_test_atoms[v]} + sub_exprs = {} + + # This is where the real magic happens: The at times extremely long and complicated SymbolicExpr sym_expr + # 0. is brought into a more readable form (sub_exprs) & + # 1. gets split into many small parts (temps), that often times appear in multiple sub_exprs, + # but now only have to be computed once, e.g. (temp_0, -F2_1*F1_1) & + # 2. those temporaries get assigned to coupling terms (rhs), i.e.: + # The coupling term corresponding to the sub-expr dx1(u[0])*dx3(v[1]) might be -temp_7*(temp_22*temp_27 + temp_33*temp_35 + temp_36*temp_37) + for u in sym_trial_atoms: + for v in sym_test_atoms: + if isinstance(u, IndexedVectorFunction) and isinstance(v, IndexedVectorFunction): + sub_expr = sym_expr[v.indices[0], u.indices[0]] + elif isinstance(u, ScalarFunction) and isinstance(v, ScalarFunction): + sub_expr = sym_expr + elif isinstance(u, ScalarFunction) and isinstance(v, IndexedVectorFunction): + sub_expr = sym_expr[v.indices[0]] + elif isinstance(u, IndexedVectorFunction) and isinstance(v, ScalarFunction): + sub_expr = sym_expr[u.indices[0]] + for ui,sui in zip(trial_atoms[u], sym_trial_atoms[u]): + trcp = trials_subs.copy() + trcp[sui] = 1 + newsub_expr = sub_expr.subs(trcp) + for vi,svi in zip(test_atoms[v],sym_test_atoms[v]): + tcp = tests_subs.copy() + tcp[svi] = 1 + expr = newsub_expr.subs(tcp) + if not expr.is_zero: + sub_exprs[ui,vi] = sympify(expr) + + temps, rhs = cse_main.cse(sub_exprs.values(), symbols=cse_main.numbered_symbols(prefix=f'temp_')) + + # ---------------------- + + # Finally, temps and rhs must be brought into a form that can be included in the assembly code, e.g. + # temp_0 = x_x1*y_x2 + # temp_1 = x_x2*z_x1 + # temp_2 = y_x1*z_x2 + # ... + # coupling_terms_u_v[k_2, q_2, k_3, q_3, 0] = temp_7*(temp_10**2*temp_9 + temp_11**2*temp_9 + temp_8**2*temp_9) + # coupling_terms_u_v[k_2, q_2, k_3, q_3, 1] = temp_18 + # coupling_terms_u_v[k_2, q_2, k_3, q_3, 2] = temp_22 + # ... + + # See above example: In our implementation of the sum factorization algorithm, we precompute arrays + # for each quadrature point in x1 direction, meaning that those arrays contain values depending on + # elements and quadrature points in x2 and x3 direction (k_2, k_3 & q_2 & q_3) + element_indices = [Symbol('k_{}'.format(i)) for i in range(2,4)] + quadrature_indices = [Symbol('q_{}'.format(i)) for i in range(2,4)] + # indices = (k_2, q_2, k_3, q_3) + indices = tuple(j for i in zip(element_indices, quadrature_indices) for j in i) + + # From the sub_exprs dictionary, we read all the appearing trial and test component combinations (blocks) that + # add a non-zero contribution to the matrix + ordered_stmts = {} + ordered_sub_exprs_keys = {} + for key in sub_exprs.keys(): + u_i, v_j = [get_atom_logical_derivatives(atom) for atom in key] + ordered_stmts[u_i, v_j] = [] + ordered_sub_exprs_keys[u_i, v_j] = [] + blocks = ordered_stmts.keys() + + block_list = list(blocks) + trial_components = [block[0] for block in block_list] + test_components = [block[1] for block in block_list] + nu = len(set(trial_components)) + nv = len(set(test_components)) + + expr = self.kernel_expr.expr + + # We store the maximum partial derivative (for a fixed direction), not including pertial derivatives + # appearing in mapping related terms (i.e., a BilinearForm on a mapped domain will have max_logical_derivative = 0 + # even though derivatives of the (spline) mapping appear in the coupling terms). + if isinstance(expr, (ImmutableDenseMatrix, Matrix)): + shape = expr.shape + logical_max_derivatives = [] + for k1 in range(shape[0]): + for k2 in range(shape[1]): + logical_max_derivatives.append(get_max_logical_partial_derivatives(expr[k1,k2])) + max_logical_derivative = max([max([value for value in dic.values()]) for dic in logical_max_derivatives]) + else: + max_logical_derivative = max([value for value in get_max_logical_partial_derivatives(expr).values()]) + + # See comment underneath this code block for more details. + # There was a test case, in which the amount of generated StencilMatrices (one for each appearing block, + # i.e., one for each trial&test component combination for which a non-zero coupling term exists) + # was larger than the amount true amount of needed StencilMatrices. + # That discrepancy appears when expr[block].is_zero wrongly does not detect that a block is zero, + # whereas the corresponding block does rightfully not appear in block_list! + if isinstance(expr, (ImmutableDenseMatrix, Matrix)): # only relevenat if either trial or test function is vector valued + g_mat_information_false = [] + shape = expr.shape + for k1 in range(shape[0]): + for k2 in range(shape[1]): + if not expr[k1,k2].is_zero: # although it might actually be zero! + if (nu == 1) and (nv > 1): + g_mat_information_false.append((k2,k1)) + else: + g_mat_information_false.append((k1,k2)) + if nu == 1: + g_mat_information_true = [(0, get_atom_logical_derivatives(block[1]).indices[0]) for block in block_list] + elif nv == 1: + g_mat_information_true = [(get_atom_logical_derivatives(block[0]).indices[0], 0) for block in block_list] + else: + g_mat_information_true = [(get_atom_logical_derivatives(block[0]).indices[0], get_atom_logical_derivatives(block[1]).indices[0]) for block in block_list] + else: + g_mat_information_false = [] + g_mat_information_true = [] + + # Julian O. 17.06.25: Back when I added this unreadable comment below I forgot to write a test for this problem. + # Eventually it might be interesting to remove everything related to `g_mat_information_false/true` + # and see where errors occur. + # + #1, 1: expr[1,1] = F0*sqrt(x1**2*(x1*cos(2*pi*x3) + 2)**2*(sin(pi*x2)**2 + cos(pi*x2)**2)**2*(sin(2*pi*x3)**2 + cos(2*pi*x3)**2)**2)*(pi*(x1*cos(2*pi*x3) + 2)* + # (-2*pi*x1*sin(pi*x2)*sin(2*pi*x3)*dx1(v1[1]) - sin(pi*x2)*cos(2*pi*x3)*dx3(v1[1]))*cos(pi*x2)*w2[1] - pi*(x1*cos(2*pi*x3) + 2)*(-2*pi*x1*sin(2*pi*x3)*cos(pi*x2)*dx1(v1[1]) - + # cos(pi*x2)*cos(2*pi*x3)*dx3(v1[1]))*sin(pi*x2)*w2[1])/(2*pi**2*x1**2*(x1*cos(2*pi*x3) + 2)**2*(sin(pi*x2)**2 + cos(pi*x2)**2)**2*(sin(2*pi*x3)**2 + cos(2*pi*x3)**2)**2) + # = 0 - but is not yet detected as 0! Hence a matrix is generated, that later is not required! + # + + # Here we create a template for the names of the coupling terms arrays, + # depending on whether or not trial and test function are scalar or vector valued + if nv > 1: + ct_str = 'coupling_terms_u_{u_i}_v_{v_j}' if nu > 1 else 'coupling_terms_u_v_{v_j}' + else: + ct_str = 'coupling_terms_u_{u_i}_v' if nu > 1 else 'coupling_terms_u_v' + + # Now we format this template based on the appearing blocks (combinations of trial and test function components) + # and transform those formatted strings into IndexedBase objects + lhs = {} + for block in blocks: + u_i = get_atom_logical_derivatives(block[0]).indices[0] if nu > 1 else 0 + v_j = get_atom_logical_derivatives(block[1]).indices[0] if nv > 1 else 0 + ct = ct_str.format(u_i=u_i, v_j=v_j) + lhs[block] = IndexedBase(f'{ct}') + + # lhs[block] will look w.g. like this coupling_terms_u_v (u, v scalar). + # Now, we add to that [k_2, q_2, k_3, q_3, count], where count enumerates the sub expressions belonging to the same block + # sub expressions corresponding to the block (u[0], v[1]) might be: (u[0], v[1]), (dx1(u[0]), v[1]), (dx2(u[0]), v[1]), ... + # and then assign the corresponding rhs, e.g. temp_7*(temp_10**2*temp_9 + temp_11**2*temp_9 + temp_8**2*temp_9), to obtain: + # coupling_terms_u_v[k_2, q_2, k_3, q_3, 4] = temp_7*(temp_10**2*temp_9 + temp_11**2*temp_9 + temp_8**2*temp_9) + counts = {block:0 for block in blocks} + for r,key in zip(rhs, sub_exprs.keys()): + u_i, v_j = [get_atom_logical_derivatives(atom) for atom in key] + count = counts[u_i, v_j] + counts[u_i, v_j] += 1 + ordered_stmts[u_i, v_j].append(Assign(lhs[u_i, v_j][(*indices, count)], r)) + ordered_sub_exprs_keys[u_i, v_j].append(key) + # ordered_stmts is a dict whose keys are combinations of trial and test functions components (e.g. u[0], v[1]), + # and whose values are a list of coupling term assignments corresponding to this block, e.g. + # (v1[0], v2[0]): [coupling_terms_u_0_v_0[k_2, q_2, k_3, q_3, 0] := -1, coupling_terms_u_0_v_0[k_2, q_2, k_3, q_3, 1] := 1] + # + # The information regarding which partial derivative combination belongs to which coupling term is stored in ordered_sub_exprs_keys. + # This dict has the same keys, but instead of coupling term assignments as values, list of tuples of partial derivative combinations are stored. + + # temps, which previously consisted of tuples like this one: (temp_0, -F2_1*F1_1), + # will now be a tuple consisting of assignments, e.g. (temp_0 := -F2_1*F1_1, ...) + temps = tuple(Assign(a,b) for a,b in temps) + + return temps, ordered_stmts, ordered_sub_exprs_keys, mapping_option, field_derivatives, g_mat_information_false, g_mat_information_true, max_logical_derivative + + #-------------------------------------------------------------------------- + def construct_arguments_generate_assembly_file(self): + """ + Collect the arguments used in the assembly method, and generate and possibly pyccelize the assembly function. + + Used only when sum factorization is enabled, else the method construct_arguments is called. + + Returns + ------- + args: tuple + The arguments passed to the assembly method. + + threads_args: None + None as openMP parallelization is not supported by this implementation. + + """ + temps, ordered_stmts, ordered_sub_exprs_keys, mapping_option, field_derivatives, g_mat_information_false, g_mat_information_true, max_logical_derivative = self.read_BilinearForm() + + # Each block corresponds to a combination of trial and test function components, and thus indeed to a "block" in the matrix. + # Not all possible combination have to exist, e.g., + # given a function space of vector valued functions V (3d) and a bilinear form a: VxV -> R, a(u, v) = (u, v)_L^2(Omega) + # there will be only 3 blocks on a logical domain (u[0]&v[0], u[1]&v[1], u[2]&v[2]), + # but up to 9 blocks on a mapped domain (e.g. u[0]&v[1], ...) + blocks = ordered_stmts.keys() + block_list = list(blocks) + trial_components = [block[0] for block in block_list] + test_components = [block[1] for block in block_list] + # dim = 1 corresponds to a scalar valued function, dim = 3 to a vector valued function + trial_dim = len(set(trial_components)) + test_dim = len(set(test_components)) + + # A reminder that this implementation only supports bilinear forms on 3d domains. + d = 3 + assert d == 3 + + # Rename - also: establish that throughout "u" corresponds to the trial function, whereas "v" corresponds to the test function + nu = trial_dim # dim of trial function; 1 (scalar) or 3 (vector) + nv = test_dim # dim of test function ; 1 (scalar) or 3 (vector) + + # Obtain the most basic information: function values, degrees, spans, ... + test_basis, test_degrees, spans, pads, test_mult = construct_test_space_arguments(self.test_basis) + trial_basis, trial_degrees, pads, trial_mult = construct_trial_space_arguments(self.trial_basis) + n_elements, quads, quad_degrees = construct_quad_grids_arguments(self.grid[0], use_weights=False) + + #! pads is being overwritten. That is because already somewhere else (__init__ of StencilMatrix via self.allocate_matrices) + # do we assert that domain and codomain (trial and test) pads coincide! + # That is not strictly necessary as Valentin at some point proved in one of his branches, but currently not implemented as + # not required. + + #! the above pads variable is multiplied by the multiplicity vector! For the remaining implementation, we need + # the pads vector un-multiplied, as obtained by : + pads = self.test_basis.space.coeff_space.pads + + # quad_degrees is the amount of quadrature points per element in each direction + # Clearly, this amount must coincide with the amount of basis function values stored per element in test_basis and trial_basis + n_element_1, n_element_2, n_element_3 = n_elements + k1, k2, k3 = quad_degrees + + # We store component wise degree and function values for trial and test function in the dictionaries + # trial_u_p, global_basis_u, test_v_p, global_basis_v + if (nu == 3) and (len(trial_basis) == 3): + # Edge Case: If the trial function space V is a VectorFunctionSpace + # but neither an Hdiv nor an Hcurl space, i.e., + # V = VectorFunctionSpace('V', domain) and not +, kind='hcurl') or +, kind='hdiv') + # then the function values in each of the three directions are identical for each of the three components. + # Hence len(trial_basis) == 3 instead of 9. + # + # global_basis_u is a dict whose values are arrays of function values of one particular trial function component, + # hence for this edge case we simply assign the same array trial_basis to each component + # Same function degree in each direction for each component -> do the same thing with trial_u_p + trial_u_p = {u:trial_degrees for u in range(nu)} + global_basis_u = {u:trial_basis for u in range(nu)} + else: + trial_u_p = {u:trial_degrees[d*u:d*(u+1)] for u in range(nu)} + global_basis_u = {u:trial_basis[d*u:d*(u+1)] for u in range(nu)} + if (nv == 3) and (len(test_basis) == 3): + # See above explanation, which also applies for the spans variable + test_v_p = {v:test_degrees for v in range(nv)} + global_basis_v = {v:test_basis for v in range(nv)} + spans = [*spans, *spans, *spans] + else: + test_v_p = {v:test_degrees[d*v:d*(v+1)] for v in range(nv)} + global_basis_v = {v:test_basis[d*v:d*(v+1)] for v in range(nv)} + + # See other method construct_arguments: + # When self._target is an Interface domain len(self._grid) == 2 + # where grid contains the QuadratureGrid of both sides of the interface + assert len(self.grid) == 1 + if self.mapping: + # We gather mapping related information in the case of a Bspline mapping + # self.mapping == False if either no or an analytical mapping + + map_coeffs = [[e._coeffs._data for e in self.mapping._fields]] + spaces = [self.mapping._fields[0].space] + map_degree = [sp.degree for sp in spaces] + map_span = [[q.spans - s for q,s in zip(sp.get_assembly_grids(*self.nquads), sp.coeff_space.starts)] for sp in spaces] + map_basis = [[q.basis for q in sp.get_assembly_grids(*self.nquads)] for sp in spaces] + points = [g.points for g in self.grid] + weights = [self.mapping.weights_field.coeffs._data] if self.is_rational_mapping else [] + + for i in range(len(self.grid)): + axis = self.grid[i].axis + # See construct_arguments - have not come across an example of when axis was not None! + assert axis is None + + map_degree = flatten(map_degree) + map_span = flatten(map_span) + map_basis = flatten(map_basis) + points = flatten(points) + mapping = [*map_coeffs[0], *weights] + else: + + mapping = [] + map_degree = [] + map_span = [] + map_basis = [] + + #---------- The following part is entirely different from the old construct_arguments method ---------- + + # Each block, say u[0]&v[1], + # consists of possibly many derivative combinations (sub-expressions) of these two components, e.g. + # dx1(u[0])&dx1(v[1]) or dx1(u[0])&dx2(v[1]) (dx1, dx2, dx3 representing respective partial derivatives). + # + # For each block, here still e.g. u[0]&v[1], + # and for each sub-expression, we store corresponding derivative information: + # get_index_logical_derivatives(dx1(u[0])) = {'x1': 1, 'x2': 0, 'x3': 0} + # get_index_logical_derivatives(dx2(v[1])) = {'x1': 0, 'x2': 1, 'x3': 0} + # Each of these 6 dicts has for each block an array of length #sub-expressions (appearing derivative combination) stored + # x2_test_keys[(u[0], v[1])][3] = 2 means, that the fourth sub-expression of block (u[0], v[1]) + # involves a second partial derivative of the test function in x2 direction + x1_trial_keys = {block:[] for block in blocks} + x1_test_keys = {block:[] for block in blocks} + x2_trial_keys = {block:[] for block in blocks} + x2_test_keys = {block:[] for block in blocks} + x3_trial_keys = {block:[] for block in blocks} + x3_test_keys = {block:[] for block in blocks} + + for block in blocks: + # alpha, beta for example being dx1(u[0]), dx2(v[1]) + for alpha, beta in ordered_sub_exprs_keys[block]: + x1_trial_keys[block].append(get_index_logical_derivatives(alpha)['x1']) + x1_test_keys [block].append(get_index_logical_derivatives(beta) ['x1']) + x2_trial_keys[block].append(get_index_logical_derivatives(alpha)['x2']) + x2_test_keys [block].append(get_index_logical_derivatives(beta) ['x2']) + x3_trial_keys[block].append(get_index_logical_derivatives(alpha)['x3']) + x3_test_keys [block].append(get_index_logical_derivatives(beta) ['x3']) + + # See sum factorization paper by Bressan & Takacs: + # coupling_terms, a3 and a2 correspond to A^{>=4}_{x1,x2,x3}, A^{>=3}_{x1,x2} and A^{>=2}_{x1} + # Here, for each block we assign a zero-array of the correct size. + coupling_terms = {} + a3 = {} + a2 = {} + + # For each block, we precompute ~enough~ products of partial derivatives of trial and basis functions in each direction + # These precomputed values will then be read rather than computed in the assembly + test_trial_1s = {} + test_trial_2s = {} + test_trial_3s = {} + + # keys_1/2/3 is a restructuring of the 6 dictionaries created above + keys_1 = {} + keys_2 = {} + keys_3 = {} + + assembly_backend = self.backend + if self._pyccelize_test_trial_computation and assembly_backend['name'] == 'pyccel': + + import os + if not os.path.isdir('__psydac__'): + os.makedirs('__psydac__') + + comm = self.comm + + if comm is not None and comm.size > 1: + if comm.rank == 0: + filename = '__psydac__/test_trial_computation.py' + code = self.test_trial_template + f = open(filename, 'w') + f.writelines(code) + f.close() + else: + filename = '__psydac__/test_trial_computation.py' + code = self.test_trial_template + f = open(filename, 'w') + f.writelines(code) + f.close() + + base_dirpath = os.getcwd() + sys.path.insert(0, base_dirpath) + + package = importlib.import_module(f'__psydac__.test_trial_computation') + kwargs = { + 'language' : 'fortran', + 'compiler_family' : assembly_backend['compiler_family'], + 'flags' : assembly_backend['flags'], + 'openmp' : True if assembly_backend['openmp'] else False, + 'verbose' : False, + 'comm' : self.comm, + } + + test_trial_func = epyccel(package.test_trial_array, **kwargs) + + for block in blocks: + # We translate a block, e.g. (u[0], v[1]) into two integers u_i=0, v_j=1. + # In the case of a scalar function (u, v instead of u[0], u[1], u[2], v[0], v[1], v[2]), store 0. + u_i = block[0].indices[0] if nu > 1 else 0 + v_j = block[1].indices[0] if nv > 1 else 0 + + # keys_2[(u[0], v[1])][3] = (1,2) means that the fourth sub-expression corresponding to the trial-test-function-component-product + # u[0] * v[1] involves a first derivative in x2 direction of the trial function and a second derivative in x2 direction of the test function + keys_1[block] = np.array([(alpha_1, beta_1) for alpha_1, beta_1 in zip(x1_trial_keys[block], x1_test_keys[block])]) + keys_2[block] = np.array([(alpha_2, beta_2) for alpha_2, beta_2 in zip(x2_trial_keys[block], x2_test_keys[block])]) + keys_3[block] = np.array([(alpha_3, beta_3) for alpha_3, beta_3 in zip(x3_trial_keys[block], x3_test_keys[block])]) + + # Those are the function values in each direction of a particular component of the trial/test function + global_basis_u_1, global_basis_u_2, global_basis_u_3 = global_basis_u[u_i] + global_basis_v_1, global_basis_v_2, global_basis_v_3 = global_basis_v[v_j] + + # Those are the Bspline degrees in each direction of a particular component of the trial/test function + trial_u_p1, trial_u_p2, trial_u_p3 = trial_u_p[u_i] + test_v_p1, test_v_p2, test_v_p3 = test_v_p [v_j] + + max_p_2 = max(test_v_p2, trial_u_p2) + max_p_3 = max(test_v_p3, trial_u_p3) + + # That's the amount of subexpressions, i.e., combinations of partial derivatives appearing for a specific combination of + # trial and test function components + n_expr = len(ordered_stmts[block]) + + # To compute enough (possibly too many, but never too few) products of trial and test functions, we read the maximum + # appearing partial derivative (for this specific block, in each direction, for both trial and test function) + max_block_trial_x1_derivative = max(x1_trial_keys[block]) + max_block_trial_x2_derivative = max(x2_trial_keys[block]) + max_block_trial_x3_derivative = max(x3_trial_keys[block]) + max_block_test_x1_derivative = max(x1_test_keys[block]) + max_block_test_x2_derivative = max(x2_test_keys[block]) + max_block_test_x3_derivative = max(x3_test_keys[block]) + + # On each Bspline cell (element / subdomain), there are (test_degree+1)*(trial_degree+1) test & trial function pairs + # of non-zero product. + # Hence, we assign zeros for each element, each quadrature point on the element, each test and trial function combination, + # and each (or even more than required) appearing partial derivative combination of these functions - in each direction + test_trial_1 = np.zeros((n_element_1, k1, test_v_p1 + 1, trial_u_p1 + 1, max_block_trial_x1_derivative+1, max_block_test_x1_derivative+1), dtype='float64') + test_trial_2 = np.zeros((n_element_2, k2, test_v_p2 + 1, trial_u_p2 + 1, max_block_trial_x2_derivative+1, max_block_test_x2_derivative+1), dtype='float64') + test_trial_3 = np.zeros((n_element_3, k3, test_v_p3 + 1, trial_u_p3 + 1, max_block_trial_x3_derivative+1, max_block_test_x3_derivative+1), dtype='float64') + + # And that's how we fill the test_trial arrays + if self._pyccelize_test_trial_computation and assembly_backend['name'] == 'pyccel': + for args in zip(n_elements, + quad_degrees, [test_v_p1, test_v_p2, test_v_p3], [trial_u_p1, trial_u_p2, trial_u_p3], + [global_basis_u_1, global_basis_u_2, global_basis_u_3], [global_basis_v_1, global_basis_v_2, global_basis_v_3], + [max_block_trial_x1_derivative, max_block_trial_x2_derivative, max_block_trial_x3_derivative], [max_block_test_x1_derivative, max_block_test_x2_derivative, max_block_test_x3_derivative], + [test_trial_1, test_trial_2, test_trial_3]): + + args = tuple(np.int64(a) if isinstance(a, int) else a for a in args) + + test_trial_func(*args) + else: + for k_1 in range(n_element_1): + for q_1 in range(k1): + for i_1 in range(test_v_p1 + 1): + for j_1 in range(trial_u_p1 + 1): + trial = global_basis_u_1[k_1, j_1, :, q_1] + test = global_basis_v_1[k_1, i_1, :, q_1] + for alpha_1 in range(max_block_trial_x1_derivative+1): + for beta_1 in range(max_block_test_x1_derivative+1): + test_trial_1[k_1, q_1, i_1, j_1, alpha_1, beta_1] = trial[alpha_1] * test[beta_1] + + for k_2 in range(n_element_2): + for q_2 in range(k2): + for i_2 in range(test_v_p2 + 1): + for j_2 in range(trial_u_p2 + 1): + trial = global_basis_u_2[k_2, j_2, :, q_2] + test = global_basis_v_2[k_2, i_2, :, q_2] + for alpha_2 in range(max_block_trial_x2_derivative+1): + for beta_2 in range(max_block_test_x2_derivative+1): + test_trial_2[k_2, q_2, i_2, j_2, alpha_2, beta_2] = trial[alpha_2] * test[beta_2] + + for k_3 in range(n_element_3): + for q_3 in range(k3): + for i_3 in range(test_v_p3 + 1): + for j_3 in range(trial_u_p3 + 1): + trial = global_basis_u_3[k_3, j_3, :, q_3] + test = global_basis_v_3[k_3, i_3, :, q_3] + for alpha_3 in range(max_block_trial_x3_derivative+1): + for beta_3 in range(max_block_test_x3_derivative+1): + test_trial_3[k_3, q_3, i_3, j_3, alpha_3, beta_3] = trial[alpha_3] * test[beta_3] + + test_trial_1s[block] = test_trial_1 + test_trial_2s[block] = test_trial_2 + test_trial_3s[block] = test_trial_3 + + # Instead of having a different a3, a2 & coupling term array for each sub-expression, we choose to have only one + # such array per block. + # a3 will store line integral values for all combinations of test and trial functions in x3 direction, hence the dimension + # (n_element_3 + test_v_p3 + (mult[2]-1)*(n_element_3-1), 2 * max_p_3 + 1) + # a2 will store surface integral values for all combinations of test and trial functions in x2 and x3 direction, hence the dimension ... + # coupling_terms stores point values of the coupling terms at all quadrature points + # but only in x2 and x3 direction, because we only "precompute" this array for a fixed quadrature point in x1 direction + + # a3[block] size explained: #sub expressions ; #test functions depending on x3 ; #complicated expression for the minimum columns needed + # to store local information correctly. 2*degree+1 in the simplest case. + n_funs_x2 = n_element_2 + test_v_p2 + (test_mult[1]-1)*(n_element_2-1) + n_funs_x3 = n_element_3 + test_v_p3 + (test_mult[2]-1)*(n_element_3-1) + n_cols_x2 = max( int(max_p_2 + 1 + np.floor(max_p_2 / test_mult[1]) * trial_mult[1]), 2*max_p_2+1 ) + n_cols_x3 = max( int(max_p_3 + 1 + np.floor(max_p_3 / test_mult[2]) * trial_mult[2]), 2*max_p_3+1 ) + + a3[block] = np.zeros((n_expr, n_funs_x3, n_cols_x3), dtype='float64') + a2[block] = np.zeros((n_expr, n_funs_x2, n_funs_x3, n_cols_x2, n_cols_x3), dtype='float64') + + coupling_terms[block] = np.zeros((n_element_2, k2, n_element_3, k3, n_expr), dtype='float64') + + # We gather the socalled new args - all other args are being obtained in a similar way using the old assembly implementation + new_args = (*list(test_trial_1s.values()), + *list(test_trial_2s.values()), + *list(test_trial_3s.values()), + *list(a3.values()), + *list(a2.values()), + *list(coupling_terms.values())) + + # This part is a bit shady. + # There has been a case, where my code wasn't running, because one instance of deep-(Psydac/Sympde/Sympy)-code + # correctly understood that a possibly complicated expression (corresponding to a block) in fact evaluates to 0, + # and hence no StencilMatrix for that particular block ever needs to be created - but a different part of + # deep-(Psydac/Sympde/Sympy)-code did not get that simplification right (yet?), and decided that the assembly code + # needs a StencilMatrix as input for this particular block. + # See readBilinearForm for additional information. + # This part of the code filters out unnecessary StencilMatrices, such that only the relevant StencilMatrices + # are being passed to the assembly function + expr = self.kernel_expr.expr + if isinstance(expr, (ImmutableDenseMatrix, Matrix)): + matrices = [] + for i, block in enumerate(g_mat_information_false): + if block in g_mat_information_true: + matrices.append(self._global_matrices[i]) + else: + matrices = self._global_matrices + + # We have gathered all args! + args = (*map_basis, *spans, *map_span, *quads, *map_degree, *n_elements, *quad_degrees, *pads, *mapping, *matrices, + *new_args) + + threads_args = () + + args = tuple(np.int64(a) if isinstance(a, int) else a for a in args) + threads_args = tuple(np.int64(a) if isinstance(a, int) else a for a in threads_args) + + #---------- We now generate the assembly file ---------- + + # file_id is a random string that has been used to name the assembly file + file_id = self.make_file(temps, ordered_stmts, field_derivatives, max_logical_derivative, test_mult, trial_mult, test_v_p, trial_u_p, keys_1, keys_2, keys_3, mapping_option) + + # Store the current directory and add it to the variable `sys.path` + # to imitate Python's import behavior + import os + base_dirpath = os.getcwd() + sys.path.insert(0, base_dirpath) + + # Import the generated assembly function + package = importlib.import_module(f'__psydac__.assemble_{file_id}') + + # The assembly function is the one that has been generated in the make_file method + assembly_function_name = f'assemble_matrix_{file_id}' + assembly_function = getattr(package, assembly_function_name) + + # If the backend is pyccel, we compile the new assembly function + assembly_backend = self.backend + if assembly_backend['name'] == 'pyccel': + kwargs = { + 'language' : 'fortran', # hardcoded for now + 'compiler_family' : assembly_backend['compiler_family'], + 'flags' : assembly_backend['flags'], + 'openmp' : True if assembly_backend['openmp'] else False, + 'verbose' : False, + # 'folder': assembly_backend['folder'], + 'comm' : self.comm, + # 'time_execution': verbose, + # 'verbose': verbose + } + new_func = epyccel(assembly_function, **kwargs) + else: + new_func = assembly_function + + # Use the new assembly function (either compiled or not) + self._func = new_func + + return args, threads_args + + #-------------------------------------------------------------------------- + @property + def test_trial_template(self): + code = '''def test_trial_array(n_element : "int64", + quad_degree : "int64", test_degree : "int64", trial_degree : "int64", + trial_basis : "float64[:,:,:,:]", test_basis : "float64[:,:,:,:]", + max_trial_derivative : "int64", max_test_derivative : "int64", + test_trial : "float64[:,:,:,:,:,:]"): + + for k in range(n_element): + for q in range(quad_degree): + for i in range(test_degree + 1): + for j in range(trial_degree + 1): + trial = trial_basis[k, j, :, q] + test = test_basis [k, i, :, q] + for alpha in range(max_trial_derivative + 1): + for beta in range(max_test_derivative + 1): + test_trial[k, q, i, j, alpha, beta] = trial[alpha] * test[beta] + + return +''' + return code diff --git a/psydac/api/fem_common.py b/psydac/api/fem_common.py new file mode 100644 index 000000000..9c6c0bd6a --- /dev/null +++ b/psydac/api/fem_common.py @@ -0,0 +1,286 @@ +from typing import Iterable + +from sympy import Expr, ImmutableDenseMatrix, Matrix + +import numpy as np + +from sympde.expr.basic import BasicForm +from sympde.expr.evaluation import KernelExpression +from sympde.topology.space import ScalarFunction, VectorFunction, IndexedVectorFunction +from sympde.topology.space import ProductSpace, VectorFunctionSpace +from sympde.topology.datatype import H1SpaceType, L2SpaceType, UndefinedSpaceType +from sympde.topology.derivatives import get_atom_logical_derivatives +from sympde.topology.derivatives import _logical_partial_derivatives + +from psydac.fem.basic import FemSpace +from psydac.linalg.stencil import StencilMatrix, StencilInterfaceMatrix +from psydac.linalg.basic import ComposedLinearOperator +from psydac.api.utilities import flatten +from psydac.api.ast.utilities import math_atoms_as_str, get_max_partial_derivatives + +# TODO [YG 01.08.2025]: Avoid importing anything from psydac.pyccel +from psydac.pyccel.ast.core import _atomic + +__all__ = ( + 'compute_max_nderiv', + 'compute_imports', + 'compute_free_arguments', + 'collect_spaces', + 'compute_diag_len', + 'construct_test_space_arguments', + 'construct_trial_space_arguments', + 'construct_quad_grids_arguments', + 'do_nothing', + 'extract_stencil_mats', + 'reset_arrays', +) + +#============================================================================== +def compute_max_nderiv(kernel_expr: KernelExpression) -> int: + """ + Compute the highest derivative order in the given kernel expression. + + Parameters + ---------- + kernel_expr : KernelExpression (from sympde.expr.evaluation) + + Returns + ------- + nderiv : int + The highest order of derivation in `terminal_expr`. + + """ + assert isinstance(kernel_expr, KernelExpression) + + terminal_expr = kernel_expr.expr + if not isinstance(terminal_expr, (ImmutableDenseMatrix, Matrix)): + terminal_expr = ImmutableDenseMatrix([[terminal_expr]]) + n_rows, n_cols = terminal_expr.shape + + atoms_types = (ScalarFunction, VectorFunction, IndexedVectorFunction) + extended_atoms_types = atoms_types + _logical_partial_derivatives + + nderiv = 0 + for i_row in range(n_rows): + for i_col in range(n_cols): + texpr = terminal_expr[i_row, i_col] + atoms = _atomic(texpr, cls=extended_atoms_types) + Fs = [get_atom_logical_derivatives(a) for a in atoms] + d = get_max_partial_derivatives(texpr, logical=True, F=Fs) + nderiv = max(nderiv, max(d.values())) + + return nderiv + +#============================================================================== +def compute_imports(expr: Expr, + spaces: Iterable[FemSpace], + *, + openmp: bool + ) -> dict[str, list[str]]: + """ + Compute all the imports to be added to the generated Python code. + + Parameters + ---------- + expr : sympy.Expr + The integrand expression of a BilinearForm, LinearForm, or + Functional. This is a pure SymPy expression where SymPDE partial + derivatives have been converted to SymPy symbols. See Notes. + + spaces : iterable of psydac.fem.FemSpace + The discrete spaces which define the finite element representation + of a BilinearForm, LinearForm, or Functional. + + openmp : bool + Whether or not OpenMP pragmas and functions are used in the code. + + Returns + ------- + imports : dict[str, list[str]] + A dictionary whose keys are the names of the Python modules to be + imported, and whose values are the names of the corresponding + objects (variables, functions, classes) to be imported from the + modules. + + Notes + ----- + Assume that we start from an object of type BilinearForm, LinearForm, + or Functional. We take the integrand and expand its vector operations + with TerminalExpr(), then pull back from physical to logical + coordinates with LogicalExpr(), and finally convert the symbolic + partial derivatives with SymbolicExpr(). Where: + - TerminalExpr is defined in sympde.expr.evaluation + - LogicalExpr is defined in sympde.topology.mapping + - SymbolicExpr is defined in sympde.topology.mapping + + The resulting expression `expr` can be passed to this function. + """ + assert isinstance(expr, Expr) + assert all(isinstance(V, FemSpace) for V in spaces) + assert isinstance(openmp, bool) + + # Determine the type of scalar quantities to be managed in the code + dtypes = [getattr(V.symbolic_space, 'codomain_type', 'real') for V in spaces] + assert all(t in ['complex', 'real'] for t in dtypes) + dtype = 'complex' if 'complex' in dtypes else 'real' + + # TODO uncomment this line when we have a SesquilinearForm defined in SymPDE + #assert isinstance(expr, SesquilinearForm) + + #... Compute the imports + math_library = 'cmath' if dtype=='complex' else 'math' # Function names are the same + math_imports = math_atoms_as_str(expr, 'math') + numpy_imports = ['array', 'zeros', 'zeros_like', 'floor'] + + imports = {'numpy': numpy_imports} + if math_imports: + imports[math_library] = math_imports + if openmp: + imports['pyccel.stdlib.internal.openmp'] = ['omp_get_thread_num'] + #... + + return imports + +#============================================================================== +def compute_free_arguments(expr: BasicForm, kernel_expr: KernelExpression) -> tuple[str]: + """ + The string representation (i.e. the names) of the free arguments in + the given BilinearForm, LinearForm, or Functional. + + Parameters + ---------- + expr : BilinearForm | LinearForm | Functional + The expression of which we want to compute the free arguments. + + kernel_expr : sympde.expr.evaluation.KernelExpression + The atomic representation of the form, which is obtained after using + LogicalExpr (if there is a mapping) and TerminalExpr on `expr`. + + Returns + ------- + tuple[str] + The string representation (i.e. the names) of the free arguments in + the given BilinearForm, LinearForm, or Functional. + + """ + assert isinstance(expr, BasicForm) + assert isinstance(kernel_expr, KernelExpression) + + free_args_dict = expr.get_free_variables() + free_args_str = tuple(str(a) for a in free_args_dict) + + return free_args_str + +#============================================================================== +def collect_spaces(space, *args): + """ + This function collect the arguments used in the assembly function + + Parameters + ---------- + space: + the symbolic space + + args : + list of discrete space components like basis values, spans, ... + + Returns + ------- + args : + list of discrete space components elements used in the asembly + + """ + + if isinstance(space, ProductSpace): + spaces = space.spaces + indices = [] + i = 0 + for space in spaces: + if isinstance(space, VectorFunctionSpace): + if isinstance(space.kind, (H1SpaceType, L2SpaceType, UndefinedSpaceType)): + indices.append(i) + else: + indices += [i+j for j in range(space.ldim)] + i = i + space.ldim + else: + indices.append(i) + i = i + 1 + args = [[e[i] for i in indices] for e in args] + + elif isinstance(space, VectorFunctionSpace): + if isinstance(space.kind, (H1SpaceType, L2SpaceType, UndefinedSpaceType)): + args = [[e[0]] for e in args] + + return args + +#============================================================================== +def compute_diag_len(p, md, mc): + n = ((np.ceil((p+1)/mc)-1)*md).astype('int') + n = n-np.minimum(0, n-p)+p+1 + return n.astype('int') + +#============================================================================== +def construct_test_space_arguments(basis_values): + space = basis_values.space + test_basis = basis_values.basis + spans = basis_values.spans + test_degrees = space.degree + pads = space.pads + multiplicity = space.multiplicity + + test_basis, test_degrees, spans = collect_spaces(space.symbolic_space, test_basis, test_degrees, spans) + + test_basis = flatten(test_basis) + test_degrees = flatten(test_degrees) + spans = flatten(spans) + pads = flatten(pads) + multiplicity = flatten(multiplicity) + pads = [p*m for p,m in zip(pads, multiplicity)] + return test_basis, test_degrees, spans, pads, multiplicity + +def construct_trial_space_arguments(basis_values): + space = basis_values.space + trial_basis = basis_values.basis + trial_degrees = space.degree + pads = space.pads + multiplicity = space.multiplicity + trial_basis, trial_degrees = collect_spaces(space.symbolic_space, trial_basis, trial_degrees) + + trial_basis = flatten(trial_basis) + trial_degrees = flatten(trial_degrees) + pads = flatten(pads) + multiplicity = flatten(multiplicity) + pads = [p*m for p,m in zip(pads, multiplicity)] + return trial_basis, trial_degrees, pads, multiplicity + +#============================================================================== +def construct_quad_grids_arguments(grid, use_weights=True): + points = grid.points + if use_weights: + weights = grid.weights + quads = flatten(list(zip(points, weights))) + else: + quads = flatten(list(zip(points))) + + nquads = flatten(grid.nquads) + n_elements = grid.n_elements + return n_elements, quads, nquads + +#============================================================================== +def do_nothing(*args): + return 0 + +#============================================================================== +def extract_stencil_mats(mats): + new_mats = [] + for M in mats: + if isinstance(M, (StencilInterfaceMatrix, StencilMatrix)): + new_mats.append(M) + elif isinstance(M, ComposedLinearOperator): + new_mats += [i for i in M.multiplicants if isinstance(i, (StencilInterfaceMatrix, StencilMatrix))] + return new_mats + +#============================================================================== +def reset_arrays(*args): + for a in args: + a[:]= 0.j if a.dtype==complex else 0. diff --git a/psydac/api/fem_sum_form.py b/psydac/api/fem_sum_form.py new file mode 100644 index 000000000..7d387b733 --- /dev/null +++ b/psydac/api/fem_sum_form.py @@ -0,0 +1,123 @@ +import numpy as np + +from sympde.expr.expr import ( + BilinearForm as sym_BilinearForm, + LinearForm as sym_LinearForm, + Functional as sym_Functional +) + +from .basic import BasicDiscrete +from .fem import DiscreteFunctional +from .fem import DiscreteLinearForm +from .fem import DiscreteBilinearForm +from .fem_bilinear_form import DiscreteBilinearForm as DiscreteBilinearForm_SF +from .fem_common import reset_arrays +from .utilities import random_string + +__all__ = ('DiscreteSumForm',) + +#============================================================================== +class DiscreteSumForm(BasicDiscrete): + + def __init__(self, a, kernel_expr, *args, **kwargs): + + # Sum factorization is only implemented for bilinear forms in 3D, in + # which case we use it by default. A 2D implementation should be the + # next step, hence we allow the user to pass `sum_factorization=True` + # even if not supported yet. In the case of linear forms or functionals + # this option is irrelevant for now, so we ignore it. + # + # In every case we remove the `sum_factorization` key from the dict + # in order to avoid errors, because none of the class constructors + # accept this argument. + sum_factorization = kwargs.pop('sum_factorization', a.ldim == 3) + + # TODO Uncomment when the SesquilinearForm exist in SymPDE + #if not isinstance(a, (sym_BilinearForm, sym_SesquilinearForm, sym_LinearForm, sym_Functional)): + # raise TypeError('> Expecting a symbolic BilinearForm, SesquilinearForm, LinearForm, Functional') + if not isinstance(a, (sym_BilinearForm, sym_LinearForm, sym_Functional)): + raise TypeError('> Expecting a symbolic BilinearForm, LinearForm, Functional') + + self._expr = a + backend = kwargs.pop('backend', None) + self._backend = backend + + folder = kwargs.get('folder', None) + self._folder = self._initialize_folder(folder) + + # create a module name if not given + tag = random_string(8) + + # ... + forms = [] + free_args = [] + self._kernel_expr = kernel_expr + operator = None + for e in kernel_expr: + if isinstance(a, sym_LinearForm): + kwargs['update_ghost_regions'] = False + ah = DiscreteLinearForm(a, e, *args, backend=backend, **kwargs) + kwargs['vector'] = ah._vector + operator = ah._vector + + # TODO Uncomment when the SesquilinearForm exist in SymPDE + # elif isinstance(a, sym_SesquilinearForm): + # kwargs['update_ghost_regions'] = False + # ah = DiscreteSesquilinearForm(a, e, *args, assembly_backend=backend, **kwargs) + # kwargs['matrix'] = ah._matrix + # operator = ah._matrix + + elif isinstance(a, sym_BilinearForm): + kwargs['update_ghost_regions'] = False + if sum_factorization: + ah = DiscreteBilinearForm_SF(a, e, *args, assembly_backend=backend, **kwargs) + else: + ah = DiscreteBilinearForm(a, e, *args, assembly_backend=backend, **kwargs) + kwargs['matrix'] = ah._matrix + operator = ah._matrix + + elif isinstance(a, sym_Functional): + ah = DiscreteFunctional(a, e, *args, backend=backend, **kwargs) + + forms.append(ah) + free_args.extend(ah.free_args) + + if isinstance(a, sym_BilinearForm): + is_broken = len(args[0].domain)>1 + if self._backend is not None and is_broken: + for mat in kwargs['matrix']._blocks.values(): + mat.set_backend(backend) + elif self._backend is not None: + kwargs['matrix'].set_backend(backend) + + self._forms = forms + self._operator = operator + self._free_args = tuple(set(free_args)) + self._is_functional = isinstance(a, sym_Functional) + # ... + + @property + def forms(self): + return self._forms + + @property + def free_args(self): + return self._free_args + + @property + def is_functional(self): + return self._is_functional + + def assemble(self, *, reset=True, **kwargs): + if not self.is_functional: + if reset : + reset_arrays(*[i for M in self.forms for i in M.global_matrices]) + + for form in self.forms: + form.assemble(reset=False, **kwargs) + self._operator.exchange_assembly_data() + return self._operator + else: + M = [form.assemble(**kwargs) for form in self.forms] + M = np.sum(M) + return M diff --git a/psydac/api/tests/test_sum_factorization_assembly_3d.py b/psydac/api/tests/test_sum_factorization_assembly_3d.py new file mode 100644 index 000000000..7aab017a3 --- /dev/null +++ b/psydac/api/tests/test_sum_factorization_assembly_3d.py @@ -0,0 +1,592 @@ +import os +import pytest +import time +import numpy as np + +from sympy import sin, sqrt, pi, Abs, cos, tan + +from sympde.calculus import dot, cross, grad, curl, div, laplace +from sympde.expr import BilinearForm, integral +from sympde.topology import element_of, elements_of, Cube, Mapping, ScalarFunctionSpace, VectorFunctionSpace, Domain, Derham + +from psydac.api.discretization import discretize +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL +from psydac.linalg.block import BlockVectorSpace +from psydac.fem.basic import FemField + +try: + mesh_dir = os.environ['PSYDAC_MESH_DIR'] + +except: + base_dir = os.path.dirname(os.path.realpath(__file__)) + base_dir = os.path.join(base_dir, '..', '..', '..') + mesh_dir = os.path.join(base_dir, 'mesh') + + +# With PR #448, matrices corresponding to bilinear forms on 3D domains are being assembled using a so called sum factorization algorithm. +# Unless explicitely using the old algorithm, this happens automatically, and hence all old tests passing should indicate that the implementation of the sum factorization algorithm has been successful. +# Nonetheless, there are various difficulties in the implementation, and possibly not all of them are accounted for in the existing tests. + +# This file is designed to test such "difficult" edge cases - for mapped (Bspline & analytical) and parametric domains: + +# Such "difficult" edge cases are: + +# 1. bilinear forms on different spaces +# 2. (FemField / analytical / ...) weight functions +# 3. high derivatives (>=2) +# 4. (multiple) free FemFields +# 5. complicated expressions +# 6. uncommen (numpy) functions that need be imported correctly in the assembly file (here in the weight function) + +# These tests also return old and new discretization and assembly times. + +# Most of the time, being close to the "old matrix" (generated using the old assembly algorithm) will be the requirement to pass a test, +# as the old implementation has not caused problems in a long time and is considered to function properly. + +# Update: Instead of testing all mapping options all of the time, we now rather randomly test one of the three options! +#@pytest.mark.parametrize('mapping', ('None', 'Analytical', 'Bspline')) +def test_assembly(): # mapping): + + rng = np.random.default_rng() # (seed=42) + mapping_options = ['None', 'Analytical', 'Bspline'] + mapping = mapping_options[int(np.floor(rng.random()*3))] + + ncells = [7, 5, 6] + degree = [2, 4, 3] + periodic = [int(np.floor(rng.random()*2))==True for _ in range(3)] + print(f'Random periodicity: {periodic}') + + trial_multiplicity = [1, 3, 2] + test_multiplicity = [2, 2, 3] + + backend = PSYDAC_BACKEND_GPYCCEL + + if mapping == 'None': + + domain = Cube('C', bounds1=(0,1), bounds2=(0,1), bounds3=(0,1)) + derham = Derham(domain) + + domain_h = discretize(domain, ncells=ncells, periodic=periodic) + derham_h = discretize(derham, domain_h, degree=degree, multiplicity=trial_multiplicity) + derham_test_h = discretize(derham, domain_h, degree=degree, multiplicity=test_multiplicity) + + elif mapping == 'Bspline': + + filename = os.path.join(mesh_dir, 'identity_3d.h5') + + domain = Domain.from_file(filename=filename) + derham = Derham(domain) + + domain_h = discretize(domain, filename=filename) + derham_h = discretize(derham, domain_h, degree=domain.mapping.get_callable_mapping().space.degree, multiplicity=trial_multiplicity) + derham_test_h = discretize(derham, domain_h, degree=domain.mapping.get_callable_mapping().space.degree, multiplicity=test_multiplicity) + + elif mapping == 'Analytical': + + class HalfSquareTorusMapping3D(Mapping): + _expressions = {'x': 'x1 * cos(x2)', + 'y': 'x1 * sin(x2)', + 'z': 'x3'} + + _ldim = 3 + _pdim = 3 + + M = HalfSquareTorusMapping3D('M') + logical_domain = Cube('C', bounds1=(0.3,1), bounds2=(0,np.pi), bounds3=(0,1)) + + domain = M(logical_domain) + derham = Derham(domain) + + domain_h = discretize(domain, ncells=ncells, periodic=periodic) + derham_h = discretize(derham, domain_h, degree=degree, multiplicity=trial_multiplicity) + derham_test_h = discretize(derham, domain_h, degree=degree, multiplicity=test_multiplicity) + + x, y, z = domain.coordinates + weight = 1 + sqrt(Abs(x*y**2 + z)) + abs(sin(x-2*y + pi + np.pi)) + + V0 = derham.V0 + V0h = derham_h.V0 + V1 = derham.V1 + V1h = derham_h.V1 + V2 = derham.V2 + V2h = derham_h.V2 + V3 = derham.V3 + V3h = derham_h.V3 + + V0h_test = derham_test_h.V0 + V1h_test = derham_test_h.V1 + V2h_test = derham_test_h.V2 + V3h_test = derham_test_h.V3 + + Vs = ScalarFunctionSpace('Vsh', domain) + Vsh = discretize(Vs, domain_h, degree=degree) + Vvc = VectorFunctionSpace('Vvc', domain, kind='hcurl') + Vvch = discretize(Vvc, domain_h, degree=degree) + Vvd = VectorFunctionSpace('Vvd', domain, kind='hdiv') + Vvdh = discretize(Vvd, domain_h, degree=degree) + + Vsh_test = discretize(Vs, domain_h, degree=degree, multiplicity=test_multiplicity) + Vvch_test = discretize(Vs, domain_h, degree=degree, multiplicity=test_multiplicity) + Vvdh_test = discretize(Vs, domain_h, degree=degree, multiplicity=test_multiplicity) + + u1, u2, F01, F02, F03 = elements_of(V0, names='u1, u2, F01, F02, F03') + v1, v2, F11, F12, F13 = elements_of(V1, names='v1, v2, F11, F12, F13') + w1, w2, F21, F22, F23 = elements_of(V2, names='w1, w2, F21, F22, F23') + f1, f2, F31, F32, F33 = elements_of(V3, names='f1, f2, F31, F32, F33') + + fs1, fs2, Fs1, Fs2, Fs3 = elements_of(Vs, names='fs1, fs2, Fs1, Fs2, Fs3') + fvc1, fvc2, Fvc1, Fvc2, Fvc3 = elements_of(Vvc, names='fvc1, fvc2, Fvc1, Fvc2, Fvc3') + fvd1, fvd2, Fvd1, Fvd2, Fvd3 = elements_of(Vvd, names='fvd1, fvd2, Fvd1, Fvd2, Fvd3') + + trial_spaces = {'V0': {'Vh':V0h, 'funcs':[u1, u2]}, + 'V1': {'Vh':V1h, 'funcs':[v1, v2]}, + 'V2': {'Vh':V2h, 'funcs':[w1, w2]}, + 'V3': {'Vh':V3h, 'funcs':[f1, f2]}, + 'Vs': {'Vh':Vsh, 'funcs':[fs1, fs2]}, + 'Vvc':{'Vh':Vvch, 'funcs':[fvc1, fvc2]}, + 'Vvd':{'Vh':Vvdh, 'funcs':[fvd1, fvd2]}} + + test_spaces = {'V0': {'Vh':V0h_test, 'funcs':[u1, u2]}, + 'V1': {'Vh':V1h_test, 'funcs':[v1, v2]}, + 'V2': {'Vh':V2h_test, 'funcs':[w1, w2]}, + 'V3': {'Vh':V3h_test, 'funcs':[f1, f2]}, + 'Vs': {'Vh':Vsh_test, 'funcs':[fs1, fs2]}, + 'Vvc':{'Vh':Vvch_test, 'funcs':[fvc1, fvc2]}, + 'Vvd':{'Vh':Vvdh_test, 'funcs':[fvd1, fvd2]}} + + F01_coeffs = V0h.coeff_space.zeros() + rng.random(size=F01_coeffs._data.shape, dtype='float64', out=F01_coeffs._data) + F11_coeffs = V1h.coeff_space.zeros() + for block in F11_coeffs.blocks: + rng.random(size=block._data.shape, dtype='float64', out=block._data) + F12_coeffs = V1h.coeff_space.zeros() + for block in F12_coeffs.blocks: + rng.random(size=block._data.shape, dtype='float64', out=block._data) + F21_coeffs = V2h.coeff_space.zeros() + for block in F21_coeffs.blocks: + rng.random(size=block._data.shape, dtype='float64', out=block._data) + Fvc1_coeffs = Vvch.coeff_space.zeros() + for block in Fvc1_coeffs.blocks: + rng.random(size=block._data.shape, dtype='float64', out=block._data) + Fs1_coeffs = Vsh.coeff_space.zeros() + rng.random(size=Fs1_coeffs._data.shape, dtype='float64', out=Fs1_coeffs._data) + Fs2_coeffs = Vsh.coeff_space.zeros() + rng.random(size=Fs2_coeffs._data.shape, dtype='float64', out=Fs2_coeffs._data) + Fvd1_coeffs = Vvdh.coeff_space.zeros() + for block in Fvd1_coeffs.blocks: + rng.random(size=block._data.shape, dtype='float64', out=block._data) + + F01_field = FemField(V0h, F01_coeffs) + F11_field = FemField(V1h, F11_coeffs) + F12_field = FemField(V1h, F12_coeffs) + F21_field = FemField(V2h, F21_coeffs) + Fvc1_field = FemField(Vvch, Fvc1_coeffs) + Fs1_field = FemField(Vsh, Fs1_coeffs) + Fs2_field = FemField(Vsh, Fs2_coeffs) + Fvd1_field = FemField(Vvdh, Fvd1_coeffs) + + bilinear_forms = { # one and two free FemFields without derivatives (with derivatives in seperate test) + # complicated expressions + 'Q' :{'trial' :'V1', 'test':'V1', + 'expr' :dot(cross(F11, v1), cross(F11, v2)), + 'fields':[F11_field, ]}, + 'equilibrium' :{'trial' :'V1', 'test':'V1', + 'expr' :dot(cross(F11, v1), cross(v2, F12)), + 'fields':[F11_field, F12_field]}, + 'Elena' :{'trial' :'V1', 'test':'V1', + 'expr' :dot(F01*v1, v2), + 'fields':[F01_field, ]}, + # weight function, free FemField, different spaces + 'dot(grad(u),v)':{'trial' :'V0', 'test':'V1', + 'expr' :dot(grad(u1), v2)*F01*weight, + 'fields':[F01_field, ]}, + 'dot(curl(v),w)':{'trial' :'V1', 'test':'V2', + 'expr' :dot(curl(v1), F21)*div(w2)*weight, + 'fields':[F21_field, ]}, + # among other difficulties: multiple scalar FemFields + 'ScalarFields' :{'trial' :'V0', 'test':'V1', + 'expr' :dot(grad(u1), curl(Fvc1))*dot(grad(Fs1), curl(v2))*Fs2*div(Fvd1), + 'fields':[Fvc1_field, Fs1_field, Fs2_field, Fvd1_field]}, + # high derivatives, not FEEC + 'bilaplace' :{'trial' :'Vs', 'test':'Vs', + 'expr' :laplace(fs1)*laplace(fs2)} + } + + # test all BFs + # bilinear_form_strings_to_test = list(bilinear_forms.keys()) + + # or only a subset + bilinear_form_strings_to_test = list(bilinear_forms.keys())[:-1] # exclude expensive bilaplace test + + bilinear_forms_to_test = {} + for name in bilinear_form_strings_to_test: + bilinear_forms_to_test[name] = bilinear_forms[name] + + int_0 = lambda expr: integral(domain, expr) + print() + + for bf_name, bf_data in bilinear_forms_to_test.items(): + + trial_space = trial_spaces[bf_data['trial']] + Vh = trial_space['Vh'] + u = trial_space['funcs'][0] + test_space = test_spaces[bf_data['test']] + Wh = test_space ['Vh'] + v = test_space['funcs'][1] + expr = bf_data['expr'] + if 'fields' in bf_data.keys(): + fields = bf_data['fields'] + + a = BilinearForm((u, v), int_0(expr)) + + t0 = time.time() + ah_old = discretize(a, domain_h, (Vh, Wh), backend=backend, sum_factorization=False) + t1 = time.time() + discretization_time_old = t1 - t0 + + t0 = time.time() + ah = discretize(a, domain_h, (Vh, Wh), backend=backend) + t1 = time.time() + discretization_time = t1 - t0 + + if bf_name == 'Q': + t0_old = time.time() + A_old = ah_old.assemble(F11=fields[0]) + t1_old = time.time() + + t0 = time.time() + A = ah.assemble(F11=fields[0]) + t1 = time.time() + elif bf_name == 'equilibrium': + t0_old = time.time() + A_old = ah_old.assemble(F11=fields[0], F12=fields[1]) + t1_old = time.time() + + t0 = time.time() + A = ah.assemble(F11=fields[0], F12=fields[1]) + t1 = time.time() + elif bf_name in ('Elena', 'dot(grad(u),v)'): + t0_old = time.time() + A_old = ah_old.assemble(F01=fields[0]) + t1_old = time.time() + + t0 = time.time() + A = ah.assemble(F01=fields[0]) + t1 = time.time() + elif bf_name == 'dot(curl(v),w)': + t0_old = time.time() + A_old = ah_old.assemble(F21=fields[0]) + t1_old = time.time() + + t0 = time.time() + A = ah.assemble(F21=fields[0]) + t1 = time.time() + elif bf_name == 'ScalarFields': + t0_old = time.time() + A_old = ah_old.assemble(Fvc1=fields[0], Fs1=fields[1], Fs2=fields[2], Fvd1=fields[3]) + t1_old = time.time() + + t0 = time.time() + A = ah.assemble(Fvc1=fields[0], Fs1=fields[1], Fs2=fields[2], Fvd1=fields[3]) + t1 = time.time() + else: + t0_old = time.time() + A_old = ah_old.assemble() + t1_old = time.time() + + t0 = time.time() + A = ah.assemble() + t1 = time.time() + + assembly_time_old = t1_old - t0_old + assembly_time = t1 - t0 + + # Testing whether two linear operators are identical by comparing their arrays is quite expensive. + # Thus we instead test whether three random domain vectors applied to both + # the old and the new matrix produce the same codomain vector. + + domain_vector1 = Vh.coeff_space.zeros() + domain_vector2 = Vh.coeff_space.zeros() + domain_vector3 = Vh.coeff_space.zeros() + domain_vectors = [domain_vector1, domain_vector2, domain_vector3] + + if isinstance(Vh.coeff_space, BlockVectorSpace): + for domain_vector in domain_vectors: + for block in domain_vector.blocks: + rng.random(size=block._data.shape, dtype='float64', out=block._data) + else: + for domain_vector in domain_vectors: + rng.random(size=domain_vector._data.shape, dtype='float64', out=domain_vector._data) + + err = [] + rel_err = [] + + for domain_vector in domain_vectors: + codomain_vector = A @ domain_vector + codomain_vector_old = A_old @ domain_vector + + norm_old = np.sqrt(codomain_vector_old.inner(codomain_vector_old)) + + diff = codomain_vector - codomain_vector_old + + err.append(np.sqrt(diff.inner(diff))) + rel_err.append(err[-1] / norm_old) + + print(f' >>> Mapping: {mapping}') + print(f' >>> BilinearForm: {bf_name}') + print(f' >>> Discretization in: Old {discretization_time_old:.3g} \t\t || New {discretization_time:.3g} \t\t || Old/New {discretization_time_old/discretization_time:.3g}') + print(f' >>> Assembly in: Old {assembly_time_old:.3g} \t \t || New {assembly_time:.3g} \t\t || Old/New {assembly_time_old/assembly_time:.3g}') + print(f' >>> Error: {max(err):.3g}') + print(f' >>> Rel. Error: {max(rel_err):.3g}') + print() + + assert max(rel_err) < 1e-12 # arbitrary rel. error bound (How to test better?) + + +# fixed by PR #507 +#@pytest.mark.xfail +def test_allocate_matrix_bug(): + """ + This test is related to Issue #504. + + The bilinear form + (V0 x V3) ni (u, f) mapsto int_{Omega} u*f + should be the transpose of the bilinear form + (V3 x V0) ni (f, u) mapsto int_{Omega} u*f + but is not. + """ + + ncells = [15, 16, 17] + degree = [4, 3, 2] + periodic = [False, True, False] + + backend = PSYDAC_BACKEND_GPYCCEL + + domain = Cube('C', bounds1=(0,1), bounds2=(0,1), bounds3=(0,1)) + derham = Derham(domain) + + domain_h = discretize(domain, ncells=ncells, periodic=periodic) + derham_h = discretize(derham, domain_h, degree=degree) + + P0, _, _, P3 = derham_h.projectors() + + V0 = derham.V0 + V0h = derham_h.V0 + V3 = derham.V3 + V3h = derham_h.V3 + + u = element_of(V0, name='u') + f = element_of(V3, name='f') + + fun = lambda x, y, z : 1 + u_coeffs = P0(fun).coeffs + f_coeffs = P3(fun).coeffs + + a0 = BilinearForm((u, f), integral(domain, u*f)) + a1 = BilinearForm((f, u), integral(domain, u*f)) + + a0h = discretize(a0, domain_h, (V0h, V3h), backend=backend, sum_factorization=False) + a1h = discretize(a1, domain_h, (V3h, V0h), backend=backend, sum_factorization=False) + + A0 = a0h.assemble() + A1 = a1h.assemble() + A1T = A1.T + + # Clearly, it should hold A1T = A0, and further ||A0|| = ||A1||. + A0arr = A0.toarray() + A1arr = A1.toarray() + A1Tarr = A1T.toarray() + + diff1 = np.linalg.norm(A0arr - A1Tarr) + diff2 = np.linalg.norm(A0arr) - np.linalg.norm(A1arr) + + print(f' || A0 - A1.T || = {diff1:.3g}') + print(f' ||A0|| - ||A1|| = {diff2:.3g}') + + # Further, the following integral should evaluate to 1. + # This however is only the case for the second integral, + # independent on whether one uses the new or old assembly algorithm. + + print(f' 1 =? {A0.dot_inner(u_coeffs, f_coeffs)}') + print(f' 1 =? {A1.dot_inner(f_coeffs, u_coeffs)}') + + assert diff1 <= 1e-12 # arbitrary error bound + assert diff2 <= 1e-12 # arbitrary error bound + +#@pytest.mark.xfail +def test_free_FemField_derivatives(): + """ + These particular bilinear forms, when using a constant 1-vector coefficient vector for the free FemFields, + causes problems in a different test file of mine. + In particular, the assembled matrices used to have really small norms (~e-13). + That is probably due to the constant 1-vector corresponding to a constant function, + which means that all appearing derivatives of free FemFields are 0. + + When using meaningful free FemFields, these dubious observations disappeared. + + """ + + ncells = [5, 2, 4] + degree = [2, 1, 3] + periodic = [False, False, False] + + backend = PSYDAC_BACKEND_GPYCCEL + + domain = Cube('C', bounds1=(0,1), bounds2=(0,1), bounds3=(0,1)) + derham = Derham(domain) + + domain_h = discretize(domain, ncells=ncells, periodic=periodic) + derham_h = discretize(derham, domain_h, degree=degree) + + P0, P1, P2, P3 = derham_h.projectors() + + V0 = derham.V0 + V0h = derham_h.V0 + V1 = derham.V1 + V1h = derham_h.V1 + V2 = derham.V2 + V2h = derham_h.V2 + + u = element_of (V0, name= 'u') + v, F1 = elements_of(V1, names='v, F1') + w1, w2, F2 = elements_of(V2, names='w1, w2, F2') + + u_func = lambda x, y, z: x + y + z + u_coeffs = P0(u_func).coeffs + + v_1 = lambda x, y, z: 1 + v_2 = lambda x, y, z: 1 + v_3 = lambda x, y, z: 1 + v_func = (v_1, v_2, v_3) + v_coeffs = P1(v_func).coeffs + + w_1 = lambda x, y, z: x + w_2 = lambda x, y, z: y + w_3 = lambda x, y, z: z + w_func = (w_1, w_2, w_3) + w_coeffs = P2(w_func).coeffs + + dubious_observations = False + + if dubious_observations: + F1_coeffs = V1h.coeff_space.zeros() + F2_coeffs = V2h.coeff_space.zeros() + for block in F1_coeffs.blocks: + block._data = np.ones(block._data.shape, dtype='float64') + for block in F2_coeffs.blocks: + block._data = np.ones(block._data.shape, dtype='float64') + F1_FF = FemField(V1h, F1_coeffs) + F2_FF = FemField(V2h, F2_coeffs) + else: + F1_1 = lambda x, y, z: z + F1_2 = lambda x, y, z: x + F1_3 = lambda x, y, z: y + F1_func = (F1_1, F1_2, F1_3) + F1_FF = P1(F1_func) + + F2_FF = FemField(V2h, w_coeffs.copy()) + + # with the above choices (dubious_observation = False): + # a0 reduces to 3*int_{Omega}x+y+z with Omega being the unit square. Expected value: 4.5 + a0 = BilinearForm((u, v), integral(domain, dot(grad(u), curl(F1)) * dot(v, F1))) + # a1 reduces to 9* ----------------------------------------- " -------------------- 13.5 + a1 = BilinearForm((u, w2), integral(domain, dot(grad(u), F2)*div(w2)*div(F2))) + # a2 reduces to 3* ----------------------------------------- " --------------------- 4.5 + a2 = BilinearForm((w1, w2), integral(domain, dot(curl(F1), w1)*div(w2))) + + a0h_old = discretize(a0, domain_h, (V0h, V1h), backend=backend, sum_factorization=False) + a1h_old = discretize(a1, domain_h, (V0h, V2h), backend=backend, sum_factorization=False) + a2h_old = discretize(a2, domain_h, (V2h, V2h), backend=backend, sum_factorization=False) + + a0h = discretize(a0, domain_h, (V0h, V1h), backend=backend) + a1h = discretize(a1, domain_h, (V0h, V2h), backend=backend) + a2h = discretize(a2, domain_h, (V2h, V2h), backend=backend) + + bfs = [(a0h_old, a0h), (a1h_old, a1h), (a2h_old, a2h)] + print() + + for i, (ah_old, ah) in enumerate(bfs): + + if i in (0, 2): + A_old = ah_old.assemble(F1=F1_FF) + A = ah.assemble(F1=F1_FF) + else: + A_old = ah_old.assemble(F2=F2_FF) + A = ah.assemble(F2=F2_FF) + + if i == 0: + value_old = A_old.dot_inner(u_coeffs, v_coeffs) + value = A.dot_inner(u_coeffs, v_coeffs) + elif i == 1: + value_old = A_old.dot_inner(u_coeffs, w_coeffs) + value = A.dot_inner(u_coeffs, w_coeffs) + else: + value_old = A_old.dot_inner(w_coeffs, w_coeffs) + value = A.dot_inner(w_coeffs, w_coeffs) + + A_old_arr = A_old.toarray() + A_arr = A.toarray() + A_old_norm = np.linalg.norm(A_old_arr) + A_norm = np.linalg.norm(A_arr) + + err = np.linalg.norm(A_old_arr - A_arr) + rel_err = err / A_old_norm + + print(f' i = {i}') + print(f' >>> Error: {err:.3g}') + print(f' >>> Rel. Error: {rel_err:.3g}') + print(f' >>> Norms: ||A_old|| = {A_old_norm:.3g} \t\t ||A|| = {A_norm:.3g}') + print() + + # arbitrary tolerance + tol = 1e-12 + if not dubious_observations: + assert abs(value-value_old) < tol + assert rel_err < tol + +def test_assembly_free_FemFields(): + + backend = PSYDAC_BACKEND_GPYCCEL + + domain = Cube(bounds1=(2626,3179), bounds2=(-138, 138), bounds3=(-760.3, 69)) + derham = Derham(domain) + V0 = derham.V0 + V1 = derham.V1 + + u = element_of(V1, name='u') + v = element_of(V1, name='v') + p = element_of(V0, name='p') + + p_call = lambda xi,yi,zi: np.cos(2*np.pi*(xi-2626)/553) + np.tan(2*np.pi*(zi+760.3)/(5*829.3)) + x,y,z = domain.coordinates + p_sym = cos(2*np.pi*(x-2626)/553) + tan(2*np.pi*(z+760.3)/(5*829.3)) + + a_sym = BilinearForm((u, v), integral(domain, p_sym*dot(u, v))) + a_fem = BilinearForm((u, v), integral(domain, p*dot(u, v))) + + ncells = (11, 1, 17) + degree = (3, 1, 4) + domain_h = discretize(domain, ncells=ncells, periodic=(False, True, False)) + derham_h = discretize(derham, domain_h, degree=degree) + V1_h = derham_h.V1 + P0 = derham_h.projectors()[0] + + p_fem = P0(p_call) + + a_sym_h = discretize(a_sym, domain_h, (V1_h, V1_h), backend=backend) + a_fem_h = discretize(a_fem, domain_h, (V1_h, V1_h), backend=backend) + + A_sym = a_sym_h.assemble() + A_fem = a_fem_h.assemble(p=p_fem) + + A_sym_sp = A_sym.tosparse() + A_fem_sp = A_fem.tosparse() + diff = A_sym_sp - A_fem_sp + norm_sym = np.linalg.norm(A_sym_sp.data) + norm_fem = np.linalg.norm(A_fem_sp.data) + error = np.linalg.norm(diff.data) + rel_err = error / norm_sym + + #print(norm_sym, norm_fem, error, rel_err) + + assert rel_err < 1e-4 diff --git a/psydac/feec/derivatives.py b/psydac/feec/derivatives.py index f7f28ccbd..1fd3de5d2 100644 --- a/psydac/feec/derivatives.py +++ b/psydac/feec/derivatives.py @@ -9,38 +9,28 @@ from psydac.fem.vector import VectorFemSpace from psydac.fem.tensor import TensorFemSpace from psydac.linalg.basic import IdentityOperator -from psydac.fem.basic import FemField, FemSpace +from psydac.fem.basic import FemField, FemSpace, FemLinearOperator from psydac.linalg.basic import LinearOperator from psydac.ddm.cart import DomainDecomposition, CartDecomposition __all__ = ( 'DirectionalDerivativeOperator', - 'DiffOperator', - 'Derivative_1D', - 'Gradient_2D', - 'Gradient_3D', - 'ScalarCurl_2D', - 'VectorCurl_2D', - 'Curl_3D', - 'Divergence_2D', - 'Divergence_3D', - 'block_tostencil' + 'Derivative1D', + 'Gradient2D', + 'Gradient3D', + 'ScalarCurl2D', + 'VectorCurl2D', + 'Curl3D', + 'Divergence2D', + 'Divergence3D', + 'BrokenGradient2D', + 'BrokenTransposedGradient2D', + 'BrokenScalarCurl2D', + 'BrokenTransposedScalarCurl2D', ) #==================================================================================================== -def block_tostencil(M): - """ - Convert a BlockLinearOperator that contains KroneckerStencilMatrix objects - to a BlockLinearOperator that contains StencilMatrix objects - """ - blocks = [list(b) for b in M.blocks] - for i1,b in enumerate(blocks): - for i2, mat in enumerate(b): - if mat is None: - continue - blocks[i1][i2] = mat.tostencil() - return BlockLinearOperator(M.domain, M.codomain, blocks=blocks) - +# Singlepatch derivative operators #==================================================================================================== class DirectionalDerivativeOperator(LinearOperator): """ @@ -361,40 +351,7 @@ def copy(self): self._diffdir, negative=self._negative, transposed=self._transposed) #==================================================================================================== -class DiffOperator: - def __init__(self, domain, codomain, matrix): - assert isinstance(domain, FemSpace) - assert isinstance(codomain, FemSpace) - assert isinstance(matrix, LinearOperator) - assert domain.coeff_space is matrix.domain - assert codomain.coeff_space is matrix.codomain - - self._domain = domain - self._codomain = codomain - self._matrix = matrix - - @property - def matrix(self): - return self._matrix - - @property - def domain(self): - return self._domain - - @property - def codomain(self): - return self._codomain - - def __call__(self, u): - assert isinstance(u, FemField) - assert u.space == self.domain - - coeffs = self.matrix.dot(u.coeffs) - - return FemField(self.codomain, coeffs=coeffs) - -#==================================================================================================== -class Derivative_1D(DiffOperator): +class Derivative1D(FemLinearOperator): """ 1D derivative. @@ -415,10 +372,10 @@ def __init__(self, H1, L2): assert H1.degree[0] == L2.degree[0] + 1 # Store data in object - super().__init__(H1, L2, DirectionalDerivativeOperator(H1.coeff_space, L2.coeff_space, 0)) + super().__init__(fem_domain = H1, fem_codomain = L2, linop = DirectionalDerivativeOperator(H1.coeff_space, L2.coeff_space, 0)) #==================================================================================================== -class Gradient_2D(DiffOperator): +class Gradient2D(FemLinearOperator): """ Gradient operator in 2D. @@ -434,7 +391,7 @@ class Gradient_2D(DiffOperator): def __init__(self, H1, Hcurl): assert isinstance( H1, TensorFemSpace); assert H1.ldim == 2 - assert isinstance(Hcurl, VectorFemSpace); assert Hcurl.ldim == 2 + assert isinstance(Hcurl, VectorFemSpace); assert Hcurl.ldim == 2 assert Hcurl.spaces[0].periodic == H1.periodic assert Hcurl.spaces[1].periodic == H1.periodic @@ -454,11 +411,10 @@ def __init__(self, H1, Hcurl): matrix = BlockLinearOperator(H1.coeff_space, Hcurl.coeff_space, blocks=blocks) # Store data in object - super().__init__(H1, Hcurl, matrix) - + super().__init__(fem_domain = H1, fem_codomain = Hcurl, linop = matrix) #==================================================================================================== -class Gradient_3D(DiffOperator): +class Gradient3D(FemLinearOperator): """ Gradient operator in 3D. @@ -497,10 +453,10 @@ def __init__(self, H1, Hcurl): matrix = BlockLinearOperator(H1.coeff_space, Hcurl.coeff_space, blocks=blocks) # Store data in object - super().__init__(H1, Hcurl, matrix) + super().__init__(fem_domain = H1, fem_codomain = Hcurl, linop = matrix) #==================================================================================================== -class ScalarCurl_2D(DiffOperator): +class ScalarCurl2D(FemLinearOperator): """ Scalar curl operator in 2D: computes a scalar field from a vector field. @@ -536,10 +492,10 @@ def __init__(self, Hcurl, L2): matrix = BlockLinearOperator(Hcurl.coeff_space, L2.coeff_space, blocks=blocks) # Store data in object - super().__init__(Hcurl, L2, matrix) + super().__init__(fem_domain = Hcurl, fem_codomain = L2, linop = matrix) #==================================================================================================== -class VectorCurl_2D(DiffOperator): +class VectorCurl2D(FemLinearOperator): """ Vector curl operator in 2D: computes a vector field from a scalar field. This is sometimes called the 'rot' operator. @@ -576,10 +532,10 @@ def __init__(self, H1, Hdiv): matrix = BlockLinearOperator(H1.coeff_space, Hdiv.coeff_space, blocks=blocks) # Store data in object - super().__init__(H1, Hdiv, matrix) + super().__init__(fem_domain = H1, fem_codomain = Hdiv, linop = matrix) #==================================================================================================== -class Curl_3D(DiffOperator): +class Curl3D(FemLinearOperator): """ Curl operator in 3D. @@ -626,10 +582,10 @@ def __init__(self, Hcurl, Hdiv): # ... # Store data in object - super().__init__(Hcurl, Hdiv, matrix) + super().__init__(fem_domain = Hcurl, fem_codomain = Hdiv, linop = matrix) #==================================================================================================== -class Divergence_2D(DiffOperator): +class Divergence2D(FemLinearOperator): """ Divergence operator in 2D. @@ -665,10 +621,10 @@ def __init__(self, Hdiv, L2): matrix = BlockLinearOperator(Hdiv.coeff_space, L2.coeff_space, blocks=blocks) # Store data in object - super().__init__(Hdiv, L2, matrix) + super().__init__(fem_domain = Hdiv, fem_codomain = L2, linop = matrix) #==================================================================================================== -class Divergence_3D(DiffOperator): +class Divergence3D(FemLinearOperator): """ Divergence operator in 3D. @@ -707,4 +663,118 @@ def __init__(self, Hdiv, L2): matrix = BlockLinearOperator(Hdiv.coeff_space, L2.coeff_space, blocks=blocks) # Store data in object - super().__init__(Hdiv, L2, matrix) + super().__init__(fem_domain = Hdiv, fem_codomain = L2, linop = matrix) + +#==================================================================================================== +# 2D Multipatch derivative operators +#==================================================================================================== +class BrokenGradient2D(FemLinearOperator): + """ + Gradient operator in a 2D multipatch domain, + acting independently on each patch. + In general, the resulting field is therefore discontinuous, or "broken". + + Parameters + ---------- + V0h : MultipatchFemSpace + Domain of the gradient operator. + + V1h : MultipatchFemSpace + Codomain of the gradient operator. + """ + def __init__(self, V0h, V1h): + + FemLinearOperator.__init__(self, fem_domain=V0h, fem_codomain=V1h) + + D0s = [Gradient2D(V0, V1) for V0, V1 in zip(V0h.spaces, V1h.spaces)] + + self._linop = BlockLinearOperator(self.linop_domain, self.linop_codomain, blocks={ + (i, i): D0i.linop for i, D0i in enumerate(D0s)}) + + def transpose(self, conjugate=False): + # todo (MCP): define as the dual differential operator + return BrokenTransposedGradient2D(self.fem_domain, self.fem_codomain) + +# ============================================================================== +class BrokenTransposedGradient2D(FemLinearOperator): + """ + Transposed gradient operator in a 2D multipatch domain, + acting independently on each patch. + In general, the resulting field is therefore discontinuous, or "broken". + + Parameters + ---------- + V0h : MultipatchFemSpace + Codomain of the transposed gradient operator. + + V1h : MultipatchFemSpace + Domain of the transposed gradient operator. + """ + def __init__(self, V0h, V1h): + + FemLinearOperator.__init__(self, fem_domain=V1h, fem_codomain=V0h) + + D0s = [Gradient2D(V0, V1) for V0, V1 in zip(V0h.spaces, V1h.spaces)] + + self._linop = BlockLinearOperator(self.linop_domain, self.linop_codomain, blocks={ + (i, i): D0i.linop.T for i, D0i in enumerate(D0s)}) + + def transpose(self, conjugate=False): + # todo (MCP): discard + return BrokenGradient2D(self.fem_codomain, self.fem_domain) + +# ============================================================================== +class BrokenScalarCurl2D(FemLinearOperator): + """ + Scalar curl operator in a 2D multipatch domain, + acting independently on each patch. + In general, the resulting field is therefore discontinuous, or "broken". + + Parameters + ---------- + V1h : MultipatchFemSpace + Domain of the scalar curl operator. + + V2h : MultipatchFemSpace + Codomain of the scalar curl operator. + """ + def __init__(self, V1h, V2h): + + FemLinearOperator.__init__(self, fem_domain=V1h, fem_codomain=V2h) + + D1s = [ScalarCurl2D(V1, V2) for V1, V2 in zip(V1h.spaces, V2h.spaces)] + + self._linop = BlockLinearOperator(self.linop_domain, self.linop_codomain, blocks={ + (i, i): D1i.linop for i, D1i in enumerate(D1s)}) + + def transpose(self, conjugate=False): + return BrokenTransposedScalarCurl2D( + V1h=self.fem_domain, V2h=self.fem_codomain) + + +# ============================================================================== +class BrokenTransposedScalarCurl2D(FemLinearOperator): + """ + Transposed scalar curl operator in a 2D multipatch domain, + acting independently on each patch. + In general, the resulting field is therefore discontinuous, or "broken". + + Parameters + ---------- + V1h : MultipatchFemSpace + Codomain of the transposed scalar curl operator. + + V2h : MultipatchFemSpace + Domain of the transposed scalar curl operator. + """ + def __init__(self, V1h, V2h): + + FemLinearOperator.__init__(self, fem_domain=V2h, fem_codomain=V1h) + + D1s = [ScalarCurl2D(V1, V2) for V1, V2 in zip(V1h.spaces, V2h.spaces)] + + self._linop = BlockLinearOperator(self.linop_domain, self.linop_codomain, blocks={ + (i, i): D1i.linop.T for i, D1i in enumerate(D1s)}) + + def transpose(self, conjugate=False): + return BrokenScalarCurl2D(V1h=self.fem_codomain, V2h=self.fem_domain) diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_geometric_projectors.py similarity index 94% rename from psydac/feec/global_projectors.py rename to psydac/feec/global_geometric_projectors.py index 9cc1451a9..21685602c 100644 --- a/psydac/feec/global_projectors.py +++ b/psydac/feec/global_geometric_projectors.py @@ -4,27 +4,28 @@ from psydac.linalg.kron import KroneckerLinearSolver, KroneckerStencilMatrix from psydac.linalg.stencil import StencilMatrix, StencilVectorSpace -from psydac.linalg.block import BlockLinearOperator +from psydac.linalg.block import BlockLinearOperator, BlockVector from psydac.core.bsplines import quadrature_grid from psydac.utilities.quadratures import gauss_legendre from psydac.fem.basic import FemField from psydac.feec import dof_kernels from psydac.fem.tensor import TensorFemSpace -from psydac.fem.vector import VectorFemSpace +from psydac.fem.vector import VectorFemSpace, MultipatchFemSpace from psydac.ddm.cart import DomainDecomposition, CartDecomposition from psydac.utilities.utils import roll_edges from abc import ABCMeta, abstractmethod -__all__ = ('GlobalProjector', 'Projector_H1', 'Projector_Hcurl', 'Projector_Hdiv', 'Projector_L2', +__all__ = ('GlobalGeometricProjector', 'GlobalGeometricProjectorH1', 'GlobalGeometricProjectorHcurl', 'GlobalGeometricProjectorHdiv', 'GlobalGeometricProjectorL2', + 'MultipatchGeometricProjector', 'evaluate_dofs_1d_0form', 'evaluate_dofs_1d_1form', 'evaluate_dofs_2d_0form', 'evaluate_dofs_2d_1form_hcurl', 'evaluate_dofs_2d_1form_hdiv', 'evaluate_dofs_2d_2form', 'evaluate_dofs_3d_0form', 'evaluate_dofs_3d_1form', 'evaluate_dofs_3d_2form', 'evaluate_dofs_3d_3form') #============================================================================== -class GlobalProjector(metaclass=ABCMeta): +class GlobalGeometricProjector(metaclass=ABCMeta): """ Projects callable functions to some scalar or vector FEM space. @@ -44,7 +45,7 @@ class GlobalProjector(metaclass=ABCMeta): space : VectorFemSpace | TensorFemSpace Some finite element space, codomain of the projection operator. The exact structure where to use histopolation and where interpolation - has to be given by a subclass of the GlobalProjector class. + has to be given by a subclass of the GlobalGeometricProjector class. As of now, it is implicitly assumed for a VectorFemSpace, that for each direction that all spaces with interpolation are the same, and all spaces with histopolation are the same (i.e. yield the same quadrature/interpolation points etc.); so use with care on an arbitrary VectorFemSpace. @@ -407,7 +408,9 @@ def __call__(self, fun, dofs_only = False): return FemField(self._space, coeffs=coeffs) #============================================================================== -class Projector_H1(GlobalProjector): +# SINGLEPATCH PROJECTORS +#============================================================================== +class GlobalGeometricProjectorH1(GlobalGeometricProjector): """ Projector from H1 to an H1-conforming finite element space (i.e. a finite dimensional subspace of H1) constructed with tensor-product B-splines in 1, @@ -457,7 +460,7 @@ def __call__(self, fun, dofs_only = False): return super().__call__(fun, dofs_only = dofs_only) #============================================================================== -class Projector_Hcurl(GlobalProjector): +class GlobalGeometricProjectorHcurl(GlobalGeometricProjector): """ Projector from H(curl) to an H(curl)-conforming finite element space, i.e. a finite dimensional subspace of H(curl), constructed with tensor-product @@ -530,7 +533,7 @@ def __call__(self, fun, dofs_only = False): return super().__call__(fun, dofs_only = dofs_only) #============================================================================== -class Projector_Hdiv(GlobalProjector): +class GlobalGeometricProjectorHdiv(GlobalGeometricProjector): """ Projector from H(div) to an H(div)-conforming finite element space, i.e. a finite dimensional subspace of H(div), constructed with tensor-product @@ -607,7 +610,7 @@ def __call__(self, fun, dofs_only = False): return super().__call__(fun, dofs_only = dofs_only) #============================================================================== -class Projector_L2(GlobalProjector): +class GlobalGeometricProjectorL2(GlobalGeometricProjector): """ Projector from L2 to an L2-conforming finite element space (i.e. a finite dimensional subspace of L2) constructed with tensor-product M-splines in 1, @@ -666,7 +669,8 @@ def __call__(self, fun, dofs_only = False): """ return super().__call__(fun, dofs_only = dofs_only) -class Projector_H1vec(GlobalProjector): +#============================================================================== +class GlobalGeometricProjectorH1vec(GlobalGeometricProjector): """ Projector from H1^3 = H1 x H1 x H1 to a conforming finite element space, i.e. a finite dimensional subspace of H1^3, constructed with tensor-product @@ -730,6 +734,42 @@ def __call__(self, fun, dofs_only = False): """ return super().__call__(fun, dofs_only = dofs_only) +#============================================================================== +# MULTIPATCH PROJECTORS (2D) +#============================================================================== +class MultipatchGeometricProjector: + """ + Global Geometric Projector base class for multipatch domains. + + Parameters + ---------- + space : MultipatchFemSpace + Multipatch finite element space, codomain of the projection operator. + Projector : type[GlobalGeometricProjector] + Class of the projector to instantiate for each patch. + nquads : Iterable[int] + Number of quadrature points per cell along each direction. + This is a parameter passed to the constructor of Projector. + """ + + def __init__(self, space, Projector, nquads=None): + assert isinstance(space, MultipatchFemSpace) + assert isinstance(Projector, type) + assert issubclass(Projector, GlobalGeometricProjector) + + self._Vh = Vh = space + self._Ps = [Projector(V, nquads=nquads) for V in Vh.spaces] + + def __call__(self, funs): + """ + project a list of functions given in the logical domain + """ + us = [P(fun) for P, fun, in zip(self._Ps, funs)] + + u_c = BlockVector(self._Vh.coeff_space, blocks=[uj.coeffs for uj in us]) + + return FemField(self._Vh, coeffs=u_c) + #============================================================================== # 1D DEGREES OF FREEDOM #============================================================================== diff --git a/psydac/feec/hodge.py b/psydac/feec/hodge.py new file mode 100644 index 000000000..26cb702b4 --- /dev/null +++ b/psydac/feec/hodge.py @@ -0,0 +1,148 @@ +import os +import numpy as np + +from sympde.topology import elements_of +from sympde.topology.space import ScalarFunction +from sympde.calculus import dot +from sympde.expr.expr import BilinearForm +from sympde.expr.expr import integral + +from psydac.api.settings import PSYDAC_BACKENDS + +# =============================================================================== +class HodgeOperator: + """ + Change of basis operator: dual basis -> primal basis + + self._linop: matrix (LinearOperator) of the primal Hodge = this is the mass matrix ! + self.dual_linop: this is the INVERSE mass matrix (LinearOperator) + + Parameters + ---------- + Vh: + The discrete space + + domain_h: + The discrete domain of the projector + + metric : + the metric of the de Rham complex + + backend_language: + The backend used to accelerate the code + + Notes + ----- + We only support the identity metric, this implies that the dual Hodge is the inverse of the primal one. + # todo: allow for non-identity metrics + """ + + def __init__(self, Vh, domain_h, metric='identity', backend_language='python'): + + self._fem_domain = Vh + self._fem_codomain = Vh + + # FemLinearOperators + self._primal_hodge = None + self._dual_hodge = None + + # LinearOperators + self._linop = None + self._dual_linop = None + + self._domain_h = domain_h + self._backend_language = backend_language + + if not (metric == 'identity'): + raise NotImplementedError('only the identity metric is available') + + self._metric = metric + + def assemble_matrix(self): + """ + the Hodge matrix is the patch-wise multi-patch mass matrix + it is not stored by default but assembled on demand + """ + from psydac.api.discretization import discretize + from psydac.fem.basic import FemLinearOperator + + if self._linop is None: + Vh = self._fem_domain + assert Vh == self._fem_codomain + + V = Vh.symbolic_space + domain = V.domain + u, v = elements_of(V, names='u, v') + + if isinstance(u, ScalarFunction): + expr = u * v + else: + expr = dot(u, v) + + a = BilinearForm((u, v), integral(domain, expr)) + ah = discretize(a, self._domain_h, [Vh, Vh], backend=PSYDAC_BACKENDS[self._backend_language]) + + self._linop = ah.assemble() # Mass matrix in stencil format + + self._primal_hodge = FemLinearOperator(self._fem_domain, self._fem_codomain, linop=self._linop) + + def assemble_dual_matrix(self, solver ='cg', **kwargs): + """ + the dual Hodge matrix is the patch-wise inverse of the multi-patch mass matrix + it is not stored by default but computed on demand, by approximate local (patch-wise) inversion of the mass matrix + """ + from psydac.linalg.solvers import inverse + from psydac.linalg.block import BlockLinearOperator + from psydac.fem.basic import FemLinearOperator + + if self._dual_linop is None: + if not self._linop: + self.assemble_matrix() + + M = self._linop # mass matrix of the (primal) basis + + if self._fem_domain.is_multipatch: + + nrows = M.n_block_rows + ncols = M.n_block_cols + + inv_M_blocks = [list(b) for b in M.blocks] + for i in range(nrows): + Mii = M[i, i] + inv_Mii = inverse(Mii, solver=solver, **kwargs) + inv_M_blocks[i][i] = inv_Mii + + self._dual_linop = BlockLinearOperator(M.codomain, M.domain, blocks=inv_M_blocks) + self._dual_hodge = FemLinearOperator(self._fem_codomain, self._fem_domain, linop=self._dual_linop) + + else: + inv_M = inverse(M, solver=solver, **kwargs) + self._dual_hodge = FemLinearOperator(self._fem_codomain, self._fem_domain, linop=self._dual_linop) + + @property + def linop(self): + if self._linop is None: + self.assemble_matrix() + + return self._linop + + @property + def dual_linop(self): + if self._dual_linop is None: + self.assemble_dual_matrix() + + return self._dual_linop + + @property + def hodge(self): + if self._linop is None: + self.assemble_matrix() + + return self._primal_hodge + + @property + def dual_hodge(self): + if self._dual_linop is None: + self.assemble_dual_matrix() + + return self._dual_hodge diff --git a/psydac/fem/basic.py b/psydac/fem/basic.py index 09a1466fb..090243205 100644 --- a/psydac/fem/basic.py +++ b/psydac/fem/basic.py @@ -7,9 +7,9 @@ """ from abc import ABCMeta, abstractmethod -from psydac.linalg.basic import Vector +from psydac.linalg.basic import Vector, LinearOperator -__all__ = ('FemSpace', 'FemField') +__all__ = ('FemSpace', 'FemField', 'FemLinearOperator') #=============================================================================== # ABSTRACT BASE CLASS: FINITE ELEMENT SPACE @@ -380,3 +380,86 @@ def __isub__(self, other): assert self._space is other._space self._coeffs -= other._coeffs return self + +#=============================================================================== +# CONCRETE CLASS: Linear Operator acting on a FEM field +#=============================================================================== +class FemLinearOperator: + """ + Linear operators with an additional FEM layer. + There is also a shorthand access to sparse matrices as they are sometimes + used in the FEEC interfaces. + Parameters + ---------- + fem_domain : psydac.fem.basic.FemSpace + The discrete space of the domain + + fem_codomain : psydac.fem.basic.FemSpace + The discrete space of the codomain + + linop : + Linear Operator. + + """ + + def __init__(self, fem_domain, fem_codomain, *, linop=None): + assert isinstance(fem_domain, FemSpace) + assert isinstance(fem_codomain, FemSpace) + if linop is not None: + assert isinstance(linop, LinearOperator) + + self._fem_domain = fem_domain + self._fem_codomain = fem_codomain + + self._linop_domain = fem_domain.coeff_space + self._linop_codomain = fem_codomain.coeff_space + + self._linop = linop + + @property + def fem_domain(self): + return self._fem_domain + + @property + def fem_codomain(self): + return self._fem_codomain + + @property + def linop_domain(self): + return self._linop_domain + + @property + def linop_codomain(self): + return self._linop_codomain + + @property + def linop(self): + return self._linop + + def toarray(self): + return self._linop.toarray() + + def tosparse(self): + return self._linop.tosparse() + + #-------------------------------------------------------------------------- + def __call__(self, u, *, out=None): + assert isinstance(u, FemField) + assert u.space == self.fem_domain + + if self._linop is not None: + coeffs = self._linop.dot(u.coeffs) + else: + raise NotImplementedError('Class does not provide a __call__ method without a linear operator') + + return FemField(self.fem_codomain, coeffs=coeffs) + + def dot(self, f_coeffs, *, out=None): + assert isinstance(f_coeffs, Vector) + assert f_coeffs.space is self._linop_domain + + if self._linop is not None: + f = FemField(self.fem_domain, coeffs=f_coeffs) + return self(f).coeffs + else: + raise NotImplementedError('Class does not provide a dot method without a linear operator') diff --git a/psydac/fem/projectors.py b/psydac/fem/projectors.py index d40db6a8c..e13185051 100644 --- a/psydac/fem/projectors.py +++ b/psydac/fem/projectors.py @@ -1,8 +1,17 @@ import numpy as np -from psydac.linalg.kron import KroneckerDenseMatrix -from psydac.core.bsplines import hrefinement_matrix -from psydac.linalg.stencil import StencilVectorSpace +from sympde.topology import element_of +from sympde.topology.space import ScalarFunction +from sympde.topology.mapping import Mapping +from sympde.calculus import dot +from sympde.expr.expr import LinearForm, integral + +from psydac.api.settings import PSYDAC_BACKENDS + +from psydac.linalg.kron import KroneckerDenseMatrix +from psydac.core.bsplines import hrefinement_matrix +from psydac.linalg.stencil import StencilVectorSpace +from psydac.fem.basic import FemSpace __all__ = ('knots_to_insert', 'knot_insertion_projection_operator') @@ -100,3 +109,52 @@ def knot_insertion_projection_operator(domain, codomain): ops.append(np.eye(d.nbasis)) return KroneckerDenseMatrix(domain.coeff_space, codomain.coeff_space, *ops) + + +def get_dual_dofs(Vh, f, domain_h, backend_language="python", return_format='stencil_array'): + """ + return the dual dofs tilde_sigma_i(f) = < Lambda_i, f >_{L2} i = 1, .. dim(Vh)) of a given function f, as a stencil array or numpy array + + Parameters + ---------- + Vh : FemSpace + The discrete space for the dual dofs + + f : + The function used for evaluation + + domain_h : + The discrete domain corresponding to Vh + + backend_language: + The backend used to accelerate the code + + return_format: + The format of the dofs, can be 'stencil_array' or 'numpy_array' + + Returns + ------- + tilde_f: + The dual dofs + """ + + from psydac.api.discretization import discretize + + assert isinstance(Vh, FemSpace) + + V = Vh.symbolic_space + v = element_of(V, name='v') + + if Vh.is_vector_valued: + expr = dot(f,v) + else: + expr = f*v + + l = LinearForm(v, integral( V.domain, expr)) + lh = discretize(l, domain_h, Vh, backend=PSYDAC_BACKENDS[backend_language]) + tilde_f = lh.assemble() + + if return_format == 'numpy_array': + return tilde_f.toarray() + else: + return tilde_f diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py index 70f849c69..60f4e7796 100644 --- a/psydac/fem/tensor.py +++ b/psydac/fem/tensor.py @@ -1217,24 +1217,77 @@ def set_refined_space(self, ncells, new_space): self._refined_space[tuple(ncells)] = new_space # ... - def plot_2d_decomposition(self, mapping=None, refine=10): + def plot_2d_decomposition(self, mapping=None, *, refine=10, fig=None, ax=None, mpi_root=0): + """ + Plot decomposition of 2D TensorFemSpace w/ mapping to 2D physical space + + Plot the domain decomposition across MPI processes of a 2D + TensorFemSpace with a mapping between 2D logical and 2D physical spaces. + This function must be called collectively, and only the root process will make + the plot. On non-root processes the arguments `fig` and `ax` must be None. + + Parameters + ---------- + mapping : BasicCallableMapping + Mapping from (eta1, eta2) to (x1, x2). + + refine : int, default=10 + Cell refinement along the logical dimensions eta1 and eta2. + + fig : plt.Figure, optional + Figure where the plot should be made. Must be None on non-root processes. + ax : plt.Axes, optional + Axes where the plot should be made. Must be None on non-root processes. + + mpi_root: int, default=0 + The rank of the MPI root process which should create the plot. + + Returns + ------- + plt.Figure + Figure where the plot was made. Coincides with `fig` if provided. + """ import matplotlib.pyplot as plt from matplotlib.patches import Polygon, Patch + from sympde.topology.mapping import BasicCallableMapping from psydac.utilities.utils import refine_array_1d + # Sanity check + assert self.ldim == 2, "Function only works in 2D" + + # Check mapping if mapping is None: mapping = lambda eta: eta else: - assert mapping.ldim == self.ldim == 2 - assert mapping.pdim == self.ldim == 2 + assert isinstance(mapping, BasicCallableMapping) + assert mapping.ldim == 2, "Domain of argument `mapping` must be 2D" + assert mapping.pdim == 2, "Codomain of argument `mapping` must be 2D" - assert refine >= 1 - N = int(refine) - V1, V2 = self.spaces + # Check refine argument + assert isinstance(refine, int), f"Argument `refine` must be int, got {type(refine)} instead" + assert refine >= 1, f"Argument `refine` must be >= 1, got {refine} instead" + # Extract information about MPI communicator mpi_comm = self.coeff_space.cart.comm mpi_rank = mpi_comm.rank + mpi_size = mpi_comm.size + + # Check mpi_root argument + assert isinstance(mpi_root, int), f"Argument `mpi_root` must be int, got {type(mpi_root)} instead" + assert mpi_root >= 0, f"Argument `mpi_root` must be >= 0, got {mpi_root} instead" + assert mpi_root < mpi_size, f"Argument `mpi_root` must be smaller than communicator size ({mpi_size}), got {mpi_root} instead" + + # Check fig and ax arguments + if mpi_rank == mpi_root: + assert isinstance(fig, plt.Figure) or fig is None, f"Argument `fig` must be matplotlib Figure, got {type(fig)} instead" + assert isinstance(ax, plt.Axes) or ax is None, f"Argument `ax` must be matplotlib Axes, got {type(ax)} instead" + else: + assert fig is None, f"Argument `fig` must be None on non-root process with rank {mpi_rank}" + assert ax is None, f"Argument `ax` must be None on non-root process with rank {mpi_rank}" + + N = refine + V1, V2 = self.spaces # Local grid, refined [sk1, sk2], [ek1, ek2] = self.local_domain @@ -1251,23 +1304,62 @@ def plot_2d_decomposition(self, mapping=None, refine=10): poly = Polygon(xy, edgecolor='None') # Gather polygons on master process - polys = mpi_comm.gather(poly) + polys = mpi_comm.gather(poly, root=mpi_root) + + # Gather (s1, s2, e1, e2) on root + if mpi_rank == mpi_root: + s1_all = np.empty(mpi_size, dtype=int) + s2_all = np.empty(mpi_size, dtype=int) + e1_all = np.empty(mpi_size, dtype=int) + e2_all = np.empty(mpi_size, dtype=int) + else: + s1_all = None + s2_all = None + e1_all = None + e2_all = None + + mpi_comm.Gather(sk1 * N, s1_all, root=mpi_root) + mpi_comm.Gather(sk2 * N, s2_all, root=mpi_root) + mpi_comm.Gather((ek1 + 1) * N, e1_all, root=mpi_root) + mpi_comm.Gather((ek2 + 1) * N, e2_all, root=mpi_root) + + # Gather pcoords on root + # TODO: use Gatherv, and NumPy arrays as buffers + gathered_pcoords = mpi_comm.gather(pcoords, root=mpi_root) #------------------------------- # Non-master processes stop here - if mpi_rank != 0: + if mpi_rank != mpi_root: return #------------------------------- - # Global grid, refined - eta1 = refine_array_1d(V1.breaks, N) - eta2 = refine_array_1d(V2.breaks, N) - pcoords = np.array([[mapping(e1, e2) for e2 in eta2] for e1 in eta1]) - xx = pcoords[:, :, 0] - yy = pcoords[:, :, 1] + # Reconstruct global grid (refined) on root process + global_shape = ((V1.breaks.size - 1) * N + 1, + (V2.breaks.size - 1) * N + 1, + 2) + pcoords_global = np.empty(global_shape) + + for rank in range(mpi_comm.size): + s1 = s1_all[rank] + e1 = e1_all[rank] + s2 = s2_all[rank] + e2 = e2_all[rank] + pcoords_global[s1:e1+1, s2:e2+1, :] = gathered_pcoords[rank] + + xx = pcoords_global[:, :, 0] + yy = pcoords_global[:, :, 1] + + # If fig or ax are given, get one from the other. Otherwise create new ones + if fig and ax: + assert ax in fig.axes, "Argument `ax` must be in `fig.axes`" + elif fig: + ax = fig.gca() + elif ax: + fig = ax.figure + else: + fig, ax = plt.subplots(1, 1) # Plot decomposed domain - fig, ax = plt.subplots(1, 1) colors = itertools.cycle(plt.rcParams['axes.prop_cycle'].by_key()['color']) handles = [] for i, (poly, color) in enumerate(zip(polys, colors)): @@ -1286,7 +1378,8 @@ def plot_2d_decomposition(self, mapping=None, refine=10): ax.set_aspect('equal') ax.legend(handles=handles, bbox_to_anchor=(1.05, 1), loc=2) fig.tight_layout() - fig.show() + + return fig # ... def __str__(self): diff --git a/psydac/fem/tests/data/decomp_analytical_1_procs.png b/psydac/fem/tests/data/decomp_analytical_1_procs.png new file mode 100644 index 000000000..f4066700f Binary files /dev/null and b/psydac/fem/tests/data/decomp_analytical_1_procs.png differ diff --git a/psydac/fem/tests/data/decomp_analytical_4_procs.png b/psydac/fem/tests/data/decomp_analytical_4_procs.png new file mode 100644 index 000000000..9eedf7529 Binary files /dev/null and b/psydac/fem/tests/data/decomp_analytical_4_procs.png differ diff --git a/psydac/fem/tests/data/decomp_spline_1_procs.png b/psydac/fem/tests/data/decomp_spline_1_procs.png new file mode 100644 index 000000000..5e4ffe257 Binary files /dev/null and b/psydac/fem/tests/data/decomp_spline_1_procs.png differ diff --git a/psydac/fem/tests/data/decomp_spline_4_procs.png b/psydac/fem/tests/data/decomp_spline_4_procs.png new file mode 100644 index 000000000..6160ff220 Binary files /dev/null and b/psydac/fem/tests/data/decomp_spline_4_procs.png differ diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index c298d687b..7b688074c 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -435,15 +435,71 @@ def H(self): def idot(self, v, out): """ - Implements out += self @ v with a temporary. - Subclasses should provide an implementation without a temporary. + Implements `out += self @ v` without a temporary, using a work array. + + This default implementation uses a local work array to store the result + of `self @ v`, and then sums it to the vector `out`. This doubles the + amount of read/write operations from/to local memory. If possible, + subclasses should provide a more efficient implementation which does + not use work arrays. + + Parameters + ---------- + v : Vector + The vector to which the linear operator `self` is applied. It must + belong to the domain of `self`. + + out : Vector + The vector to be incremented by `self @ v`. It must belong to the + codomain of `self`. """ - assert isinstance(v, Vector) - assert v.space == self.domain + assert isinstance( v, Vector) assert isinstance(out, Vector) - assert out.space == self.codomain - out += self.dot(v) + assert v.space is self.domain + assert out.space is self.codomain + + if not hasattr(self, '_work'): + self._work = self.codomain.zeros() + + self.dot(v, out=self._work) + out += self._work + + def dot_inner(self, v, w): + """ + Compute the inner product of (self @ v) with w, without a temporary. + + This is equivalent to self.dot(v).inner(w), but avoids the creation of + a temporary vector because the result of self.dot(v) is stored in a + local work array. If self is a positive-definite operator, this + operation is a (weighted) inner product. + + Parameters + ---------- + v : Vector + The vector to which the linear operator (self) is applied. It must + belong to the domain of self. + + w : Vector + The second vector in the inner product. It must belong to the + codomain of self. + + Returns + ------- + float | complex + The result of the inner product between (self @ v) and w. If the + field of self is real, this is a real number. If the field of self + is complex, this is a complex number. + """ + assert isinstance(v, Vector) + assert isinstance(w, Vector) + assert v.space is self.domain + assert w.space is self.codomain + + if not hasattr(self, '_work'): + self._work = self.codomain.zeros() + + return self.dot(v, out=self._work).inner(w) def dot_inner(self, v, w): """ @@ -695,12 +751,15 @@ def operator(self): def dtype(self): return None + def set_scalar(self, c): + """ Modifies the scalar with which this LinearOperator is multiplied. E.g. for updating the stepsize.""" + self._scalar = c + def toarray(self): - return self._scalar*self._operator.toarray() + return self._scalar * self._operator.toarray() def tosparse(self): - from scipy.sparse import csr_matrix - return self._scalar*csr_matrix(self._operator.toarray()) + return self._scalar * self._operator.tosparse().tocsr() def transpose(self, conjugate=False): return ScaledLinearOperator(domain=self.codomain, codomain=self.domain, c=self._scalar if not conjugate else np.conjugate(self._scalar), A=self._operator.transpose(conjugate=conjugate)) @@ -758,7 +817,11 @@ def __init__(self, domain, codomain, *args): self._domain = domain self._codomain = codomain self._addends = addends + self._out = codomain.zeros() + #------------------------------------- + # Abstract interface + #------------------------------------- @property def domain(self): """ The domain of the linear operator, element of class ``VectorSpace``. """ @@ -769,26 +832,39 @@ def codomain(self): """ The codomain of the linear operator, element of class ``VectorSpace``. """ return self._codomain - @property - def addends(self): - """ A tuple containing the addends of the linear operator, elements of class ``LinearOperator``. """ - return self._addends - @property def dtype(self): return None + def tosparse(self): + from scipy.sparse import csr_matrix + out = csr_matrix(self.shape, dtype=self.dtype) + for a in self._addends: + out += a.tosparse() + return out + def toarray(self): out = np.zeros(self.shape, dtype=self.dtype) for a in self._addends: out += a.toarray() return out - def tosparse(self): - from scipy.sparse import csr_matrix - out = csr_matrix(self.shape, dtype=self.dtype) - for a in self._addends: - out += a.tosparse() + def dot(self, v, out=None): + """ Evaluates SumLinearOperator object at a vector v element of domain. """ + + assert isinstance(v, Vector) + assert v.space is self.domain + + if out is not None: + assert isinstance(out, Vector) + assert out.space is self.codomain + out *= 0 + else: + out = self.codomain.zeros() + + for A in self._addends: + A.idot(v, out) + return out def transpose(self, conjugate=False): @@ -797,6 +873,14 @@ def transpose(self, conjugate=False): t_addends = (*t_addends, a.transpose(conjugate=conjugate)) return SumLinearOperator(self.codomain, self.domain, *t_addends) + #-------------------------------------- + # Other properties/methods + #-------------------------------------- + @property + def addends(self): + """ A tuple containing the addends of the linear operator, elements of class ``LinearOperator``. """ + return self._addends + @staticmethod def simplify(addends): """ Simplifies a sum of linear operators by combining addends of the same class. """ @@ -819,23 +903,6 @@ def simplify(addends): out = (*out, A) return out - def dot(self, v, out=None): - """ Evaluates SumLinearOperator object at a vector v element of domain. """ - assert isinstance(v, Vector) - assert v.space == self.domain - if out is not None: - assert isinstance(out, Vector) - assert out.space == self.codomain - out *= 0 - for a in self._addends: - a.idot(v, out) - return out - else: - out = self.codomain.zeros() - for a in self._addends: - a.idot(v, out=out) - return out - #=============================================================================== class ComposedLinearOperator(LinearOperator): r""" @@ -1294,7 +1361,7 @@ def dot(self, v, out=None, **kwargs): self._dot(v, out=out, **kwargs) else: # provided dot product does not take an out argument: we simply copy the result into out - self._dot(v).copy(out=out, **kwargs) + self._dot(v, **kwargs).copy(out=out) return out diff --git a/psydac/linalg/sparse.py b/psydac/linalg/sparse.py new file mode 100644 index 000000000..71fa29fed --- /dev/null +++ b/psydac/linalg/sparse.py @@ -0,0 +1,114 @@ +# coding: utf-8 + +from scipy.sparse import sparray, csr_array, bsr_array +from scipy.sparse import spmatrix, csr_matrix, bsr_matrix + +from psydac.linalg.basic import LinearOperator +from psydac.linalg.basic import VectorSpace, Vector, LinearOperator +from psydac.linalg.stencil import StencilVector +from psydac.linalg.block import BlockVector + +__all__ = ( + 'SparseMatrixLinearOperator', +) + +class SparseMatrixLinearOperator(LinearOperator): + """ + LinearOperator representation of a sparse matrix. + + Parameters + ---------- + domain : VectorSpace + The domain of the operator. + + codomain : VectorSpace + The codomain of the operator. + + sparse_matrix : scipy.sparse.sparray | scipy.sparse.spmatrix + The sparse SciPy matrix representing the operator. Recommended formats are + CSR and BSR. Any other format will be converted to CSR (csr_array). + """ + + def __init__(self, domain, codomain, sparse_matrix): + + assert isinstance(domain, VectorSpace) + assert isinstance(codomain, VectorSpace) + assert isinstance(sparse_matrix, (sparray, spmatrix)) + + if not isinstance(sparse_matrix, + (csr_array, csr_matrix, + bsr_array, bsr_matrix)): + sparse_matrix = sparse_matrix.tocsr() + + if domain.parallel: + raise NotImplementedError('Parallel SparseMatrixLinearOperator not supported yet.') + + self._domain = domain + self._codomain = codomain + self._matrix = sparse_matrix + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._codomain + + @property + def dtype(self): + return self._matrix.dtype + + def toarray(self): + return self._matrix.toarray() + + def tosparse(self): + return self._matrix + + def transpose(self, conjugate=False): + if conjugate: + return SparseMatrixLinearOperator(self.codomain, self.domain, self._matrix.getH().tocsr()) + else: + return SparseMatrixLinearOperator(self.codomain, self.domain, self._matrix.T.tocsr()) + + def dot(self, v, out=None): + assert isinstance(v, Vector) + assert v.space is self.domain + + if out is not None: + assert isinstance(out, Vector) + assert out.space is self.codomain + out *= 0 + else: + out = self.codomain.zeros() + + self._dot_recursive(v, out=out) + + return out + + def _dot_recursive(self, v, out, ind_V=0, ind_W=0): + V = v.space + W = out.space + + if isinstance(v, StencilVector): + index_global_W = tuple(slice(s, e+1) for s, e in zip(W.starts, W.ends)) + index_global_V = tuple(slice(s, e+1) for s, e in zip(V.starts, V.ends)) + + dim_W = W.dimension + dim_V = V.dimension + + out[index_global_W].flat += self._matrix[ind_W:ind_W+dim_W, ind_V:ind_V+dim_V] @ v[index_global_V].flat + + elif isinstance(v, BlockVector): + + offset_i = ind_W + for (i, Wi) in enumerate(W.spaces): + + offset_j = ind_V + for (j, Vj) in enumerate(V.spaces): + + self._dot_recursive(v[j], out[i], ind_V=offset_j, ind_W=offset_i) + + offset_j += Vj.dimension + + offset_i += Wi.dimension diff --git a/psydac/linalg/tests/test_block.py b/psydac/linalg/tests/test_block.py index 00badb580..65894f7b1 100644 --- a/psydac/linalg/tests/test_block.py +++ b/psydac/linalg/tests/test_block.py @@ -11,6 +11,7 @@ from psydac.linalg.block import BlockVectorSpace, BlockVector from psydac.linalg.block import BlockLinearOperator from psydac.linalg.utilities import array_to_psydac, petsc_to_psydac +from psydac.linalg.sparse import SparseMatrixLinearOperator from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL from psydac.ddm.cart import DomainDecomposition, CartDecomposition @@ -728,6 +729,87 @@ def test_block_linear_operator_serial_dot( dtype, n1, n2, p1, p2, P1, P2 ): @pytest.mark.parametrize( 'dtype', [float] ) @pytest.mark.parametrize( 'n1', [8, 16] ) @pytest.mark.parametrize( 'n2', [8, 12] ) +@pytest.mark.parametrize( 'p1', [1, 3] ) +@pytest.mark.parametrize( 'p2', [1, 2] ) +@pytest.mark.parametrize( 'P1', [True, False] ) +@pytest.mark.parametrize( 'P2', [True] ) + +def test_sparse_matrix_linear_operator_serial_dot( dtype, n1, n2, p1, p2, P1, P2 ): + # set seed for reproducibility + seed(n1*n2*p1*p2) + + D = DomainDecomposition([n1,n2], periods=[P1,P2]) + + # Partition the points + npts = [n1,n2] + global_starts, global_ends = compute_global_starts_ends(D, npts) + + cart = CartDecomposition(D, npts, global_starts, global_ends, pads=[p1,p2], shifts=[1,1]) + + # Create vector spaces, stencil matrices, and stencil vectors + V = StencilVectorSpace( cart, dtype=dtype ) + M1 = StencilMatrix( V, V) + M2 = StencilMatrix( V, V ) + M3 = StencilMatrix( V, V ) + x1 = StencilVector( V ) + x2 = StencilVector( V ) + + # Fill in stencil matrices based on diagonal index + if dtype==complex: + f=lambda k1,k2: 10j*k1+k2 + else: + f=lambda k1,k2: 10*k1+k2 + + for k1 in range(-p1,p1+1): + for k2 in range(-p2,p2+1): + M1[:,:,k1,k2] = f(k1,k2) + M2[:,:,k1,k2] = f(k1,k2)+2. + M3[:,:,k1,k2] = f(k1,k2)+5. + + M1.remove_spurious_entries() + M2.remove_spurious_entries() + M3.remove_spurious_entries() + + # Fill in vector with random values, then update ghost regions + for i1 in range(n1): + for i2 in range(n2): + x1[i1,i2] = 2.0*random() - 1.0 + x2[i1,i2] = 5.0*random() - 1.0 + x1.update_ghost_regions() + x2.update_ghost_regions() + + W = BlockVectorSpace(V, V) + + # Construct a BlockLinearOperator object containing M1, M2, M, using 3 ways + # |M1 M2| + # L = | | + # |M3 0 | + + dict_blocks = {(0,0):M1, (0,1):M2, (1,0):M3} + + L = BlockLinearOperator( W, W, blocks=dict_blocks ) + Lm = SparseMatrixLinearOperator(W, W, L.tosparse().tocsr()) + + # Construct a BlockVector object containing x1 and x2 + # |x1| + # X = | | + # |x2| + + X = BlockVector(W) + X[0] = x1 + X[1] = x2 + + # Compute BlockLinearOperator product + Y = L.dot(X) + + Ym = Lm.dot(X) + + # Check data in 1D array + assert np.allclose( Ym.toarray(), Y.toarray(), rtol=1e-12, atol=1e-12 ) +#=============================================================================== +@pytest.mark.parametrize( 'dtype', [float, complex] ) +@pytest.mark.parametrize( 'n1', [8, 16] ) +@pytest.mark.parametrize( 'n2', [8, 12] ) @pytest.mark.parametrize( 'p1', [1, 2] ) @pytest.mark.parametrize( 'p2', [1, 3] ) @pytest.mark.parametrize( 'P1', [True, False] ) diff --git a/psydac/linalg/utilities.py b/psydac/linalg/utilities.py index 57a9a6b86..fb13a2e6d 100644 --- a/psydac/linalg/utilities.py +++ b/psydac/linalg/utilities.py @@ -4,14 +4,14 @@ from math import sqrt from psydac.linalg.basic import Vector -from psydac.linalg.stencil import StencilVectorSpace, StencilVector +from psydac.linalg.stencil import StencilVector, StencilVectorSpace from psydac.linalg.block import BlockVector, BlockVectorSpace from psydac.linalg.topetsc import petsc_local_to_psydac, get_npts_per_block __all__ = ( 'array_to_psydac', 'petsc_to_psydac', - '_sym_ortho' + '_sym_ortho', ) #============================================================================== diff --git a/pyproject.toml b/pyproject.toml index e7717ce09..080d99941 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "psydac" -version = "2.5.0.dev0" +version = "2.6.0.dev0" description = "Python package for isogeometric analysis (IGA)" readme = "README.md" requires-python = ">= 3.10" @@ -43,6 +43,7 @@ test = [ "pytest-cov >= 5.0.0", 'pytest >= 4.5', 'pytest-xdist >= 1.16', + 'Pillow', # Python Imaging Library (PIL) fork ] mpi = [ 'mpi4py >= 4',