diff --git a/.github/workflows/reusable-profiling-clusters.yml b/.github/workflows/reusable-profiling-clusters.yml index 1459915d6..c01ef9870 100644 --- a/.github/workflows/reusable-profiling-clusters.yml +++ b/.github/workflows/reusable-profiling-clusters.yml @@ -104,7 +104,7 @@ jobs: - name: Submit SLURM job, package and push profiling data run: | source .venv-fortran/bin/activate - python profiling/submit_diocotron_job.py \ + python profiling/examples/ToyGyrokinetic/submit_diocotron_strong_scaling.py \ --language ${{ matrix.compile-language }} \ --compiler GNU \ --upload diff --git a/profiling/examples/Poisson/cube_strong_scaling/params_poisson.py b/profiling/examples/Poisson/cube_strong_scaling/params_poisson.py index 379059283..40f4962fb 100644 --- a/profiling/examples/Poisson/cube_strong_scaling/params_poisson.py +++ b/profiling/examples/Poisson/cube_strong_scaling/params_poisson.py @@ -6,9 +6,10 @@ name = "Poisson strong scaling on 3D cube, no preconditioner" description = """ -Strong scaling test for Poisson equation on a 3D cube. -The manufactured solution is a Gaussian bump times sin(pi*x/Lx), exciting the full spectrum. -Homogeneous Dirichlet boundary conditions are set in direction x. +Strong scaling test for Poisson equation on a 3D cube (Cuboid). +The manufactured solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5), +exciting the full spectrum. +Homogeneous Dirichlet boundary conditions are set in direction eta1. The linear system is solved with unpreconditioned CG. """ @@ -128,30 +129,40 @@ from struphy.initial.base import GenericPerturbation -# Gaussian bump times sin(pi*x/Lx) (exact Dirichlet in x); the bump is narrow enough -# to be periodic in y and z up to exp(-(Ly/2)**2/w**2) ~ 1e-11. Being localized, it -# excites the whole spectrum of the Laplacian, so CG has to work for convergence. -w = 0.3 -x0, y0, z0 = Lx / 2, 0.0, Lz / 2 +# The exact solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5). +# Its value at the Dirichlet boundary eta1 = 0, 1 is exp(-0.25 / w**2) ~ 1e-11, and it is periodic +# in eta2 and eta3 up to the same accuracy. Being localized, it excites the whole spectrum of the +# Laplacian, so CG has to work for convergence. Both the solution and the right-hand side +# -Laplace(phi) are 0-forms given as functions of the logical coordinates (given_in_basis="0"); +# the pushed-forward solution at x = F(eta) is phi(eta). +w = 0.1 +eta0 = (0.5, 0.5, 0.5) -def exact_solution(x, y, z): - gauss = np.exp(-((x - x0) ** 2 + (y - y0) ** 2 + (z - z0) ** 2) / w**2) - return np.sin(np.pi / Lx * x) * gauss +def gaussian_derivatives(e1, e2, e3): + """The Gaussian and its first and second derivatives with respect to eta1, eta2, eta3.""" + eta = (e1, e2, e3) + g = np.exp(-sum((e - c) ** 2 for e, c in zip(eta, eta0)) / w**2) + d = [-2 * (e - c) / w**2 * g for e, c in zip(eta, eta0)] + dd = [(4 * (e - c) ** 2 / w**4 - 2 / w**2) * g for e, c in zip(eta, eta0)] + return g, d, dd -def rhs_fun(x, y, z): - # -Laplace(s*g) = -(s'' g + 2 s' dg/dx + s Laplace(g)) - r2 = (x - x0) ** 2 + (y - y0) ** 2 + (z - z0) ** 2 - gauss = np.exp(-r2 / w**2) - s = np.sin(np.pi / Lx * x) - ds = np.pi / Lx * np.cos(np.pi / Lx * x) - lap_gauss = (4 * r2 / w**4 - 6 / w**2) * gauss - dgauss_dx = -2 * (x - x0) / w**2 * gauss - return (np.pi / Lx) ** 2 * s * gauss - 2 * ds * dgauss_dx - s * lap_gauss +def exact_solution(e1, e2, e3): + return gaussian_derivatives(e1, e2, e3)[0] -rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="physical") +def laplacian(e1, e2, e3): + """Physical Laplacian of the exact solution, in logical coordinates (Cuboid: x_i = l_i + L_i * eta_i).""" + g, d, dd = gaussian_derivatives(e1, e2, e3) + return dd[0] / Lx**2 + dd[1] / Ly**2 + dd[2] / Lz**2 + + +def rhs_fun(e1, e2, e3): + return -laplacian(e1, e2, e3) + + +rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="0") model.em_fields.source.add_perturbation(rhs_perturbation) @@ -178,126 +189,64 @@ def rhs_fun(x, y, z): from matplotlib import pyplot as plt - def plot_slices(num, exact, name, slice_pt_x=0, slice_pt_y=0, slice_pt_z=0): + def slice_axes(k, X, Y, Z, E1, E2, E3): + """Axes (horizontal, vertical, labels) for plotting the slice eta_k = 0.5 in physical coordinates.""" + if k == 0: + return Y, Z, "y", "z" + elif k == 1: + return X, Z, "x", "z" + else: + return X, Y, "x", "y" + def plot_slices(num, exact, name): + """Numerical solution, exact solution and error on the slices eta_k = 0.5 (k = 1, 2, 3).""" fig = plt.figure(figsize=(16, 12)) - - plt.subplot(3, 3, 1) - plt.pcolor(y[slice_pt_x, :, :], z[slice_pt_x, :, :], num[slice_pt_x, :, :]) - plt.colorbar() - plt.xlabel("y") - plt.ylabel("z") - plt.title("{} from struphy, slice at x = {:.2f}".format(name, x[slice_pt_x, 0, 0])) - - plt.subplot(3, 3, 4) - plt.pcolor( - y[slice_pt_x, :, :], - z[slice_pt_x, :, :], - exact(x[slice_pt_x, :, :], y[slice_pt_x, :, :], z[slice_pt_x, :, :]), - ) - plt.colorbar() - plt.xlabel("y") - plt.ylabel("z") - plt.title("{} exact, slice at x = {:.2f}".format(name, x[slice_pt_x, 0, 0])) - - plt.subplot(3, 3, 7) - plt.pcolor( - y[slice_pt_x, :, :], - z[slice_pt_x, :, :], - np.abs(num[slice_pt_x, :, :] - exact(x[slice_pt_x, :, :], y[slice_pt_x, :, :], z[slice_pt_x, :, :])), - ) - plt.colorbar() - plt.xlabel("y") - plt.ylabel("z") - plt.title("{} error, slice at x = {:.2f}".format(name, x[slice_pt_x, 0, 0])) - - plt.subplot(3, 3, 2) - plt.pcolor(x[:, slice_pt_y, :], z[:, slice_pt_y, :], num[:, slice_pt_y, :]) - plt.colorbar() - plt.xlabel("x") - plt.ylabel("z") - plt.title("{} from struphy, slice at y = {:.2f}".format(name, y[0, slice_pt_y, 0])) - - plt.subplot(3, 3, 5) - plt.pcolor( - x[:, slice_pt_y, :], - z[:, slice_pt_y, :], - exact(x[:, slice_pt_y, :], y[:, slice_pt_y, :], z[:, slice_pt_y, :]), - ) - plt.colorbar() - plt.xlabel("x") - plt.ylabel("z") - plt.title("{} exact, slice at y = {:.2f}".format(name, y[0, slice_pt_y, 0])) - - plt.subplot(3, 3, 8) - plt.pcolor( - x[:, slice_pt_y, :], - z[:, slice_pt_y, :], - np.abs(num[:, slice_pt_y, :] - exact(x[:, slice_pt_y, :], y[:, slice_pt_y, :], z[:, slice_pt_y, :])), - ) - plt.colorbar() - plt.xlabel("x") - plt.ylabel("z") - plt.title("{} error, slice at y = {:.2f}".format(name, y[0, slice_pt_y, 0])) - - plt.subplot(3, 3, 3) - plt.pcolor(x[:, :, slice_pt_z], y[:, :, slice_pt_z], num[:, :, slice_pt_z]) - plt.colorbar() - plt.xlabel("x") - plt.ylabel("y") - plt.title("{} from struphy, slice at z = {:.2f}".format(name, z[0, 0, slice_pt_z])) - - plt.subplot(3, 3, 6) - plt.pcolor( - x[:, :, slice_pt_z], - y[:, :, slice_pt_z], - exact(x[:, :, slice_pt_z], y[:, :, slice_pt_z], z[:, :, slice_pt_z]), - ) - plt.colorbar() - plt.xlabel("x") - plt.ylabel("y") - plt.title("{} exact, slice at z = {:.2f}".format(name, z[0, 0, slice_pt_z])) - - plt.subplot(3, 3, 9) - plt.pcolor( - x[:, :, slice_pt_z], - y[:, :, slice_pt_z], - np.abs(num[:, :, slice_pt_z] - exact(x[:, :, slice_pt_z], y[:, :, slice_pt_z], z[:, :, slice_pt_z])), - ) - plt.colorbar() - plt.xlabel("x") - plt.ylabel("y") - plt.title("{} error, slice at z = {:.2f}".format(name, z[0, 0, slice_pt_z])) - + for k in range(3): + idx = tuple(num.shape[k] // 2 if i == k else slice(None) for i in range(3)) + h, v, hlabel, vlabel = slice_axes(k, X[idx], Y[idx], Z[idx], E1[idx], E2[idx], E3[idx]) + eta_k = (E1, E2, E3)[k][idx].flat[0] + for row, (data, what) in enumerate( + [(num[idx], "from struphy"), (exact[idx], "exact"), (np.abs(num[idx] - exact[idx]), "error")] + ): + plt.subplot(3, 3, 3 * row + k + 1) + plt.pcolormesh(h, v, data, shading="gouraud") + plt.colorbar() + plt.xlabel(hlabel) + plt.ylabel(vlabel) + plt.title(f"{name} {what}, slice at eta{k + 1} = {eta_k:.2f}") + fig.tight_layout() return fig if sim.rank == 0: # Raw FEEC fields are evaluated directly from the saved spline coefficients at the # cell centres of the simulation grid (serial Derham, safe on rank 0 only). - rhs_data = out.evaluate("em_fields/source", t=0) - rhs = rhs_data.values + rhs = out.evaluate("em_fields/source", t=0).values phi_data = out.evaluate("em_fields/phi", t=-1) phi = phi_data.values - x, y, z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) - - slice_pt_x = x.shape[0] // 2 - slice_pt_y = y.shape[1] // 2 - slice_pt_z = 0 - fig_rhs = plot_slices(rhs, rhs_fun, "RHS", slice_pt_x=slice_pt_x, slice_pt_y=slice_pt_y, slice_pt_z=slice_pt_z) - fig_phi = plot_slices( - phi, exact_solution, "Phi", slice_pt_x=slice_pt_x, slice_pt_y=slice_pt_y, slice_pt_z=slice_pt_z + # logical coordinates of the evaluation points and their images x = F(eta) + E1, E2, E3 = np.meshgrid( + *(phi_data.coords[c].values for c in ("eta1", "eta2", "eta3")), + indexing="ij", ) + X, Y, Z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) + + # pushed-forward exact solutions: the 0-forms at x = F(eta) are their values at eta + rhs_exact = rhs_fun(E1, E2, E3) + phi_exact = exact_solution(E1, E2, E3) + + fig_rhs = plot_slices(rhs, rhs_exact, "RHS") + fig_phi = plot_slices(phi, phi_exact, "Phi") - rel_err_rhs = np.max(np.abs(rhs - rhs_fun(x, y, z))) / np.max(np.abs(rhs_fun(x, y, z))) - rel_err_phi = np.max(np.abs(phi - exact_solution(x, y, z))) / np.max(np.abs(exact_solution(x, y, z))) + rel_err_rhs = np.max(np.abs(rhs - rhs_exact)) / np.max(np.abs(rhs_exact)) + rel_err_phi = np.max(np.abs(phi - phi_exact)) / np.max(np.abs(phi_exact)) print(f"Max relative error in RHS: {rel_err_rhs:.2e}") print(f"Max relative error in Phi: {rel_err_phi:.2e}") - assert rel_err_rhs < 6e-3, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." - assert rel_err_phi < 5e-3, ( + assert rel_err_rhs < 1.7e-2, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." + assert rel_err_phi < 1.2e-2, ( f"The computed solution does not match the exact solution, max rel error = {rel_err_phi}." ) diff --git a/profiling/examples/Poisson/cube_strong_scaling_massmatrix/params_poisson.py b/profiling/examples/Poisson/cube_strong_scaling_massmatrix/params_poisson.py new file mode 100644 index 000000000..a5871ffd7 --- /dev/null +++ b/profiling/examples/Poisson/cube_strong_scaling_massmatrix/params_poisson.py @@ -0,0 +1,270 @@ +# ----------------------------- +# Description of the simulation +# ----------------------------- +# Please fill in a verbal description of the simulation. +# It will be printed at the beginning of the simulation and can be used to keep track of the different runs. + +name = "Poisson strong scaling on 3D cube, mass-matrix preconditioner" +description = """ +Strong scaling test for Poisson equation on a 3D cube (Cuboid). +The manufactured solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5), +exciting the full spectrum. +Homogeneous Dirichlet boundary conditions are set in direction eta1. +The linear system is solved with CG preconditioned by the Kronecker approximation of the inverse 0-form mass matrix (MassMatrixPreconditioner). +""" + +import logging + +from struphy import set_logging_level + +set_logging_level(logging.INFO) + +import argparse + +# ------------------ +# Import Struphy API +# ------------------ +from struphy import ( + BaseUnits, + DerhamOptions, + EnvironmentOptions, + FieldsBackground, + ProfilingOptions, + Simulation, + Time, + domains, + equils, + grids, + perturbations, +) + +# --------------------- +# Instance of the model +# --------------------- +from struphy.models import Poisson + +# Units +base_units = BaseUnits() + +# Model instance +model = Poisson(base_units=base_units) + +# List all variables and decide whether to save their data +model.em_fields.phi.save_data = True +model.em_fields.source.save_data = True + +# -------------------------- +# Instance of the simulation +# -------------------------- + +# Environment options +# `--id` distinguishes runs that share a rank count but differ in something else; the +# profiling driver passes its launch counter (see `ProfilingJob.build_commands`). +# Unknown flags are ignored so the driver can forward other parameters as well. +parser = argparse.ArgumentParser() +parser.add_argument("--id", type=int, default=0, help="Run id, used to name the output folder.") +args, _ = parser.parse_known_args() + +env = EnvironmentOptions( + sim_folder=f"sim_{args.id:02d}", + restart=False, +) + +# Time stepping +time_opts = Time() + +# Geometry +Lx = 2.0 +Ly = 3.0 +Lz = 4.0 +domain = domains.Cuboid(r1=Lx, l2=-Ly / 2, r2=Ly / 2, r3=Lz) + +# Fluid equilibrium (can be used as part of initial conditions) +equil = None + +# Grid +grid = grids.TensorProductGrid(num_elements=(64, 64, 64), mpi_dims_mask=(True, True, True)) + +# Derham options +derham_opts = DerhamOptions(degree=(1, 2, 3), bcs=(("dirichlet", "dirichlet"), None, None)) + +# Profilinig options +profiling_opts = ProfilingOptions( + use_line_profiler=True, +) + +# Simulation object +sim = Simulation( + model=model, + name=name, + description=description, + params_path=__file__, + env=env, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + profiling_opts=profiling_opts, +) + +# ------------------ +# Propagator options +# ------------------ + +from struphy.linear_algebra.solver import SolverParameters + +solver_params = SolverParameters(tol=1e-8, maxiter=3000, info=True, recycle=False) +model.propagators.poisson.options = model.propagators.poisson.Options( + stab_eps=0.0, + solver="pcg", + precond="MassMatrixPreconditioner", + solver_params=solver_params, +) + +# ------------------ +# Initial conditions +# ------------------ +import numpy as np + +from struphy.initial.base import GenericPerturbation + +# The exact solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5). +# Its value at the Dirichlet boundary eta1 = 0, 1 is exp(-0.25 / w**2) ~ 1e-11, and it is periodic +# in eta2 and eta3 up to the same accuracy. Being localized, it excites the whole spectrum of the +# Laplacian, so CG has to work for convergence. Both the solution and the right-hand side +# -Laplace(phi) are 0-forms given as functions of the logical coordinates (given_in_basis="0"); +# the pushed-forward solution at x = F(eta) is phi(eta). +w = 0.1 +eta0 = (0.5, 0.5, 0.5) + + +def gaussian_derivatives(e1, e2, e3): + """The Gaussian and its first and second derivatives with respect to eta1, eta2, eta3.""" + eta = (e1, e2, e3) + g = np.exp(-sum((e - c) ** 2 for e, c in zip(eta, eta0)) / w**2) + d = [-2 * (e - c) / w**2 * g for e, c in zip(eta, eta0)] + dd = [(4 * (e - c) ** 2 / w**4 - 2 / w**2) * g for e, c in zip(eta, eta0)] + return g, d, dd + + +def exact_solution(e1, e2, e3): + return gaussian_derivatives(e1, e2, e3)[0] + + +def laplacian(e1, e2, e3): + """Physical Laplacian of the exact solution, in logical coordinates (Cuboid: x_i = l_i + L_i * eta_i).""" + g, d, dd = gaussian_derivatives(e1, e2, e3) + return dd[0] / Lx**2 + dd[1] / Ly**2 + dd[2] / Lz**2 + + +def rhs_fun(e1, e2, e3): + return -laplacian(e1, e2, e3) + + +rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="0") + +model.em_fields.source.add_perturbation(rhs_perturbation) + + +if __name__ == "__main__": + # switches for debugging and testing + estimate_mem = False + run = True + pproc = True + save_figs = True + + if estimate_mem: + sim.estimate_mem(print_report=True) + exit() + + if run: + out = sim.run(profiling_activated=True, one_time_step=True) + else: + # use this for working with existing sim data from a previous run + out = sim.output + + if pproc: + out.pproc(create_vtk=True, parallel=True) + + from matplotlib import pyplot as plt + + def slice_axes(k, X, Y, Z, E1, E2, E3): + """Axes (horizontal, vertical, labels) for plotting the slice eta_k = 0.5 in physical coordinates.""" + if k == 0: + return Y, Z, "y", "z" + elif k == 1: + return X, Z, "x", "z" + else: + return X, Y, "x", "y" + + def plot_slices(num, exact, name): + """Numerical solution, exact solution and error on the slices eta_k = 0.5 (k = 1, 2, 3).""" + fig = plt.figure(figsize=(16, 12)) + for k in range(3): + idx = tuple(num.shape[k] // 2 if i == k else slice(None) for i in range(3)) + h, v, hlabel, vlabel = slice_axes(k, X[idx], Y[idx], Z[idx], E1[idx], E2[idx], E3[idx]) + eta_k = (E1, E2, E3)[k][idx].flat[0] + for row, (data, what) in enumerate( + [(num[idx], "from struphy"), (exact[idx], "exact"), (np.abs(num[idx] - exact[idx]), "error")] + ): + plt.subplot(3, 3, 3 * row + k + 1) + plt.pcolormesh(h, v, data, shading="gouraud") + plt.colorbar() + plt.xlabel(hlabel) + plt.ylabel(vlabel) + plt.title(f"{name} {what}, slice at eta{k + 1} = {eta_k:.2f}") + fig.tight_layout() + return fig + + if sim.rank == 0: + # Raw FEEC fields are evaluated directly from the saved spline coefficients at the + # cell centres of the simulation grid (serial Derham, safe on rank 0 only). + rhs = out.evaluate("em_fields/source", t=0).values + + phi_data = out.evaluate("em_fields/phi", t=-1) + phi = phi_data.values + + # logical coordinates of the evaluation points and their images x = F(eta) + E1, E2, E3 = np.meshgrid( + *(phi_data.coords[c].values for c in ("eta1", "eta2", "eta3")), + indexing="ij", + ) + X, Y, Z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) + + # pushed-forward exact solutions: the 0-forms at x = F(eta) are their values at eta + rhs_exact = rhs_fun(E1, E2, E3) + phi_exact = exact_solution(E1, E2, E3) + + fig_rhs = plot_slices(rhs, rhs_exact, "RHS") + fig_phi = plot_slices(phi, phi_exact, "Phi") + + rel_err_rhs = np.max(np.abs(rhs - rhs_exact)) / np.max(np.abs(rhs_exact)) + rel_err_phi = np.max(np.abs(phi - phi_exact)) / np.max(np.abs(phi_exact)) + + print(f"Max relative error in RHS: {rel_err_rhs:.2e}") + print(f"Max relative error in Phi: {rel_err_phi:.2e}") + + assert rel_err_rhs < 1.7e-2, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." + assert rel_err_phi < 1.2e-2, ( + f"The computed solution does not match the exact solution, max rel error = {rel_err_phi}." + ) + + import os + + # `out.path_out` is the out's (absolute) output folder; `sim_folder` alone is a bare + # name resolved against the CWD. The profiling packaging picks these files up from + # here and uploads them as `results-out`. + results_dir = os.path.join(out.path_out, "results") + os.makedirs(results_dir, exist_ok=True) + + np.save(os.path.join(results_dir, "rel_err_rhs.npy"), rel_err_rhs) + np.save(os.path.join(results_dir, "rel_err_phi.npy"), rel_err_phi) + np.save(os.path.join(results_dir, "resolution.npy"), sim.grid.num_elements) + np.save(os.path.join(results_dir, "spline_degree.npy"), sim.derham_opts.degree) + + if save_figs: + fig_rhs.savefig(os.path.join(results_dir, "rhs_slices.png")) + fig_phi.savefig(os.path.join(results_dir, "phi_slices.png")) + else: + plt.show() diff --git a/profiling/examples/Poisson/cube_strong_scaling_multigrid/params_poisson.py b/profiling/examples/Poisson/cube_strong_scaling_multigrid/params_poisson.py index a60554827..ff030f266 100644 --- a/profiling/examples/Poisson/cube_strong_scaling_multigrid/params_poisson.py +++ b/profiling/examples/Poisson/cube_strong_scaling_multigrid/params_poisson.py @@ -6,9 +6,10 @@ name = "Poisson strong scaling on 3D cube, multigrid preconditioner" description = """ -Strong scaling test for Poisson equation on a 3D cube. -The manufactured solution is a Gaussian bump times sin(pi*x/Lx), exciting the full spectrum. -Homogeneous Dirichlet boundary conditions are set in direction x. +Strong scaling test for Poisson equation on a 3D cube (Cuboid). +The manufactured solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5), +exciting the full spectrum. +Homogeneous Dirichlet boundary conditions are set in direction eta1. The linear system is solved with CG preconditioned by geometric multigrid. """ @@ -114,7 +115,7 @@ from struphy.linear_algebra.solver import SolverParameters -solver_params = SolverParameters(tol=1e-8, maxiter=100, info=True, recycle=False) +solver_params = SolverParameters(tol=1e-8, maxiter=3000, info=True, recycle=False) model.propagators.poisson.options = model.propagators.poisson.Options( stab_eps=0.0, solver="pcg", @@ -130,30 +131,40 @@ from struphy.initial.base import GenericPerturbation -# Gaussian bump times sin(pi*x/Lx) (exact Dirichlet in x); the bump is narrow enough -# to be periodic in y and z up to exp(-(Ly/2)**2/w**2) ~ 1e-11. Being localized, it -# excites the whole spectrum of the Laplacian, so CG has to work for convergence. -w = 0.3 -x0, y0, z0 = Lx / 2, 0.0, Lz / 2 +# The exact solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5). +# Its value at the Dirichlet boundary eta1 = 0, 1 is exp(-0.25 / w**2) ~ 1e-11, and it is periodic +# in eta2 and eta3 up to the same accuracy. Being localized, it excites the whole spectrum of the +# Laplacian, so CG has to work for convergence. Both the solution and the right-hand side +# -Laplace(phi) are 0-forms given as functions of the logical coordinates (given_in_basis="0"); +# the pushed-forward solution at x = F(eta) is phi(eta). +w = 0.1 +eta0 = (0.5, 0.5, 0.5) -def exact_solution(x, y, z): - gauss = np.exp(-((x - x0) ** 2 + (y - y0) ** 2 + (z - z0) ** 2) / w**2) - return np.sin(np.pi / Lx * x) * gauss +def gaussian_derivatives(e1, e2, e3): + """The Gaussian and its first and second derivatives with respect to eta1, eta2, eta3.""" + eta = (e1, e2, e3) + g = np.exp(-sum((e - c) ** 2 for e, c in zip(eta, eta0)) / w**2) + d = [-2 * (e - c) / w**2 * g for e, c in zip(eta, eta0)] + dd = [(4 * (e - c) ** 2 / w**4 - 2 / w**2) * g for e, c in zip(eta, eta0)] + return g, d, dd -def rhs_fun(x, y, z): - # -Laplace(s*g) = -(s'' g + 2 s' dg/dx + s Laplace(g)) - r2 = (x - x0) ** 2 + (y - y0) ** 2 + (z - z0) ** 2 - gauss = np.exp(-r2 / w**2) - s = np.sin(np.pi / Lx * x) - ds = np.pi / Lx * np.cos(np.pi / Lx * x) - lap_gauss = (4 * r2 / w**4 - 6 / w**2) * gauss - dgauss_dx = -2 * (x - x0) / w**2 * gauss - return (np.pi / Lx) ** 2 * s * gauss - 2 * ds * dgauss_dx - s * lap_gauss +def exact_solution(e1, e2, e3): + return gaussian_derivatives(e1, e2, e3)[0] -rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="physical") +def laplacian(e1, e2, e3): + """Physical Laplacian of the exact solution, in logical coordinates (Cuboid: x_i = l_i + L_i * eta_i).""" + g, d, dd = gaussian_derivatives(e1, e2, e3) + return dd[0] / Lx**2 + dd[1] / Ly**2 + dd[2] / Lz**2 + + +def rhs_fun(e1, e2, e3): + return -laplacian(e1, e2, e3) + + +rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="0") model.em_fields.source.add_perturbation(rhs_perturbation) @@ -180,126 +191,64 @@ def rhs_fun(x, y, z): from matplotlib import pyplot as plt - def plot_slices(num, exact, name, slice_pt_x=0, slice_pt_y=0, slice_pt_z=0): + def slice_axes(k, X, Y, Z, E1, E2, E3): + """Axes (horizontal, vertical, labels) for plotting the slice eta_k = 0.5 in physical coordinates.""" + if k == 0: + return Y, Z, "y", "z" + elif k == 1: + return X, Z, "x", "z" + else: + return X, Y, "x", "y" + def plot_slices(num, exact, name): + """Numerical solution, exact solution and error on the slices eta_k = 0.5 (k = 1, 2, 3).""" fig = plt.figure(figsize=(16, 12)) - - plt.subplot(3, 3, 1) - plt.pcolor(y[slice_pt_x, :, :], z[slice_pt_x, :, :], num[slice_pt_x, :, :]) - plt.colorbar() - plt.xlabel("y") - plt.ylabel("z") - plt.title("{} from struphy, slice at x = {:.2f}".format(name, x[slice_pt_x, 0, 0])) - - plt.subplot(3, 3, 4) - plt.pcolor( - y[slice_pt_x, :, :], - z[slice_pt_x, :, :], - exact(x[slice_pt_x, :, :], y[slice_pt_x, :, :], z[slice_pt_x, :, :]), - ) - plt.colorbar() - plt.xlabel("y") - plt.ylabel("z") - plt.title("{} exact, slice at x = {:.2f}".format(name, x[slice_pt_x, 0, 0])) - - plt.subplot(3, 3, 7) - plt.pcolor( - y[slice_pt_x, :, :], - z[slice_pt_x, :, :], - np.abs(num[slice_pt_x, :, :] - exact(x[slice_pt_x, :, :], y[slice_pt_x, :, :], z[slice_pt_x, :, :])), - ) - plt.colorbar() - plt.xlabel("y") - plt.ylabel("z") - plt.title("{} error, slice at x = {:.2f}".format(name, x[slice_pt_x, 0, 0])) - - plt.subplot(3, 3, 2) - plt.pcolor(x[:, slice_pt_y, :], z[:, slice_pt_y, :], num[:, slice_pt_y, :]) - plt.colorbar() - plt.xlabel("x") - plt.ylabel("z") - plt.title("{} from struphy, slice at y = {:.2f}".format(name, y[0, slice_pt_y, 0])) - - plt.subplot(3, 3, 5) - plt.pcolor( - x[:, slice_pt_y, :], - z[:, slice_pt_y, :], - exact(x[:, slice_pt_y, :], y[:, slice_pt_y, :], z[:, slice_pt_y, :]), - ) - plt.colorbar() - plt.xlabel("x") - plt.ylabel("z") - plt.title("{} exact, slice at y = {:.2f}".format(name, y[0, slice_pt_y, 0])) - - plt.subplot(3, 3, 8) - plt.pcolor( - x[:, slice_pt_y, :], - z[:, slice_pt_y, :], - np.abs(num[:, slice_pt_y, :] - exact(x[:, slice_pt_y, :], y[:, slice_pt_y, :], z[:, slice_pt_y, :])), - ) - plt.colorbar() - plt.xlabel("x") - plt.ylabel("z") - plt.title("{} error, slice at y = {:.2f}".format(name, y[0, slice_pt_y, 0])) - - plt.subplot(3, 3, 3) - plt.pcolor(x[:, :, slice_pt_z], y[:, :, slice_pt_z], num[:, :, slice_pt_z]) - plt.colorbar() - plt.xlabel("x") - plt.ylabel("y") - plt.title("{} from struphy, slice at z = {:.2f}".format(name, z[0, 0, slice_pt_z])) - - plt.subplot(3, 3, 6) - plt.pcolor( - x[:, :, slice_pt_z], - y[:, :, slice_pt_z], - exact(x[:, :, slice_pt_z], y[:, :, slice_pt_z], z[:, :, slice_pt_z]), - ) - plt.colorbar() - plt.xlabel("x") - plt.ylabel("y") - plt.title("{} exact, slice at z = {:.2f}".format(name, z[0, 0, slice_pt_z])) - - plt.subplot(3, 3, 9) - plt.pcolor( - x[:, :, slice_pt_z], - y[:, :, slice_pt_z], - np.abs(num[:, :, slice_pt_z] - exact(x[:, :, slice_pt_z], y[:, :, slice_pt_z], z[:, :, slice_pt_z])), - ) - plt.colorbar() - plt.xlabel("x") - plt.ylabel("y") - plt.title("{} error, slice at z = {:.2f}".format(name, z[0, 0, slice_pt_z])) - + for k in range(3): + idx = tuple(num.shape[k] // 2 if i == k else slice(None) for i in range(3)) + h, v, hlabel, vlabel = slice_axes(k, X[idx], Y[idx], Z[idx], E1[idx], E2[idx], E3[idx]) + eta_k = (E1, E2, E3)[k][idx].flat[0] + for row, (data, what) in enumerate( + [(num[idx], "from struphy"), (exact[idx], "exact"), (np.abs(num[idx] - exact[idx]), "error")] + ): + plt.subplot(3, 3, 3 * row + k + 1) + plt.pcolormesh(h, v, data, shading="gouraud") + plt.colorbar() + plt.xlabel(hlabel) + plt.ylabel(vlabel) + plt.title(f"{name} {what}, slice at eta{k + 1} = {eta_k:.2f}") + fig.tight_layout() return fig if sim.rank == 0: # Raw FEEC fields are evaluated directly from the saved spline coefficients at the # cell centres of the simulation grid (serial Derham, safe on rank 0 only). - rhs_data = out.evaluate("em_fields/source", t=0) - rhs = rhs_data.values + rhs = out.evaluate("em_fields/source", t=0).values phi_data = out.evaluate("em_fields/phi", t=-1) phi = phi_data.values - x, y, z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) - - slice_pt_x = x.shape[0] // 2 - slice_pt_y = y.shape[1] // 2 - slice_pt_z = 0 - fig_rhs = plot_slices(rhs, rhs_fun, "RHS", slice_pt_x=slice_pt_x, slice_pt_y=slice_pt_y, slice_pt_z=slice_pt_z) - fig_phi = plot_slices( - phi, exact_solution, "Phi", slice_pt_x=slice_pt_x, slice_pt_y=slice_pt_y, slice_pt_z=slice_pt_z + # logical coordinates of the evaluation points and their images x = F(eta) + E1, E2, E3 = np.meshgrid( + *(phi_data.coords[c].values for c in ("eta1", "eta2", "eta3")), + indexing="ij", ) + X, Y, Z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) + + # pushed-forward exact solutions: the 0-forms at x = F(eta) are their values at eta + rhs_exact = rhs_fun(E1, E2, E3) + phi_exact = exact_solution(E1, E2, E3) + + fig_rhs = plot_slices(rhs, rhs_exact, "RHS") + fig_phi = plot_slices(phi, phi_exact, "Phi") - rel_err_rhs = np.max(np.abs(rhs - rhs_fun(x, y, z))) / np.max(np.abs(rhs_fun(x, y, z))) - rel_err_phi = np.max(np.abs(phi - exact_solution(x, y, z))) / np.max(np.abs(exact_solution(x, y, z))) + rel_err_rhs = np.max(np.abs(rhs - rhs_exact)) / np.max(np.abs(rhs_exact)) + rel_err_phi = np.max(np.abs(phi - phi_exact)) / np.max(np.abs(phi_exact)) print(f"Max relative error in RHS: {rel_err_rhs:.2e}") print(f"Max relative error in Phi: {rel_err_phi:.2e}") - assert rel_err_rhs < 6e-3, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." - assert rel_err_phi < 5e-3, ( + assert rel_err_rhs < 1.7e-2, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." + assert rel_err_phi < 1.2e-2, ( f"The computed solution does not match the exact solution, max rel error = {rel_err_phi}." ) diff --git a/profiling/examples/Poisson/cube_strong_scaling_stiffness/params_poisson.py b/profiling/examples/Poisson/cube_strong_scaling_stiffness/params_poisson.py new file mode 100644 index 000000000..4fe500f7e --- /dev/null +++ b/profiling/examples/Poisson/cube_strong_scaling_stiffness/params_poisson.py @@ -0,0 +1,270 @@ +# ----------------------------- +# Description of the simulation +# ----------------------------- +# Please fill in a verbal description of the simulation. +# It will be printed at the beginning of the simulation and can be used to keep track of the different runs. + +name = "Poisson strong scaling on 3D cube, stiffness preconditioner" +description = """ +Strong scaling test for Poisson equation on a 3D cube (Cuboid). +The manufactured solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5), +exciting the full spectrum. +Homogeneous Dirichlet boundary conditions are set in direction eta1. +The linear system is solved with CG preconditioned by the Kronecker approximation of the inverse Laplacian (StiffnessPreconditioner). +""" + +import logging + +from struphy import set_logging_level + +set_logging_level(logging.INFO) + +import argparse + +# ------------------ +# Import Struphy API +# ------------------ +from struphy import ( + BaseUnits, + DerhamOptions, + EnvironmentOptions, + FieldsBackground, + ProfilingOptions, + Simulation, + Time, + domains, + equils, + grids, + perturbations, +) + +# --------------------- +# Instance of the model +# --------------------- +from struphy.models import Poisson + +# Units +base_units = BaseUnits() + +# Model instance +model = Poisson(base_units=base_units) + +# List all variables and decide whether to save their data +model.em_fields.phi.save_data = True +model.em_fields.source.save_data = True + +# -------------------------- +# Instance of the simulation +# -------------------------- + +# Environment options +# `--id` distinguishes runs that share a rank count but differ in something else; the +# profiling driver passes its launch counter (see `ProfilingJob.build_commands`). +# Unknown flags are ignored so the driver can forward other parameters as well. +parser = argparse.ArgumentParser() +parser.add_argument("--id", type=int, default=0, help="Run id, used to name the output folder.") +args, _ = parser.parse_known_args() + +env = EnvironmentOptions( + sim_folder=f"sim_{args.id:02d}", + restart=False, +) + +# Time stepping +time_opts = Time() + +# Geometry +Lx = 2.0 +Ly = 3.0 +Lz = 4.0 +domain = domains.Cuboid(r1=Lx, l2=-Ly / 2, r2=Ly / 2, r3=Lz) + +# Fluid equilibrium (can be used as part of initial conditions) +equil = None + +# Grid +grid = grids.TensorProductGrid(num_elements=(64, 64, 64), mpi_dims_mask=(True, True, True)) + +# Derham options +derham_opts = DerhamOptions(degree=(1, 2, 3), bcs=(("dirichlet", "dirichlet"), None, None)) + +# Profilinig options +profiling_opts = ProfilingOptions( + use_line_profiler=True, +) + +# Simulation object +sim = Simulation( + model=model, + name=name, + description=description, + params_path=__file__, + env=env, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + profiling_opts=profiling_opts, +) + +# ------------------ +# Propagator options +# ------------------ + +from struphy.linear_algebra.solver import SolverParameters + +solver_params = SolverParameters(tol=1e-8, maxiter=3000, info=True, recycle=False) +model.propagators.poisson.options = model.propagators.poisson.Options( + stab_eps=0.0, + solver="pcg", + precond="StiffnessPreconditioner", + solver_params=solver_params, +) + +# ------------------ +# Initial conditions +# ------------------ +import numpy as np + +from struphy.initial.base import GenericPerturbation + +# The exact solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5). +# Its value at the Dirichlet boundary eta1 = 0, 1 is exp(-0.25 / w**2) ~ 1e-11, and it is periodic +# in eta2 and eta3 up to the same accuracy. Being localized, it excites the whole spectrum of the +# Laplacian, so CG has to work for convergence. Both the solution and the right-hand side +# -Laplace(phi) are 0-forms given as functions of the logical coordinates (given_in_basis="0"); +# the pushed-forward solution at x = F(eta) is phi(eta). +w = 0.1 +eta0 = (0.5, 0.5, 0.5) + + +def gaussian_derivatives(e1, e2, e3): + """The Gaussian and its first and second derivatives with respect to eta1, eta2, eta3.""" + eta = (e1, e2, e3) + g = np.exp(-sum((e - c) ** 2 for e, c in zip(eta, eta0)) / w**2) + d = [-2 * (e - c) / w**2 * g for e, c in zip(eta, eta0)] + dd = [(4 * (e - c) ** 2 / w**4 - 2 / w**2) * g for e, c in zip(eta, eta0)] + return g, d, dd + + +def exact_solution(e1, e2, e3): + return gaussian_derivatives(e1, e2, e3)[0] + + +def laplacian(e1, e2, e3): + """Physical Laplacian of the exact solution, in logical coordinates (Cuboid: x_i = l_i + L_i * eta_i).""" + g, d, dd = gaussian_derivatives(e1, e2, e3) + return dd[0] / Lx**2 + dd[1] / Ly**2 + dd[2] / Lz**2 + + +def rhs_fun(e1, e2, e3): + return -laplacian(e1, e2, e3) + + +rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="0") + +model.em_fields.source.add_perturbation(rhs_perturbation) + + +if __name__ == "__main__": + # switches for debugging and testing + estimate_mem = False + run = True + pproc = True + save_figs = True + + if estimate_mem: + sim.estimate_mem(print_report=True) + exit() + + if run: + out = sim.run(profiling_activated=True, one_time_step=True) + else: + # use this for working with existing sim data from a previous run + out = sim.output + + if pproc: + out.pproc(create_vtk=True, parallel=True) + + from matplotlib import pyplot as plt + + def slice_axes(k, X, Y, Z, E1, E2, E3): + """Axes (horizontal, vertical, labels) for plotting the slice eta_k = 0.5 in physical coordinates.""" + if k == 0: + return Y, Z, "y", "z" + elif k == 1: + return X, Z, "x", "z" + else: + return X, Y, "x", "y" + + def plot_slices(num, exact, name): + """Numerical solution, exact solution and error on the slices eta_k = 0.5 (k = 1, 2, 3).""" + fig = plt.figure(figsize=(16, 12)) + for k in range(3): + idx = tuple(num.shape[k] // 2 if i == k else slice(None) for i in range(3)) + h, v, hlabel, vlabel = slice_axes(k, X[idx], Y[idx], Z[idx], E1[idx], E2[idx], E3[idx]) + eta_k = (E1, E2, E3)[k][idx].flat[0] + for row, (data, what) in enumerate( + [(num[idx], "from struphy"), (exact[idx], "exact"), (np.abs(num[idx] - exact[idx]), "error")] + ): + plt.subplot(3, 3, 3 * row + k + 1) + plt.pcolormesh(h, v, data, shading="gouraud") + plt.colorbar() + plt.xlabel(hlabel) + plt.ylabel(vlabel) + plt.title(f"{name} {what}, slice at eta{k + 1} = {eta_k:.2f}") + fig.tight_layout() + return fig + + if sim.rank == 0: + # Raw FEEC fields are evaluated directly from the saved spline coefficients at the + # cell centres of the simulation grid (serial Derham, safe on rank 0 only). + rhs = out.evaluate("em_fields/source", t=0).values + + phi_data = out.evaluate("em_fields/phi", t=-1) + phi = phi_data.values + + # logical coordinates of the evaluation points and their images x = F(eta) + E1, E2, E3 = np.meshgrid( + *(phi_data.coords[c].values for c in ("eta1", "eta2", "eta3")), + indexing="ij", + ) + X, Y, Z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) + + # pushed-forward exact solutions: the 0-forms at x = F(eta) are their values at eta + rhs_exact = rhs_fun(E1, E2, E3) + phi_exact = exact_solution(E1, E2, E3) + + fig_rhs = plot_slices(rhs, rhs_exact, "RHS") + fig_phi = plot_slices(phi, phi_exact, "Phi") + + rel_err_rhs = np.max(np.abs(rhs - rhs_exact)) / np.max(np.abs(rhs_exact)) + rel_err_phi = np.max(np.abs(phi - phi_exact)) / np.max(np.abs(phi_exact)) + + print(f"Max relative error in RHS: {rel_err_rhs:.2e}") + print(f"Max relative error in Phi: {rel_err_phi:.2e}") + + assert rel_err_rhs < 1.7e-2, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." + assert rel_err_phi < 1.2e-2, ( + f"The computed solution does not match the exact solution, max rel error = {rel_err_phi}." + ) + + import os + + # `out.path_out` is the out's (absolute) output folder; `sim_folder` alone is a bare + # name resolved against the CWD. The profiling packaging picks these files up from + # here and uploads them as `results-out`. + results_dir = os.path.join(out.path_out, "results") + os.makedirs(results_dir, exist_ok=True) + + np.save(os.path.join(results_dir, "rel_err_rhs.npy"), rel_err_rhs) + np.save(os.path.join(results_dir, "rel_err_phi.npy"), rel_err_phi) + np.save(os.path.join(results_dir, "resolution.npy"), sim.grid.num_elements) + np.save(os.path.join(results_dir, "spline_degree.npy"), sim.derham_opts.degree) + + if save_figs: + fig_rhs.savefig(os.path.join(results_dir, "rhs_slices.png")) + fig_phi.savefig(os.path.join(results_dir, "phi_slices.png")) + else: + plt.show() diff --git a/profiling/examples/Poisson/hollow_cylinder_strong_scaling/params_poisson.py b/profiling/examples/Poisson/hollow_cylinder_strong_scaling/params_poisson.py new file mode 100644 index 000000000..4af93b5d5 --- /dev/null +++ b/profiling/examples/Poisson/hollow_cylinder_strong_scaling/params_poisson.py @@ -0,0 +1,277 @@ +# ----------------------------- +# Description of the simulation +# ----------------------------- +# Please fill in a verbal description of the simulation. +# It will be printed at the beginning of the simulation and can be used to keep track of the different runs. + +name = "Poisson strong scaling on hollow cylinder, no preconditioner" +description = """ +Strong scaling test for Poisson equation on a hollow cylinder (HollowCylinder, a1=1e-2, a2=1, Lz=4). +The manufactured solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5), +exciting the full spectrum. +Homogeneous Dirichlet boundary conditions are set in direction eta1. +The linear system is solved with unpreconditioned CG. +""" + +import logging + +from struphy import set_logging_level + +set_logging_level(logging.INFO) + +import argparse + +# ------------------ +# Import Struphy API +# ------------------ +from struphy import ( + BaseUnits, + DerhamOptions, + EnvironmentOptions, + FieldsBackground, + ProfilingOptions, + Simulation, + Time, + domains, + equils, + grids, + perturbations, +) + +# --------------------- +# Instance of the model +# --------------------- +from struphy.models import Poisson + +# Units +base_units = BaseUnits() + +# Model instance +model = Poisson(base_units=base_units) + +# List all variables and decide whether to save their data +model.em_fields.phi.save_data = True +model.em_fields.source.save_data = True + +# -------------------------- +# Instance of the simulation +# -------------------------- + +# Environment options +# `--id` distinguishes runs that share a rank count but differ in something else; the +# profiling driver passes its launch counter (see `ProfilingJob.build_commands`). +# Unknown flags are ignored so the driver can forward other parameters as well. +parser = argparse.ArgumentParser() +parser.add_argument("--id", type=int, default=0, help="Run id, used to name the output folder.") +args, _ = parser.parse_known_args() + +env = EnvironmentOptions( + sim_folder=f"sim_{args.id:02d}", + restart=False, +) + +# Time stepping +time_opts = Time() + +# Geometry +a1 = 1e-2 +a2 = 1.0 +Lz = 4.0 +domain = domains.HollowCylinder(a1=a1, a2=a2, Lz=Lz) + +# Fluid equilibrium (can be used as part of initial conditions) +equil = None + +# Grid +grid = grids.TensorProductGrid(num_elements=(64, 64, 64), mpi_dims_mask=(True, True, True)) + +# Derham options +derham_opts = DerhamOptions(degree=(1, 2, 3), bcs=(("dirichlet", "dirichlet"), None, None)) + +# Profilinig options +profiling_opts = ProfilingOptions( + use_line_profiler=True, +) + +# Simulation object +sim = Simulation( + model=model, + name=name, + description=description, + params_path=__file__, + env=env, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + profiling_opts=profiling_opts, +) + +# ------------------ +# Propagator options +# ------------------ + +from struphy.linear_algebra.solver import SolverParameters + +solver_params = SolverParameters(tol=1e-8, maxiter=3000, info=True, recycle=False) +model.propagators.poisson.options = model.propagators.poisson.Options( + stab_eps=0.0, + solver="pcg", + precond=None, + solver_params=solver_params, +) + +# ------------------ +# Initial conditions +# ------------------ +import numpy as np + +from struphy.initial.base import GenericPerturbation + +# The exact solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5). +# Its value at the Dirichlet boundary eta1 = 0, 1 is exp(-0.25 / w**2) ~ 1e-11, and it is periodic +# in eta2 and eta3 up to the same accuracy. Being localized, it excites the whole spectrum of the +# Laplacian, so CG has to work for convergence. Both the solution and the right-hand side +# -Laplace(phi) are 0-forms given as functions of the logical coordinates (given_in_basis="0"); +# the pushed-forward solution at x = F(eta) is phi(eta). +w = 0.1 +eta0 = (0.5, 0.5, 0.5) + + +def gaussian_derivatives(e1, e2, e3): + """The Gaussian and its first and second derivatives with respect to eta1, eta2, eta3.""" + eta = (e1, e2, e3) + g = np.exp(-sum((e - c) ** 2 for e, c in zip(eta, eta0)) / w**2) + d = [-2 * (e - c) / w**2 * g for e, c in zip(eta, eta0)] + dd = [(4 * (e - c) ** 2 / w**4 - 2 / w**2) * g for e, c in zip(eta, eta0)] + return g, d, dd + + +def exact_solution(e1, e2, e3): + return gaussian_derivatives(e1, e2, e3)[0] + + +def laplacian(e1, e2, e3): + """Physical Laplacian of the exact solution, in logical coordinates. + + HollowCylinder: r = a1 + (a2 - a1) * eta1, theta = 2*pi*eta2, z = Lz * eta3, and + Laplace(phi) = phi_rr + phi_r / r + phi_theta,theta / r**2 + phi_zz. + """ + g, d, dd = gaussian_derivatives(e1, e2, e3) + da = a2 - a1 + r = a1 + da * e1 + return dd[0] / da**2 + d[0] / (da * r) + dd[1] / (2 * np.pi * r) ** 2 + dd[2] / Lz**2 + + +def rhs_fun(e1, e2, e3): + return -laplacian(e1, e2, e3) + + +rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="0") + +model.em_fields.source.add_perturbation(rhs_perturbation) + + +if __name__ == "__main__": + # switches for debugging and testing + estimate_mem = False + run = True + pproc = True + save_figs = True + + if estimate_mem: + sim.estimate_mem(print_report=True) + exit() + + if run: + out = sim.run(profiling_activated=True, one_time_step=True) + else: + # use this for working with existing sim data from a previous run + out = sim.output + + if pproc: + out.pproc(create_vtk=True, parallel=True) + + from matplotlib import pyplot as plt + + def slice_axes(k, X, Y, Z, E1, E2, E3): + """Axes (horizontal, vertical, labels) for plotting the slice eta_k = 0.5: the cross-section (x, y) + at eta3 = 0.5, the half plane (r, z) at eta2 = 0.5 and the logical (eta2, eta3) at eta1 = 0.5.""" + if k == 0: + return E2, E3, "eta2", "eta3" + elif k == 1: + return np.sqrt(X**2 + Y**2), Z, "r", "z" + else: + return X, Y, "x", "y" + + def plot_slices(num, exact, name): + """Numerical solution, exact solution and error on the slices eta_k = 0.5 (k = 1, 2, 3).""" + fig = plt.figure(figsize=(16, 12)) + for k in range(3): + idx = tuple(num.shape[k] // 2 if i == k else slice(None) for i in range(3)) + h, v, hlabel, vlabel = slice_axes(k, X[idx], Y[idx], Z[idx], E1[idx], E2[idx], E3[idx]) + eta_k = (E1, E2, E3)[k][idx].flat[0] + for row, (data, what) in enumerate( + [(num[idx], "from struphy"), (exact[idx], "exact"), (np.abs(num[idx] - exact[idx]), "error")] + ): + plt.subplot(3, 3, 3 * row + k + 1) + plt.pcolormesh(h, v, data, shading="gouraud") + plt.colorbar() + plt.xlabel(hlabel) + plt.ylabel(vlabel) + plt.title(f"{name} {what}, slice at eta{k + 1} = {eta_k:.2f}") + fig.tight_layout() + return fig + + if sim.rank == 0: + # Raw FEEC fields are evaluated directly from the saved spline coefficients at the + # cell centres of the simulation grid (serial Derham, safe on rank 0 only). + rhs = out.evaluate("em_fields/source", t=0).values + + phi_data = out.evaluate("em_fields/phi", t=-1) + phi = phi_data.values + + # logical coordinates of the evaluation points and their images x = F(eta) + E1, E2, E3 = np.meshgrid( + *(phi_data.coords[c].values for c in ("eta1", "eta2", "eta3")), + indexing="ij", + ) + X, Y, Z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) + + # pushed-forward exact solutions: the 0-forms at x = F(eta) are their values at eta + rhs_exact = rhs_fun(E1, E2, E3) + phi_exact = exact_solution(E1, E2, E3) + + fig_rhs = plot_slices(rhs, rhs_exact, "RHS") + fig_phi = plot_slices(phi, phi_exact, "Phi") + + rel_err_rhs = np.max(np.abs(rhs - rhs_exact)) / np.max(np.abs(rhs_exact)) + rel_err_phi = np.max(np.abs(phi - phi_exact)) / np.max(np.abs(phi_exact)) + + print(f"Max relative error in RHS: {rel_err_rhs:.2e}") + print(f"Max relative error in Phi: {rel_err_phi:.2e}") + + assert rel_err_rhs < 2.2e-2, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." + assert rel_err_phi < 1.3e-2, ( + f"The computed solution does not match the exact solution, max rel error = {rel_err_phi}." + ) + + import os + + # `out.path_out` is the out's (absolute) output folder; `sim_folder` alone is a bare + # name resolved against the CWD. The profiling packaging picks these files up from + # here and uploads them as `results-out`. + results_dir = os.path.join(out.path_out, "results") + os.makedirs(results_dir, exist_ok=True) + + np.save(os.path.join(results_dir, "rel_err_rhs.npy"), rel_err_rhs) + np.save(os.path.join(results_dir, "rel_err_phi.npy"), rel_err_phi) + np.save(os.path.join(results_dir, "resolution.npy"), sim.grid.num_elements) + np.save(os.path.join(results_dir, "spline_degree.npy"), sim.derham_opts.degree) + + if save_figs: + fig_rhs.savefig(os.path.join(results_dir, "rhs_slices.png")) + fig_phi.savefig(os.path.join(results_dir, "phi_slices.png")) + else: + plt.show() diff --git a/profiling/examples/Poisson/hollow_cylinder_strong_scaling_massmatrix/params_poisson.py b/profiling/examples/Poisson/hollow_cylinder_strong_scaling_massmatrix/params_poisson.py new file mode 100644 index 000000000..c0fb9b8c8 --- /dev/null +++ b/profiling/examples/Poisson/hollow_cylinder_strong_scaling_massmatrix/params_poisson.py @@ -0,0 +1,277 @@ +# ----------------------------- +# Description of the simulation +# ----------------------------- +# Please fill in a verbal description of the simulation. +# It will be printed at the beginning of the simulation and can be used to keep track of the different runs. + +name = "Poisson strong scaling on hollow cylinder, mass-matrix preconditioner" +description = """ +Strong scaling test for Poisson equation on a hollow cylinder (HollowCylinder, a1=1e-2, a2=1, Lz=4). +The manufactured solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5), +exciting the full spectrum. +Homogeneous Dirichlet boundary conditions are set in direction eta1. +The linear system is solved with CG preconditioned by the Kronecker approximation of the inverse 0-form mass matrix (MassMatrixPreconditioner). +""" + +import logging + +from struphy import set_logging_level + +set_logging_level(logging.INFO) + +import argparse + +# ------------------ +# Import Struphy API +# ------------------ +from struphy import ( + BaseUnits, + DerhamOptions, + EnvironmentOptions, + FieldsBackground, + ProfilingOptions, + Simulation, + Time, + domains, + equils, + grids, + perturbations, +) + +# --------------------- +# Instance of the model +# --------------------- +from struphy.models import Poisson + +# Units +base_units = BaseUnits() + +# Model instance +model = Poisson(base_units=base_units) + +# List all variables and decide whether to save their data +model.em_fields.phi.save_data = True +model.em_fields.source.save_data = True + +# -------------------------- +# Instance of the simulation +# -------------------------- + +# Environment options +# `--id` distinguishes runs that share a rank count but differ in something else; the +# profiling driver passes its launch counter (see `ProfilingJob.build_commands`). +# Unknown flags are ignored so the driver can forward other parameters as well. +parser = argparse.ArgumentParser() +parser.add_argument("--id", type=int, default=0, help="Run id, used to name the output folder.") +args, _ = parser.parse_known_args() + +env = EnvironmentOptions( + sim_folder=f"sim_{args.id:02d}", + restart=False, +) + +# Time stepping +time_opts = Time() + +# Geometry +a1 = 1e-2 +a2 = 1.0 +Lz = 4.0 +domain = domains.HollowCylinder(a1=a1, a2=a2, Lz=Lz) + +# Fluid equilibrium (can be used as part of initial conditions) +equil = None + +# Grid +grid = grids.TensorProductGrid(num_elements=(64, 64, 64), mpi_dims_mask=(True, True, True)) + +# Derham options +derham_opts = DerhamOptions(degree=(1, 2, 3), bcs=(("dirichlet", "dirichlet"), None, None)) + +# Profilinig options +profiling_opts = ProfilingOptions( + use_line_profiler=True, +) + +# Simulation object +sim = Simulation( + model=model, + name=name, + description=description, + params_path=__file__, + env=env, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + profiling_opts=profiling_opts, +) + +# ------------------ +# Propagator options +# ------------------ + +from struphy.linear_algebra.solver import SolverParameters + +solver_params = SolverParameters(tol=1e-8, maxiter=3000, info=True, recycle=False) +model.propagators.poisson.options = model.propagators.poisson.Options( + stab_eps=0.0, + solver="pcg", + precond="MassMatrixPreconditioner", + solver_params=solver_params, +) + +# ------------------ +# Initial conditions +# ------------------ +import numpy as np + +from struphy.initial.base import GenericPerturbation + +# The exact solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5). +# Its value at the Dirichlet boundary eta1 = 0, 1 is exp(-0.25 / w**2) ~ 1e-11, and it is periodic +# in eta2 and eta3 up to the same accuracy. Being localized, it excites the whole spectrum of the +# Laplacian, so CG has to work for convergence. Both the solution and the right-hand side +# -Laplace(phi) are 0-forms given as functions of the logical coordinates (given_in_basis="0"); +# the pushed-forward solution at x = F(eta) is phi(eta). +w = 0.1 +eta0 = (0.5, 0.5, 0.5) + + +def gaussian_derivatives(e1, e2, e3): + """The Gaussian and its first and second derivatives with respect to eta1, eta2, eta3.""" + eta = (e1, e2, e3) + g = np.exp(-sum((e - c) ** 2 for e, c in zip(eta, eta0)) / w**2) + d = [-2 * (e - c) / w**2 * g for e, c in zip(eta, eta0)] + dd = [(4 * (e - c) ** 2 / w**4 - 2 / w**2) * g for e, c in zip(eta, eta0)] + return g, d, dd + + +def exact_solution(e1, e2, e3): + return gaussian_derivatives(e1, e2, e3)[0] + + +def laplacian(e1, e2, e3): + """Physical Laplacian of the exact solution, in logical coordinates. + + HollowCylinder: r = a1 + (a2 - a1) * eta1, theta = 2*pi*eta2, z = Lz * eta3, and + Laplace(phi) = phi_rr + phi_r / r + phi_theta,theta / r**2 + phi_zz. + """ + g, d, dd = gaussian_derivatives(e1, e2, e3) + da = a2 - a1 + r = a1 + da * e1 + return dd[0] / da**2 + d[0] / (da * r) + dd[1] / (2 * np.pi * r) ** 2 + dd[2] / Lz**2 + + +def rhs_fun(e1, e2, e3): + return -laplacian(e1, e2, e3) + + +rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="0") + +model.em_fields.source.add_perturbation(rhs_perturbation) + + +if __name__ == "__main__": + # switches for debugging and testing + estimate_mem = False + run = True + pproc = True + save_figs = True + + if estimate_mem: + sim.estimate_mem(print_report=True) + exit() + + if run: + out = sim.run(profiling_activated=True, one_time_step=True) + else: + # use this for working with existing sim data from a previous run + out = sim.output + + if pproc: + out.pproc(create_vtk=True, parallel=True) + + from matplotlib import pyplot as plt + + def slice_axes(k, X, Y, Z, E1, E2, E3): + """Axes (horizontal, vertical, labels) for plotting the slice eta_k = 0.5: the cross-section (x, y) + at eta3 = 0.5, the half plane (r, z) at eta2 = 0.5 and the logical (eta2, eta3) at eta1 = 0.5.""" + if k == 0: + return E2, E3, "eta2", "eta3" + elif k == 1: + return np.sqrt(X**2 + Y**2), Z, "r", "z" + else: + return X, Y, "x", "y" + + def plot_slices(num, exact, name): + """Numerical solution, exact solution and error on the slices eta_k = 0.5 (k = 1, 2, 3).""" + fig = plt.figure(figsize=(16, 12)) + for k in range(3): + idx = tuple(num.shape[k] // 2 if i == k else slice(None) for i in range(3)) + h, v, hlabel, vlabel = slice_axes(k, X[idx], Y[idx], Z[idx], E1[idx], E2[idx], E3[idx]) + eta_k = (E1, E2, E3)[k][idx].flat[0] + for row, (data, what) in enumerate( + [(num[idx], "from struphy"), (exact[idx], "exact"), (np.abs(num[idx] - exact[idx]), "error")] + ): + plt.subplot(3, 3, 3 * row + k + 1) + plt.pcolormesh(h, v, data, shading="gouraud") + plt.colorbar() + plt.xlabel(hlabel) + plt.ylabel(vlabel) + plt.title(f"{name} {what}, slice at eta{k + 1} = {eta_k:.2f}") + fig.tight_layout() + return fig + + if sim.rank == 0: + # Raw FEEC fields are evaluated directly from the saved spline coefficients at the + # cell centres of the simulation grid (serial Derham, safe on rank 0 only). + rhs = out.evaluate("em_fields/source", t=0).values + + phi_data = out.evaluate("em_fields/phi", t=-1) + phi = phi_data.values + + # logical coordinates of the evaluation points and their images x = F(eta) + E1, E2, E3 = np.meshgrid( + *(phi_data.coords[c].values for c in ("eta1", "eta2", "eta3")), + indexing="ij", + ) + X, Y, Z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) + + # pushed-forward exact solutions: the 0-forms at x = F(eta) are their values at eta + rhs_exact = rhs_fun(E1, E2, E3) + phi_exact = exact_solution(E1, E2, E3) + + fig_rhs = plot_slices(rhs, rhs_exact, "RHS") + fig_phi = plot_slices(phi, phi_exact, "Phi") + + rel_err_rhs = np.max(np.abs(rhs - rhs_exact)) / np.max(np.abs(rhs_exact)) + rel_err_phi = np.max(np.abs(phi - phi_exact)) / np.max(np.abs(phi_exact)) + + print(f"Max relative error in RHS: {rel_err_rhs:.2e}") + print(f"Max relative error in Phi: {rel_err_phi:.2e}") + + assert rel_err_rhs < 2.2e-2, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." + assert rel_err_phi < 1.3e-2, ( + f"The computed solution does not match the exact solution, max rel error = {rel_err_phi}." + ) + + import os + + # `out.path_out` is the out's (absolute) output folder; `sim_folder` alone is a bare + # name resolved against the CWD. The profiling packaging picks these files up from + # here and uploads them as `results-out`. + results_dir = os.path.join(out.path_out, "results") + os.makedirs(results_dir, exist_ok=True) + + np.save(os.path.join(results_dir, "rel_err_rhs.npy"), rel_err_rhs) + np.save(os.path.join(results_dir, "rel_err_phi.npy"), rel_err_phi) + np.save(os.path.join(results_dir, "resolution.npy"), sim.grid.num_elements) + np.save(os.path.join(results_dir, "spline_degree.npy"), sim.derham_opts.degree) + + if save_figs: + fig_rhs.savefig(os.path.join(results_dir, "rhs_slices.png")) + fig_phi.savefig(os.path.join(results_dir, "phi_slices.png")) + else: + plt.show() diff --git a/profiling/examples/Poisson/hollow_cylinder_strong_scaling_multigrid/params_poisson.py b/profiling/examples/Poisson/hollow_cylinder_strong_scaling_multigrid/params_poisson.py new file mode 100644 index 000000000..a43e58040 --- /dev/null +++ b/profiling/examples/Poisson/hollow_cylinder_strong_scaling_multigrid/params_poisson.py @@ -0,0 +1,279 @@ +# ----------------------------- +# Description of the simulation +# ----------------------------- +# Please fill in a verbal description of the simulation. +# It will be printed at the beginning of the simulation and can be used to keep track of the different runs. + +name = "Poisson strong scaling on hollow cylinder, multigrid preconditioner" +description = """ +Strong scaling test for Poisson equation on a hollow cylinder (HollowCylinder, a1=1e-2, a2=1, Lz=4). +The manufactured solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5), +exciting the full spectrum. +Homogeneous Dirichlet boundary conditions are set in direction eta1. +The linear system is solved with CG preconditioned by geometric multigrid. +""" + +import logging + +from struphy import set_logging_level + +set_logging_level(logging.INFO) + +import argparse + +# ------------------ +# Import Struphy API +# ------------------ +from struphy import ( + BaseUnits, + DerhamOptions, + EnvironmentOptions, + FieldsBackground, + MultiGridOptions, + ProfilingOptions, + Simulation, + Time, + domains, + equils, + grids, + perturbations, +) + +# --------------------- +# Instance of the model +# --------------------- +from struphy.models import Poisson + +# Units +base_units = BaseUnits() + +# Model instance +model = Poisson(base_units=base_units) + +# List all variables and decide whether to save their data +model.em_fields.phi.save_data = True +model.em_fields.source.save_data = True + +# -------------------------- +# Instance of the simulation +# -------------------------- + +# Environment options +# `--id` distinguishes runs that share a rank count but differ in something else; the +# profiling driver passes its launch counter (see `ProfilingJob.build_commands`). +# Unknown flags are ignored so the driver can forward other parameters as well. +parser = argparse.ArgumentParser() +parser.add_argument("--id", type=int, default=0, help="Run id, used to name the output folder.") +args, _ = parser.parse_known_args() + +env = EnvironmentOptions( + sim_folder=f"sim_{args.id:02d}", + restart=False, +) + +# Time stepping +time_opts = Time() + +# Geometry +a1 = 1e-2 +a2 = 1.0 +Lz = 4.0 +domain = domains.HollowCylinder(a1=a1, a2=a2, Lz=Lz) + +# Fluid equilibrium (can be used as part of initial conditions) +equil = None + +# Grid +grid = grids.TensorProductGrid(num_elements=(64, 64, 64), mpi_dims_mask=(True, True, True)) + +# Derham options +derham_opts = DerhamOptions(degree=(1, 2, 3), bcs=(("dirichlet", "dirichlet"), None, None)) + +# Profilinig options +profiling_opts = ProfilingOptions( + use_line_profiler=True, +) + +# Simulation object +sim = Simulation( + model=model, + name=name, + description=description, + params_path=__file__, + env=env, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + profiling_opts=profiling_opts, +) + +# ------------------ +# Propagator options +# ------------------ + +from struphy.linear_algebra.solver import SolverParameters + +solver_params = SolverParameters(tol=1e-8, maxiter=3000, info=True, recycle=False) +model.propagators.poisson.options = model.propagators.poisson.Options( + stab_eps=0.0, + solver="pcg", + precond="MultiGrid", + multigrid=MultiGridOptions(), + solver_params=solver_params, +) + +# ------------------ +# Initial conditions +# ------------------ +import numpy as np + +from struphy.initial.base import GenericPerturbation + +# The exact solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5). +# Its value at the Dirichlet boundary eta1 = 0, 1 is exp(-0.25 / w**2) ~ 1e-11, and it is periodic +# in eta2 and eta3 up to the same accuracy. Being localized, it excites the whole spectrum of the +# Laplacian, so CG has to work for convergence. Both the solution and the right-hand side +# -Laplace(phi) are 0-forms given as functions of the logical coordinates (given_in_basis="0"); +# the pushed-forward solution at x = F(eta) is phi(eta). +w = 0.1 +eta0 = (0.5, 0.5, 0.5) + + +def gaussian_derivatives(e1, e2, e3): + """The Gaussian and its first and second derivatives with respect to eta1, eta2, eta3.""" + eta = (e1, e2, e3) + g = np.exp(-sum((e - c) ** 2 for e, c in zip(eta, eta0)) / w**2) + d = [-2 * (e - c) / w**2 * g for e, c in zip(eta, eta0)] + dd = [(4 * (e - c) ** 2 / w**4 - 2 / w**2) * g for e, c in zip(eta, eta0)] + return g, d, dd + + +def exact_solution(e1, e2, e3): + return gaussian_derivatives(e1, e2, e3)[0] + + +def laplacian(e1, e2, e3): + """Physical Laplacian of the exact solution, in logical coordinates. + + HollowCylinder: r = a1 + (a2 - a1) * eta1, theta = 2*pi*eta2, z = Lz * eta3, and + Laplace(phi) = phi_rr + phi_r / r + phi_theta,theta / r**2 + phi_zz. + """ + g, d, dd = gaussian_derivatives(e1, e2, e3) + da = a2 - a1 + r = a1 + da * e1 + return dd[0] / da**2 + d[0] / (da * r) + dd[1] / (2 * np.pi * r) ** 2 + dd[2] / Lz**2 + + +def rhs_fun(e1, e2, e3): + return -laplacian(e1, e2, e3) + + +rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="0") + +model.em_fields.source.add_perturbation(rhs_perturbation) + + +if __name__ == "__main__": + # switches for debugging and testing + estimate_mem = False + run = True + pproc = True + save_figs = True + + if estimate_mem: + sim.estimate_mem(print_report=True) + exit() + + if run: + out = sim.run(profiling_activated=True, one_time_step=True) + else: + # use this for working with existing sim data from a previous run + out = sim.output + + if pproc: + out.pproc(create_vtk=True, parallel=True) + + from matplotlib import pyplot as plt + + def slice_axes(k, X, Y, Z, E1, E2, E3): + """Axes (horizontal, vertical, labels) for plotting the slice eta_k = 0.5: the cross-section (x, y) + at eta3 = 0.5, the half plane (r, z) at eta2 = 0.5 and the logical (eta2, eta3) at eta1 = 0.5.""" + if k == 0: + return E2, E3, "eta2", "eta3" + elif k == 1: + return np.sqrt(X**2 + Y**2), Z, "r", "z" + else: + return X, Y, "x", "y" + + def plot_slices(num, exact, name): + """Numerical solution, exact solution and error on the slices eta_k = 0.5 (k = 1, 2, 3).""" + fig = plt.figure(figsize=(16, 12)) + for k in range(3): + idx = tuple(num.shape[k] // 2 if i == k else slice(None) for i in range(3)) + h, v, hlabel, vlabel = slice_axes(k, X[idx], Y[idx], Z[idx], E1[idx], E2[idx], E3[idx]) + eta_k = (E1, E2, E3)[k][idx].flat[0] + for row, (data, what) in enumerate( + [(num[idx], "from struphy"), (exact[idx], "exact"), (np.abs(num[idx] - exact[idx]), "error")] + ): + plt.subplot(3, 3, 3 * row + k + 1) + plt.pcolormesh(h, v, data, shading="gouraud") + plt.colorbar() + plt.xlabel(hlabel) + plt.ylabel(vlabel) + plt.title(f"{name} {what}, slice at eta{k + 1} = {eta_k:.2f}") + fig.tight_layout() + return fig + + if sim.rank == 0: + # Raw FEEC fields are evaluated directly from the saved spline coefficients at the + # cell centres of the simulation grid (serial Derham, safe on rank 0 only). + rhs = out.evaluate("em_fields/source", t=0).values + + phi_data = out.evaluate("em_fields/phi", t=-1) + phi = phi_data.values + + # logical coordinates of the evaluation points and their images x = F(eta) + E1, E2, E3 = np.meshgrid( + *(phi_data.coords[c].values for c in ("eta1", "eta2", "eta3")), + indexing="ij", + ) + X, Y, Z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) + + # pushed-forward exact solutions: the 0-forms at x = F(eta) are their values at eta + rhs_exact = rhs_fun(E1, E2, E3) + phi_exact = exact_solution(E1, E2, E3) + + fig_rhs = plot_slices(rhs, rhs_exact, "RHS") + fig_phi = plot_slices(phi, phi_exact, "Phi") + + rel_err_rhs = np.max(np.abs(rhs - rhs_exact)) / np.max(np.abs(rhs_exact)) + rel_err_phi = np.max(np.abs(phi - phi_exact)) / np.max(np.abs(phi_exact)) + + print(f"Max relative error in RHS: {rel_err_rhs:.2e}") + print(f"Max relative error in Phi: {rel_err_phi:.2e}") + + assert rel_err_rhs < 2.2e-2, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." + assert rel_err_phi < 1.3e-2, ( + f"The computed solution does not match the exact solution, max rel error = {rel_err_phi}." + ) + + import os + + # `out.path_out` is the out's (absolute) output folder; `sim_folder` alone is a bare + # name resolved against the CWD. The profiling packaging picks these files up from + # here and uploads them as `results-out`. + results_dir = os.path.join(out.path_out, "results") + os.makedirs(results_dir, exist_ok=True) + + np.save(os.path.join(results_dir, "rel_err_rhs.npy"), rel_err_rhs) + np.save(os.path.join(results_dir, "rel_err_phi.npy"), rel_err_phi) + np.save(os.path.join(results_dir, "resolution.npy"), sim.grid.num_elements) + np.save(os.path.join(results_dir, "spline_degree.npy"), sim.derham_opts.degree) + + if save_figs: + fig_rhs.savefig(os.path.join(results_dir, "rhs_slices.png")) + fig_phi.savefig(os.path.join(results_dir, "phi_slices.png")) + else: + plt.show() diff --git a/profiling/examples/Poisson/hollow_cylinder_strong_scaling_stiffness/params_poisson.py b/profiling/examples/Poisson/hollow_cylinder_strong_scaling_stiffness/params_poisson.py new file mode 100644 index 000000000..20d04710c --- /dev/null +++ b/profiling/examples/Poisson/hollow_cylinder_strong_scaling_stiffness/params_poisson.py @@ -0,0 +1,277 @@ +# ----------------------------- +# Description of the simulation +# ----------------------------- +# Please fill in a verbal description of the simulation. +# It will be printed at the beginning of the simulation and can be used to keep track of the different runs. + +name = "Poisson strong scaling on hollow cylinder, stiffness preconditioner" +description = """ +Strong scaling test for Poisson equation on a hollow cylinder (HollowCylinder, a1=1e-2, a2=1, Lz=4). +The manufactured solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5), +exciting the full spectrum. +Homogeneous Dirichlet boundary conditions are set in direction eta1. +The linear system is solved with CG preconditioned by the Kronecker approximation of the inverse Laplacian (StiffnessPreconditioner). +""" + +import logging + +from struphy import set_logging_level + +set_logging_level(logging.INFO) + +import argparse + +# ------------------ +# Import Struphy API +# ------------------ +from struphy import ( + BaseUnits, + DerhamOptions, + EnvironmentOptions, + FieldsBackground, + ProfilingOptions, + Simulation, + Time, + domains, + equils, + grids, + perturbations, +) + +# --------------------- +# Instance of the model +# --------------------- +from struphy.models import Poisson + +# Units +base_units = BaseUnits() + +# Model instance +model = Poisson(base_units=base_units) + +# List all variables and decide whether to save their data +model.em_fields.phi.save_data = True +model.em_fields.source.save_data = True + +# -------------------------- +# Instance of the simulation +# -------------------------- + +# Environment options +# `--id` distinguishes runs that share a rank count but differ in something else; the +# profiling driver passes its launch counter (see `ProfilingJob.build_commands`). +# Unknown flags are ignored so the driver can forward other parameters as well. +parser = argparse.ArgumentParser() +parser.add_argument("--id", type=int, default=0, help="Run id, used to name the output folder.") +args, _ = parser.parse_known_args() + +env = EnvironmentOptions( + sim_folder=f"sim_{args.id:02d}", + restart=False, +) + +# Time stepping +time_opts = Time() + +# Geometry +a1 = 1e-2 +a2 = 1.0 +Lz = 4.0 +domain = domains.HollowCylinder(a1=a1, a2=a2, Lz=Lz) + +# Fluid equilibrium (can be used as part of initial conditions) +equil = None + +# Grid +grid = grids.TensorProductGrid(num_elements=(64, 64, 64), mpi_dims_mask=(True, True, True)) + +# Derham options +derham_opts = DerhamOptions(degree=(1, 2, 3), bcs=(("dirichlet", "dirichlet"), None, None)) + +# Profilinig options +profiling_opts = ProfilingOptions( + use_line_profiler=True, +) + +# Simulation object +sim = Simulation( + model=model, + name=name, + description=description, + params_path=__file__, + env=env, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + profiling_opts=profiling_opts, +) + +# ------------------ +# Propagator options +# ------------------ + +from struphy.linear_algebra.solver import SolverParameters + +solver_params = SolverParameters(tol=1e-8, maxiter=3000, info=True, recycle=False) +model.propagators.poisson.options = model.propagators.poisson.Options( + stab_eps=0.0, + solver="pcg", + precond="StiffnessPreconditioner", + solver_params=solver_params, +) + +# ------------------ +# Initial conditions +# ------------------ +import numpy as np + +from struphy.initial.base import GenericPerturbation + +# The exact solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5). +# Its value at the Dirichlet boundary eta1 = 0, 1 is exp(-0.25 / w**2) ~ 1e-11, and it is periodic +# in eta2 and eta3 up to the same accuracy. Being localized, it excites the whole spectrum of the +# Laplacian, so CG has to work for convergence. Both the solution and the right-hand side +# -Laplace(phi) are 0-forms given as functions of the logical coordinates (given_in_basis="0"); +# the pushed-forward solution at x = F(eta) is phi(eta). +w = 0.1 +eta0 = (0.5, 0.5, 0.5) + + +def gaussian_derivatives(e1, e2, e3): + """The Gaussian and its first and second derivatives with respect to eta1, eta2, eta3.""" + eta = (e1, e2, e3) + g = np.exp(-sum((e - c) ** 2 for e, c in zip(eta, eta0)) / w**2) + d = [-2 * (e - c) / w**2 * g for e, c in zip(eta, eta0)] + dd = [(4 * (e - c) ** 2 / w**4 - 2 / w**2) * g for e, c in zip(eta, eta0)] + return g, d, dd + + +def exact_solution(e1, e2, e3): + return gaussian_derivatives(e1, e2, e3)[0] + + +def laplacian(e1, e2, e3): + """Physical Laplacian of the exact solution, in logical coordinates. + + HollowCylinder: r = a1 + (a2 - a1) * eta1, theta = 2*pi*eta2, z = Lz * eta3, and + Laplace(phi) = phi_rr + phi_r / r + phi_theta,theta / r**2 + phi_zz. + """ + g, d, dd = gaussian_derivatives(e1, e2, e3) + da = a2 - a1 + r = a1 + da * e1 + return dd[0] / da**2 + d[0] / (da * r) + dd[1] / (2 * np.pi * r) ** 2 + dd[2] / Lz**2 + + +def rhs_fun(e1, e2, e3): + return -laplacian(e1, e2, e3) + + +rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="0") + +model.em_fields.source.add_perturbation(rhs_perturbation) + + +if __name__ == "__main__": + # switches for debugging and testing + estimate_mem = False + run = True + pproc = True + save_figs = True + + if estimate_mem: + sim.estimate_mem(print_report=True) + exit() + + if run: + out = sim.run(profiling_activated=True, one_time_step=True) + else: + # use this for working with existing sim data from a previous run + out = sim.output + + if pproc: + out.pproc(create_vtk=True, parallel=True) + + from matplotlib import pyplot as plt + + def slice_axes(k, X, Y, Z, E1, E2, E3): + """Axes (horizontal, vertical, labels) for plotting the slice eta_k = 0.5: the cross-section (x, y) + at eta3 = 0.5, the half plane (r, z) at eta2 = 0.5 and the logical (eta2, eta3) at eta1 = 0.5.""" + if k == 0: + return E2, E3, "eta2", "eta3" + elif k == 1: + return np.sqrt(X**2 + Y**2), Z, "r", "z" + else: + return X, Y, "x", "y" + + def plot_slices(num, exact, name): + """Numerical solution, exact solution and error on the slices eta_k = 0.5 (k = 1, 2, 3).""" + fig = plt.figure(figsize=(16, 12)) + for k in range(3): + idx = tuple(num.shape[k] // 2 if i == k else slice(None) for i in range(3)) + h, v, hlabel, vlabel = slice_axes(k, X[idx], Y[idx], Z[idx], E1[idx], E2[idx], E3[idx]) + eta_k = (E1, E2, E3)[k][idx].flat[0] + for row, (data, what) in enumerate( + [(num[idx], "from struphy"), (exact[idx], "exact"), (np.abs(num[idx] - exact[idx]), "error")] + ): + plt.subplot(3, 3, 3 * row + k + 1) + plt.pcolormesh(h, v, data, shading="gouraud") + plt.colorbar() + plt.xlabel(hlabel) + plt.ylabel(vlabel) + plt.title(f"{name} {what}, slice at eta{k + 1} = {eta_k:.2f}") + fig.tight_layout() + return fig + + if sim.rank == 0: + # Raw FEEC fields are evaluated directly from the saved spline coefficients at the + # cell centres of the simulation grid (serial Derham, safe on rank 0 only). + rhs = out.evaluate("em_fields/source", t=0).values + + phi_data = out.evaluate("em_fields/phi", t=-1) + phi = phi_data.values + + # logical coordinates of the evaluation points and their images x = F(eta) + E1, E2, E3 = np.meshgrid( + *(phi_data.coords[c].values for c in ("eta1", "eta2", "eta3")), + indexing="ij", + ) + X, Y, Z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) + + # pushed-forward exact solutions: the 0-forms at x = F(eta) are their values at eta + rhs_exact = rhs_fun(E1, E2, E3) + phi_exact = exact_solution(E1, E2, E3) + + fig_rhs = plot_slices(rhs, rhs_exact, "RHS") + fig_phi = plot_slices(phi, phi_exact, "Phi") + + rel_err_rhs = np.max(np.abs(rhs - rhs_exact)) / np.max(np.abs(rhs_exact)) + rel_err_phi = np.max(np.abs(phi - phi_exact)) / np.max(np.abs(phi_exact)) + + print(f"Max relative error in RHS: {rel_err_rhs:.2e}") + print(f"Max relative error in Phi: {rel_err_phi:.2e}") + + assert rel_err_rhs < 2.2e-2, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." + assert rel_err_phi < 1.3e-2, ( + f"The computed solution does not match the exact solution, max rel error = {rel_err_phi}." + ) + + import os + + # `out.path_out` is the out's (absolute) output folder; `sim_folder` alone is a bare + # name resolved against the CWD. The profiling packaging picks these files up from + # here and uploads them as `results-out`. + results_dir = os.path.join(out.path_out, "results") + os.makedirs(results_dir, exist_ok=True) + + np.save(os.path.join(results_dir, "rel_err_rhs.npy"), rel_err_rhs) + np.save(os.path.join(results_dir, "rel_err_phi.npy"), rel_err_phi) + np.save(os.path.join(results_dir, "resolution.npy"), sim.grid.num_elements) + np.save(os.path.join(results_dir, "spline_degree.npy"), sim.derham_opts.degree) + + if save_figs: + fig_rhs.savefig(os.path.join(results_dir, "rhs_slices.png")) + fig_phi.savefig(os.path.join(results_dir, "phi_slices.png")) + else: + plt.show() diff --git a/profiling/examples/Poisson/hollow_torus_strong_scaling/params_poisson.py b/profiling/examples/Poisson/hollow_torus_strong_scaling/params_poisson.py new file mode 100644 index 000000000..5a8584d6f --- /dev/null +++ b/profiling/examples/Poisson/hollow_torus_strong_scaling/params_poisson.py @@ -0,0 +1,289 @@ +# ----------------------------- +# Description of the simulation +# ----------------------------- +# Please fill in a verbal description of the simulation. +# It will be printed at the beginning of the simulation and can be used to keep track of the different runs. + +name = "Poisson strong scaling on hollow torus, no preconditioner" +description = """ +Strong scaling test for Poisson equation on a full hollow torus (HollowTorus, a1=1e-2, a2=1, R0=4). +The manufactured solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5), +exciting the full spectrum. +Homogeneous Dirichlet boundary conditions are set in direction eta1. +The linear system is solved with unpreconditioned CG. +""" + +import logging + +from struphy import set_logging_level + +set_logging_level(logging.INFO) + +import argparse + +# ------------------ +# Import Struphy API +# ------------------ +from struphy import ( + BaseUnits, + DerhamOptions, + EnvironmentOptions, + FieldsBackground, + ProfilingOptions, + Simulation, + Time, + domains, + equils, + grids, + perturbations, +) + +# --------------------- +# Instance of the model +# --------------------- +from struphy.models import Poisson + +# Units +base_units = BaseUnits() + +# Model instance +model = Poisson(base_units=base_units) + +# List all variables and decide whether to save their data +model.em_fields.phi.save_data = True +model.em_fields.source.save_data = True + +# -------------------------- +# Instance of the simulation +# -------------------------- + +# Environment options +# `--id` distinguishes runs that share a rank count but differ in something else; the +# profiling driver passes its launch counter (see `ProfilingJob.build_commands`). +# Unknown flags are ignored so the driver can forward other parameters as well. +parser = argparse.ArgumentParser() +parser.add_argument("--id", type=int, default=0, help="Run id, used to name the output folder.") +args, _ = parser.parse_known_args() + +env = EnvironmentOptions( + sim_folder=f"sim_{args.id:02d}", + restart=False, +) + +# Time stepping +time_opts = Time() + +# Geometry +a1 = 1e-2 +a2 = 1.0 +R0 = 4.0 +tor_period = 1 +domain = domains.HollowTorus(a1=a1, a2=a2, R0=R0, sfl=False, pol_period=1, tor_period=tor_period) + +# Fluid equilibrium (can be used as part of initial conditions) +equil = None + +# Grid +grid = grids.TensorProductGrid(num_elements=(64, 64, 64), mpi_dims_mask=(True, True, True)) + +# Derham options +derham_opts = DerhamOptions(degree=(1, 2, 3), bcs=(("dirichlet", "dirichlet"), None, None)) + +# Profilinig options +profiling_opts = ProfilingOptions( + use_line_profiler=True, +) + +# Simulation object +sim = Simulation( + model=model, + name=name, + description=description, + params_path=__file__, + env=env, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + profiling_opts=profiling_opts, +) + +# ------------------ +# Propagator options +# ------------------ + +from struphy.linear_algebra.solver import SolverParameters + +solver_params = SolverParameters(tol=1e-8, maxiter=3000, info=True, recycle=False) +model.propagators.poisson.options = model.propagators.poisson.Options( + stab_eps=0.0, + solver="pcg", + precond=None, + solver_params=solver_params, +) + +# ------------------ +# Initial conditions +# ------------------ +import numpy as np + +from struphy.initial.base import GenericPerturbation + +# The exact solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5). +# Its value at the Dirichlet boundary eta1 = 0, 1 is exp(-0.25 / w**2) ~ 1e-11, and it is periodic +# in eta2 and eta3 up to the same accuracy. Being localized, it excites the whole spectrum of the +# Laplacian, so CG has to work for convergence. Both the solution and the right-hand side +# -Laplace(phi) are 0-forms given as functions of the logical coordinates (given_in_basis="0"); +# the pushed-forward solution at x = F(eta) is phi(eta). +w = 0.1 +eta0 = (0.5, 0.5, 0.5) + + +def gaussian_derivatives(e1, e2, e3): + """The Gaussian and its first and second derivatives with respect to eta1, eta2, eta3.""" + eta = (e1, e2, e3) + g = np.exp(-sum((e - c) ** 2 for e, c in zip(eta, eta0)) / w**2) + d = [-2 * (e - c) / w**2 * g for e, c in zip(eta, eta0)] + dd = [(4 * (e - c) ** 2 / w**4 - 2 / w**2) * g for e, c in zip(eta, eta0)] + return g, d, dd + + +def exact_solution(e1, e2, e3): + return gaussian_derivatives(e1, e2, e3)[0] + + +def laplacian(e1, e2, e3): + """Physical Laplacian of the exact solution, in logical coordinates. + + HollowTorus (sfl=False): r = a1 + (a2 - a1) * eta1, theta = 2*pi*eta2, phi = 2*pi*eta3 / tor_period, + R = R0 + r*cos(theta), metric dr**2 + r**2 dtheta**2 + R**2 dphi**2, and + Laplace(u) = u_rr + (R + r*cos(theta)) / (r*R) * u_r + u_theta,theta / r**2 - sin(theta) / (r*R) * u_theta + + u_phi,phi / R**2. + """ + g, d, dd = gaussian_derivatives(e1, e2, e3) + da = a2 - a1 + r = a1 + da * e1 + theta = 2 * np.pi * e2 + R = R0 + r * np.cos(theta) + dphi = 2 * np.pi / tor_period + return ( + dd[0] / da**2 + + (R + r * np.cos(theta)) / (r * R) * d[0] / da + + dd[1] / (2 * np.pi * r) ** 2 + - np.sin(theta) / (r * R) * d[1] / (2 * np.pi) + + dd[2] / (dphi * R) ** 2 + ) + + +def rhs_fun(e1, e2, e3): + return -laplacian(e1, e2, e3) + + +rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="0") + +model.em_fields.source.add_perturbation(rhs_perturbation) + + +if __name__ == "__main__": + # switches for debugging and testing + estimate_mem = False + run = True + pproc = True + save_figs = True + + if estimate_mem: + sim.estimate_mem(print_report=True) + exit() + + if run: + out = sim.run(profiling_activated=True, one_time_step=True) + else: + # use this for working with existing sim data from a previous run + out = sim.output + + if pproc: + out.pproc(create_vtk=True, parallel=True) + + from matplotlib import pyplot as plt + + def slice_axes(k, X, Y, Z, E1, E2, E3): + """Axes (horizontal, vertical, labels) for plotting the slice eta_k = 0.5: the poloidal plane (R, z) + at eta3 = 0.5 and the logical planes at eta1 = 0.5 and eta2 = 0.5.""" + if k == 0: + return E2, E3, "eta2", "eta3" + elif k == 1: + return E1, E3, "eta1", "eta3" + else: + return np.sqrt(X**2 + Y**2), Z, "R", "z" + + def plot_slices(num, exact, name): + """Numerical solution, exact solution and error on the slices eta_k = 0.5 (k = 1, 2, 3).""" + fig = plt.figure(figsize=(16, 12)) + for k in range(3): + idx = tuple(num.shape[k] // 2 if i == k else slice(None) for i in range(3)) + h, v, hlabel, vlabel = slice_axes(k, X[idx], Y[idx], Z[idx], E1[idx], E2[idx], E3[idx]) + eta_k = (E1, E2, E3)[k][idx].flat[0] + for row, (data, what) in enumerate( + [(num[idx], "from struphy"), (exact[idx], "exact"), (np.abs(num[idx] - exact[idx]), "error")] + ): + plt.subplot(3, 3, 3 * row + k + 1) + plt.pcolormesh(h, v, data, shading="gouraud") + plt.colorbar() + plt.xlabel(hlabel) + plt.ylabel(vlabel) + plt.title(f"{name} {what}, slice at eta{k + 1} = {eta_k:.2f}") + fig.tight_layout() + return fig + + if sim.rank == 0: + # Raw FEEC fields are evaluated directly from the saved spline coefficients at the + # cell centres of the simulation grid (serial Derham, safe on rank 0 only). + rhs = out.evaluate("em_fields/source", t=0).values + + phi_data = out.evaluate("em_fields/phi", t=-1) + phi = phi_data.values + + # logical coordinates of the evaluation points and their images x = F(eta) + E1, E2, E3 = np.meshgrid( + *(phi_data.coords[c].values for c in ("eta1", "eta2", "eta3")), + indexing="ij", + ) + X, Y, Z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) + + # pushed-forward exact solutions: the 0-forms at x = F(eta) are their values at eta + rhs_exact = rhs_fun(E1, E2, E3) + phi_exact = exact_solution(E1, E2, E3) + + fig_rhs = plot_slices(rhs, rhs_exact, "RHS") + fig_phi = plot_slices(phi, phi_exact, "Phi") + + rel_err_rhs = np.max(np.abs(rhs - rhs_exact)) / np.max(np.abs(rhs_exact)) + rel_err_phi = np.max(np.abs(phi - phi_exact)) / np.max(np.abs(phi_exact)) + + print(f"Max relative error in RHS: {rel_err_rhs:.2e}") + print(f"Max relative error in Phi: {rel_err_phi:.2e}") + + assert rel_err_rhs < 2.2e-2, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." + assert rel_err_phi < 1.3e-2, ( + f"The computed solution does not match the exact solution, max rel error = {rel_err_phi}." + ) + + import os + + # `out.path_out` is the out's (absolute) output folder; `sim_folder` alone is a bare + # name resolved against the CWD. The profiling packaging picks these files up from + # here and uploads them as `results-out`. + results_dir = os.path.join(out.path_out, "results") + os.makedirs(results_dir, exist_ok=True) + + np.save(os.path.join(results_dir, "rel_err_rhs.npy"), rel_err_rhs) + np.save(os.path.join(results_dir, "rel_err_phi.npy"), rel_err_phi) + np.save(os.path.join(results_dir, "resolution.npy"), sim.grid.num_elements) + np.save(os.path.join(results_dir, "spline_degree.npy"), sim.derham_opts.degree) + + if save_figs: + fig_rhs.savefig(os.path.join(results_dir, "rhs_slices.png")) + fig_phi.savefig(os.path.join(results_dir, "phi_slices.png")) + else: + plt.show() diff --git a/profiling/examples/Poisson/hollow_torus_strong_scaling_massmatrix/params_poisson.py b/profiling/examples/Poisson/hollow_torus_strong_scaling_massmatrix/params_poisson.py new file mode 100644 index 000000000..67552ef1b --- /dev/null +++ b/profiling/examples/Poisson/hollow_torus_strong_scaling_massmatrix/params_poisson.py @@ -0,0 +1,289 @@ +# ----------------------------- +# Description of the simulation +# ----------------------------- +# Please fill in a verbal description of the simulation. +# It will be printed at the beginning of the simulation and can be used to keep track of the different runs. + +name = "Poisson strong scaling on hollow torus, mass-matrix preconditioner" +description = """ +Strong scaling test for Poisson equation on a full hollow torus (HollowTorus, a1=1e-2, a2=1, R0=4). +The manufactured solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5), +exciting the full spectrum. +Homogeneous Dirichlet boundary conditions are set in direction eta1. +The linear system is solved with CG preconditioned by the Kronecker approximation of the inverse 0-form mass matrix (MassMatrixPreconditioner). +""" + +import logging + +from struphy import set_logging_level + +set_logging_level(logging.INFO) + +import argparse + +# ------------------ +# Import Struphy API +# ------------------ +from struphy import ( + BaseUnits, + DerhamOptions, + EnvironmentOptions, + FieldsBackground, + ProfilingOptions, + Simulation, + Time, + domains, + equils, + grids, + perturbations, +) + +# --------------------- +# Instance of the model +# --------------------- +from struphy.models import Poisson + +# Units +base_units = BaseUnits() + +# Model instance +model = Poisson(base_units=base_units) + +# List all variables and decide whether to save their data +model.em_fields.phi.save_data = True +model.em_fields.source.save_data = True + +# -------------------------- +# Instance of the simulation +# -------------------------- + +# Environment options +# `--id` distinguishes runs that share a rank count but differ in something else; the +# profiling driver passes its launch counter (see `ProfilingJob.build_commands`). +# Unknown flags are ignored so the driver can forward other parameters as well. +parser = argparse.ArgumentParser() +parser.add_argument("--id", type=int, default=0, help="Run id, used to name the output folder.") +args, _ = parser.parse_known_args() + +env = EnvironmentOptions( + sim_folder=f"sim_{args.id:02d}", + restart=False, +) + +# Time stepping +time_opts = Time() + +# Geometry +a1 = 1e-2 +a2 = 1.0 +R0 = 4.0 +tor_period = 1 +domain = domains.HollowTorus(a1=a1, a2=a2, R0=R0, sfl=False, pol_period=1, tor_period=tor_period) + +# Fluid equilibrium (can be used as part of initial conditions) +equil = None + +# Grid +grid = grids.TensorProductGrid(num_elements=(64, 64, 64), mpi_dims_mask=(True, True, True)) + +# Derham options +derham_opts = DerhamOptions(degree=(1, 2, 3), bcs=(("dirichlet", "dirichlet"), None, None)) + +# Profilinig options +profiling_opts = ProfilingOptions( + use_line_profiler=True, +) + +# Simulation object +sim = Simulation( + model=model, + name=name, + description=description, + params_path=__file__, + env=env, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + profiling_opts=profiling_opts, +) + +# ------------------ +# Propagator options +# ------------------ + +from struphy.linear_algebra.solver import SolverParameters + +solver_params = SolverParameters(tol=1e-8, maxiter=3000, info=True, recycle=False) +model.propagators.poisson.options = model.propagators.poisson.Options( + stab_eps=0.0, + solver="pcg", + precond="MassMatrixPreconditioner", + solver_params=solver_params, +) + +# ------------------ +# Initial conditions +# ------------------ +import numpy as np + +from struphy.initial.base import GenericPerturbation + +# The exact solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5). +# Its value at the Dirichlet boundary eta1 = 0, 1 is exp(-0.25 / w**2) ~ 1e-11, and it is periodic +# in eta2 and eta3 up to the same accuracy. Being localized, it excites the whole spectrum of the +# Laplacian, so CG has to work for convergence. Both the solution and the right-hand side +# -Laplace(phi) are 0-forms given as functions of the logical coordinates (given_in_basis="0"); +# the pushed-forward solution at x = F(eta) is phi(eta). +w = 0.1 +eta0 = (0.5, 0.5, 0.5) + + +def gaussian_derivatives(e1, e2, e3): + """The Gaussian and its first and second derivatives with respect to eta1, eta2, eta3.""" + eta = (e1, e2, e3) + g = np.exp(-sum((e - c) ** 2 for e, c in zip(eta, eta0)) / w**2) + d = [-2 * (e - c) / w**2 * g for e, c in zip(eta, eta0)] + dd = [(4 * (e - c) ** 2 / w**4 - 2 / w**2) * g for e, c in zip(eta, eta0)] + return g, d, dd + + +def exact_solution(e1, e2, e3): + return gaussian_derivatives(e1, e2, e3)[0] + + +def laplacian(e1, e2, e3): + """Physical Laplacian of the exact solution, in logical coordinates. + + HollowTorus (sfl=False): r = a1 + (a2 - a1) * eta1, theta = 2*pi*eta2, phi = 2*pi*eta3 / tor_period, + R = R0 + r*cos(theta), metric dr**2 + r**2 dtheta**2 + R**2 dphi**2, and + Laplace(u) = u_rr + (R + r*cos(theta)) / (r*R) * u_r + u_theta,theta / r**2 - sin(theta) / (r*R) * u_theta + + u_phi,phi / R**2. + """ + g, d, dd = gaussian_derivatives(e1, e2, e3) + da = a2 - a1 + r = a1 + da * e1 + theta = 2 * np.pi * e2 + R = R0 + r * np.cos(theta) + dphi = 2 * np.pi / tor_period + return ( + dd[0] / da**2 + + (R + r * np.cos(theta)) / (r * R) * d[0] / da + + dd[1] / (2 * np.pi * r) ** 2 + - np.sin(theta) / (r * R) * d[1] / (2 * np.pi) + + dd[2] / (dphi * R) ** 2 + ) + + +def rhs_fun(e1, e2, e3): + return -laplacian(e1, e2, e3) + + +rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="0") + +model.em_fields.source.add_perturbation(rhs_perturbation) + + +if __name__ == "__main__": + # switches for debugging and testing + estimate_mem = False + run = True + pproc = True + save_figs = True + + if estimate_mem: + sim.estimate_mem(print_report=True) + exit() + + if run: + out = sim.run(profiling_activated=True, one_time_step=True) + else: + # use this for working with existing sim data from a previous run + out = sim.output + + if pproc: + out.pproc(create_vtk=True, parallel=True) + + from matplotlib import pyplot as plt + + def slice_axes(k, X, Y, Z, E1, E2, E3): + """Axes (horizontal, vertical, labels) for plotting the slice eta_k = 0.5: the poloidal plane (R, z) + at eta3 = 0.5 and the logical planes at eta1 = 0.5 and eta2 = 0.5.""" + if k == 0: + return E2, E3, "eta2", "eta3" + elif k == 1: + return E1, E3, "eta1", "eta3" + else: + return np.sqrt(X**2 + Y**2), Z, "R", "z" + + def plot_slices(num, exact, name): + """Numerical solution, exact solution and error on the slices eta_k = 0.5 (k = 1, 2, 3).""" + fig = plt.figure(figsize=(16, 12)) + for k in range(3): + idx = tuple(num.shape[k] // 2 if i == k else slice(None) for i in range(3)) + h, v, hlabel, vlabel = slice_axes(k, X[idx], Y[idx], Z[idx], E1[idx], E2[idx], E3[idx]) + eta_k = (E1, E2, E3)[k][idx].flat[0] + for row, (data, what) in enumerate( + [(num[idx], "from struphy"), (exact[idx], "exact"), (np.abs(num[idx] - exact[idx]), "error")] + ): + plt.subplot(3, 3, 3 * row + k + 1) + plt.pcolormesh(h, v, data, shading="gouraud") + plt.colorbar() + plt.xlabel(hlabel) + plt.ylabel(vlabel) + plt.title(f"{name} {what}, slice at eta{k + 1} = {eta_k:.2f}") + fig.tight_layout() + return fig + + if sim.rank == 0: + # Raw FEEC fields are evaluated directly from the saved spline coefficients at the + # cell centres of the simulation grid (serial Derham, safe on rank 0 only). + rhs = out.evaluate("em_fields/source", t=0).values + + phi_data = out.evaluate("em_fields/phi", t=-1) + phi = phi_data.values + + # logical coordinates of the evaluation points and their images x = F(eta) + E1, E2, E3 = np.meshgrid( + *(phi_data.coords[c].values for c in ("eta1", "eta2", "eta3")), + indexing="ij", + ) + X, Y, Z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) + + # pushed-forward exact solutions: the 0-forms at x = F(eta) are their values at eta + rhs_exact = rhs_fun(E1, E2, E3) + phi_exact = exact_solution(E1, E2, E3) + + fig_rhs = plot_slices(rhs, rhs_exact, "RHS") + fig_phi = plot_slices(phi, phi_exact, "Phi") + + rel_err_rhs = np.max(np.abs(rhs - rhs_exact)) / np.max(np.abs(rhs_exact)) + rel_err_phi = np.max(np.abs(phi - phi_exact)) / np.max(np.abs(phi_exact)) + + print(f"Max relative error in RHS: {rel_err_rhs:.2e}") + print(f"Max relative error in Phi: {rel_err_phi:.2e}") + + assert rel_err_rhs < 2.2e-2, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." + assert rel_err_phi < 1.3e-2, ( + f"The computed solution does not match the exact solution, max rel error = {rel_err_phi}." + ) + + import os + + # `out.path_out` is the out's (absolute) output folder; `sim_folder` alone is a bare + # name resolved against the CWD. The profiling packaging picks these files up from + # here and uploads them as `results-out`. + results_dir = os.path.join(out.path_out, "results") + os.makedirs(results_dir, exist_ok=True) + + np.save(os.path.join(results_dir, "rel_err_rhs.npy"), rel_err_rhs) + np.save(os.path.join(results_dir, "rel_err_phi.npy"), rel_err_phi) + np.save(os.path.join(results_dir, "resolution.npy"), sim.grid.num_elements) + np.save(os.path.join(results_dir, "spline_degree.npy"), sim.derham_opts.degree) + + if save_figs: + fig_rhs.savefig(os.path.join(results_dir, "rhs_slices.png")) + fig_phi.savefig(os.path.join(results_dir, "phi_slices.png")) + else: + plt.show() diff --git a/profiling/examples/Poisson/hollow_torus_strong_scaling_multigrid/params_poisson.py b/profiling/examples/Poisson/hollow_torus_strong_scaling_multigrid/params_poisson.py new file mode 100644 index 000000000..54c8f42d3 --- /dev/null +++ b/profiling/examples/Poisson/hollow_torus_strong_scaling_multigrid/params_poisson.py @@ -0,0 +1,291 @@ +# ----------------------------- +# Description of the simulation +# ----------------------------- +# Please fill in a verbal description of the simulation. +# It will be printed at the beginning of the simulation and can be used to keep track of the different runs. + +name = "Poisson strong scaling on hollow torus, multigrid preconditioner" +description = """ +Strong scaling test for Poisson equation on a full hollow torus (HollowTorus, a1=1e-2, a2=1, R0=4). +The manufactured solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5), +exciting the full spectrum. +Homogeneous Dirichlet boundary conditions are set in direction eta1. +The linear system is solved with CG preconditioned by geometric multigrid. +""" + +import logging + +from struphy import set_logging_level + +set_logging_level(logging.INFO) + +import argparse + +# ------------------ +# Import Struphy API +# ------------------ +from struphy import ( + BaseUnits, + DerhamOptions, + EnvironmentOptions, + FieldsBackground, + MultiGridOptions, + ProfilingOptions, + Simulation, + Time, + domains, + equils, + grids, + perturbations, +) + +# --------------------- +# Instance of the model +# --------------------- +from struphy.models import Poisson + +# Units +base_units = BaseUnits() + +# Model instance +model = Poisson(base_units=base_units) + +# List all variables and decide whether to save their data +model.em_fields.phi.save_data = True +model.em_fields.source.save_data = True + +# -------------------------- +# Instance of the simulation +# -------------------------- + +# Environment options +# `--id` distinguishes runs that share a rank count but differ in something else; the +# profiling driver passes its launch counter (see `ProfilingJob.build_commands`). +# Unknown flags are ignored so the driver can forward other parameters as well. +parser = argparse.ArgumentParser() +parser.add_argument("--id", type=int, default=0, help="Run id, used to name the output folder.") +args, _ = parser.parse_known_args() + +env = EnvironmentOptions( + sim_folder=f"sim_{args.id:02d}", + restart=False, +) + +# Time stepping +time_opts = Time() + +# Geometry +a1 = 1e-2 +a2 = 1.0 +R0 = 4.0 +tor_period = 1 +domain = domains.HollowTorus(a1=a1, a2=a2, R0=R0, sfl=False, pol_period=1, tor_period=tor_period) + +# Fluid equilibrium (can be used as part of initial conditions) +equil = None + +# Grid +grid = grids.TensorProductGrid(num_elements=(64, 64, 64), mpi_dims_mask=(True, True, True)) + +# Derham options +derham_opts = DerhamOptions(degree=(1, 2, 3), bcs=(("dirichlet", "dirichlet"), None, None)) + +# Profilinig options +profiling_opts = ProfilingOptions( + use_line_profiler=True, +) + +# Simulation object +sim = Simulation( + model=model, + name=name, + description=description, + params_path=__file__, + env=env, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + profiling_opts=profiling_opts, +) + +# ------------------ +# Propagator options +# ------------------ + +from struphy.linear_algebra.solver import SolverParameters + +solver_params = SolverParameters(tol=1e-8, maxiter=3000, info=True, recycle=False) +model.propagators.poisson.options = model.propagators.poisson.Options( + stab_eps=0.0, + solver="pcg", + precond="MultiGrid", + multigrid=MultiGridOptions(), + solver_params=solver_params, +) + +# ------------------ +# Initial conditions +# ------------------ +import numpy as np + +from struphy.initial.base import GenericPerturbation + +# The exact solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5). +# Its value at the Dirichlet boundary eta1 = 0, 1 is exp(-0.25 / w**2) ~ 1e-11, and it is periodic +# in eta2 and eta3 up to the same accuracy. Being localized, it excites the whole spectrum of the +# Laplacian, so CG has to work for convergence. Both the solution and the right-hand side +# -Laplace(phi) are 0-forms given as functions of the logical coordinates (given_in_basis="0"); +# the pushed-forward solution at x = F(eta) is phi(eta). +w = 0.1 +eta0 = (0.5, 0.5, 0.5) + + +def gaussian_derivatives(e1, e2, e3): + """The Gaussian and its first and second derivatives with respect to eta1, eta2, eta3.""" + eta = (e1, e2, e3) + g = np.exp(-sum((e - c) ** 2 for e, c in zip(eta, eta0)) / w**2) + d = [-2 * (e - c) / w**2 * g for e, c in zip(eta, eta0)] + dd = [(4 * (e - c) ** 2 / w**4 - 2 / w**2) * g for e, c in zip(eta, eta0)] + return g, d, dd + + +def exact_solution(e1, e2, e3): + return gaussian_derivatives(e1, e2, e3)[0] + + +def laplacian(e1, e2, e3): + """Physical Laplacian of the exact solution, in logical coordinates. + + HollowTorus (sfl=False): r = a1 + (a2 - a1) * eta1, theta = 2*pi*eta2, phi = 2*pi*eta3 / tor_period, + R = R0 + r*cos(theta), metric dr**2 + r**2 dtheta**2 + R**2 dphi**2, and + Laplace(u) = u_rr + (R + r*cos(theta)) / (r*R) * u_r + u_theta,theta / r**2 - sin(theta) / (r*R) * u_theta + + u_phi,phi / R**2. + """ + g, d, dd = gaussian_derivatives(e1, e2, e3) + da = a2 - a1 + r = a1 + da * e1 + theta = 2 * np.pi * e2 + R = R0 + r * np.cos(theta) + dphi = 2 * np.pi / tor_period + return ( + dd[0] / da**2 + + (R + r * np.cos(theta)) / (r * R) * d[0] / da + + dd[1] / (2 * np.pi * r) ** 2 + - np.sin(theta) / (r * R) * d[1] / (2 * np.pi) + + dd[2] / (dphi * R) ** 2 + ) + + +def rhs_fun(e1, e2, e3): + return -laplacian(e1, e2, e3) + + +rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="0") + +model.em_fields.source.add_perturbation(rhs_perturbation) + + +if __name__ == "__main__": + # switches for debugging and testing + estimate_mem = False + run = True + pproc = True + save_figs = True + + if estimate_mem: + sim.estimate_mem(print_report=True) + exit() + + if run: + out = sim.run(profiling_activated=True, one_time_step=True) + else: + # use this for working with existing sim data from a previous run + out = sim.output + + if pproc: + out.pproc(create_vtk=True, parallel=True) + + from matplotlib import pyplot as plt + + def slice_axes(k, X, Y, Z, E1, E2, E3): + """Axes (horizontal, vertical, labels) for plotting the slice eta_k = 0.5: the poloidal plane (R, z) + at eta3 = 0.5 and the logical planes at eta1 = 0.5 and eta2 = 0.5.""" + if k == 0: + return E2, E3, "eta2", "eta3" + elif k == 1: + return E1, E3, "eta1", "eta3" + else: + return np.sqrt(X**2 + Y**2), Z, "R", "z" + + def plot_slices(num, exact, name): + """Numerical solution, exact solution and error on the slices eta_k = 0.5 (k = 1, 2, 3).""" + fig = plt.figure(figsize=(16, 12)) + for k in range(3): + idx = tuple(num.shape[k] // 2 if i == k else slice(None) for i in range(3)) + h, v, hlabel, vlabel = slice_axes(k, X[idx], Y[idx], Z[idx], E1[idx], E2[idx], E3[idx]) + eta_k = (E1, E2, E3)[k][idx].flat[0] + for row, (data, what) in enumerate( + [(num[idx], "from struphy"), (exact[idx], "exact"), (np.abs(num[idx] - exact[idx]), "error")] + ): + plt.subplot(3, 3, 3 * row + k + 1) + plt.pcolormesh(h, v, data, shading="gouraud") + plt.colorbar() + plt.xlabel(hlabel) + plt.ylabel(vlabel) + plt.title(f"{name} {what}, slice at eta{k + 1} = {eta_k:.2f}") + fig.tight_layout() + return fig + + if sim.rank == 0: + # Raw FEEC fields are evaluated directly from the saved spline coefficients at the + # cell centres of the simulation grid (serial Derham, safe on rank 0 only). + rhs = out.evaluate("em_fields/source", t=0).values + + phi_data = out.evaluate("em_fields/phi", t=-1) + phi = phi_data.values + + # logical coordinates of the evaluation points and their images x = F(eta) + E1, E2, E3 = np.meshgrid( + *(phi_data.coords[c].values for c in ("eta1", "eta2", "eta3")), + indexing="ij", + ) + X, Y, Z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) + + # pushed-forward exact solutions: the 0-forms at x = F(eta) are their values at eta + rhs_exact = rhs_fun(E1, E2, E3) + phi_exact = exact_solution(E1, E2, E3) + + fig_rhs = plot_slices(rhs, rhs_exact, "RHS") + fig_phi = plot_slices(phi, phi_exact, "Phi") + + rel_err_rhs = np.max(np.abs(rhs - rhs_exact)) / np.max(np.abs(rhs_exact)) + rel_err_phi = np.max(np.abs(phi - phi_exact)) / np.max(np.abs(phi_exact)) + + print(f"Max relative error in RHS: {rel_err_rhs:.2e}") + print(f"Max relative error in Phi: {rel_err_phi:.2e}") + + assert rel_err_rhs < 2.2e-2, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." + assert rel_err_phi < 1.3e-2, ( + f"The computed solution does not match the exact solution, max rel error = {rel_err_phi}." + ) + + import os + + # `out.path_out` is the out's (absolute) output folder; `sim_folder` alone is a bare + # name resolved against the CWD. The profiling packaging picks these files up from + # here and uploads them as `results-out`. + results_dir = os.path.join(out.path_out, "results") + os.makedirs(results_dir, exist_ok=True) + + np.save(os.path.join(results_dir, "rel_err_rhs.npy"), rel_err_rhs) + np.save(os.path.join(results_dir, "rel_err_phi.npy"), rel_err_phi) + np.save(os.path.join(results_dir, "resolution.npy"), sim.grid.num_elements) + np.save(os.path.join(results_dir, "spline_degree.npy"), sim.derham_opts.degree) + + if save_figs: + fig_rhs.savefig(os.path.join(results_dir, "rhs_slices.png")) + fig_phi.savefig(os.path.join(results_dir, "phi_slices.png")) + else: + plt.show() diff --git a/profiling/examples/Poisson/hollow_torus_strong_scaling_stiffness/params_poisson.py b/profiling/examples/Poisson/hollow_torus_strong_scaling_stiffness/params_poisson.py new file mode 100644 index 000000000..784face9c --- /dev/null +++ b/profiling/examples/Poisson/hollow_torus_strong_scaling_stiffness/params_poisson.py @@ -0,0 +1,289 @@ +# ----------------------------- +# Description of the simulation +# ----------------------------- +# Please fill in a verbal description of the simulation. +# It will be printed at the beginning of the simulation and can be used to keep track of the different runs. + +name = "Poisson strong scaling on hollow torus, stiffness preconditioner" +description = """ +Strong scaling test for Poisson equation on a full hollow torus (HollowTorus, a1=1e-2, a2=1, R0=4). +The manufactured solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5), +exciting the full spectrum. +Homogeneous Dirichlet boundary conditions are set in direction eta1. +The linear system is solved with CG preconditioned by the Kronecker approximation of the inverse Laplacian (StiffnessPreconditioner). +""" + +import logging + +from struphy import set_logging_level + +set_logging_level(logging.INFO) + +import argparse + +# ------------------ +# Import Struphy API +# ------------------ +from struphy import ( + BaseUnits, + DerhamOptions, + EnvironmentOptions, + FieldsBackground, + ProfilingOptions, + Simulation, + Time, + domains, + equils, + grids, + perturbations, +) + +# --------------------- +# Instance of the model +# --------------------- +from struphy.models import Poisson + +# Units +base_units = BaseUnits() + +# Model instance +model = Poisson(base_units=base_units) + +# List all variables and decide whether to save their data +model.em_fields.phi.save_data = True +model.em_fields.source.save_data = True + +# -------------------------- +# Instance of the simulation +# -------------------------- + +# Environment options +# `--id` distinguishes runs that share a rank count but differ in something else; the +# profiling driver passes its launch counter (see `ProfilingJob.build_commands`). +# Unknown flags are ignored so the driver can forward other parameters as well. +parser = argparse.ArgumentParser() +parser.add_argument("--id", type=int, default=0, help="Run id, used to name the output folder.") +args, _ = parser.parse_known_args() + +env = EnvironmentOptions( + sim_folder=f"sim_{args.id:02d}", + restart=False, +) + +# Time stepping +time_opts = Time() + +# Geometry +a1 = 1e-2 +a2 = 1.0 +R0 = 4.0 +tor_period = 1 +domain = domains.HollowTorus(a1=a1, a2=a2, R0=R0, sfl=False, pol_period=1, tor_period=tor_period) + +# Fluid equilibrium (can be used as part of initial conditions) +equil = None + +# Grid +grid = grids.TensorProductGrid(num_elements=(64, 64, 64), mpi_dims_mask=(True, True, True)) + +# Derham options +derham_opts = DerhamOptions(degree=(1, 2, 3), bcs=(("dirichlet", "dirichlet"), None, None)) + +# Profilinig options +profiling_opts = ProfilingOptions( + use_line_profiler=True, +) + +# Simulation object +sim = Simulation( + model=model, + name=name, + description=description, + params_path=__file__, + env=env, + time_opts=time_opts, + domain=domain, + equil=equil, + grid=grid, + derham_opts=derham_opts, + profiling_opts=profiling_opts, +) + +# ------------------ +# Propagator options +# ------------------ + +from struphy.linear_algebra.solver import SolverParameters + +solver_params = SolverParameters(tol=1e-8, maxiter=3000, info=True, recycle=False) +model.propagators.poisson.options = model.propagators.poisson.Options( + stab_eps=0.0, + solver="pcg", + precond="StiffnessPreconditioner", + solver_params=solver_params, +) + +# ------------------ +# Initial conditions +# ------------------ +import numpy as np + +from struphy.initial.base import GenericPerturbation + +# The exact solution is a Gaussian 0-form on the logical domain, centered at eta = (0.5, 0.5, 0.5). +# Its value at the Dirichlet boundary eta1 = 0, 1 is exp(-0.25 / w**2) ~ 1e-11, and it is periodic +# in eta2 and eta3 up to the same accuracy. Being localized, it excites the whole spectrum of the +# Laplacian, so CG has to work for convergence. Both the solution and the right-hand side +# -Laplace(phi) are 0-forms given as functions of the logical coordinates (given_in_basis="0"); +# the pushed-forward solution at x = F(eta) is phi(eta). +w = 0.1 +eta0 = (0.5, 0.5, 0.5) + + +def gaussian_derivatives(e1, e2, e3): + """The Gaussian and its first and second derivatives with respect to eta1, eta2, eta3.""" + eta = (e1, e2, e3) + g = np.exp(-sum((e - c) ** 2 for e, c in zip(eta, eta0)) / w**2) + d = [-2 * (e - c) / w**2 * g for e, c in zip(eta, eta0)] + dd = [(4 * (e - c) ** 2 / w**4 - 2 / w**2) * g for e, c in zip(eta, eta0)] + return g, d, dd + + +def exact_solution(e1, e2, e3): + return gaussian_derivatives(e1, e2, e3)[0] + + +def laplacian(e1, e2, e3): + """Physical Laplacian of the exact solution, in logical coordinates. + + HollowTorus (sfl=False): r = a1 + (a2 - a1) * eta1, theta = 2*pi*eta2, phi = 2*pi*eta3 / tor_period, + R = R0 + r*cos(theta), metric dr**2 + r**2 dtheta**2 + R**2 dphi**2, and + Laplace(u) = u_rr + (R + r*cos(theta)) / (r*R) * u_r + u_theta,theta / r**2 - sin(theta) / (r*R) * u_theta + + u_phi,phi / R**2. + """ + g, d, dd = gaussian_derivatives(e1, e2, e3) + da = a2 - a1 + r = a1 + da * e1 + theta = 2 * np.pi * e2 + R = R0 + r * np.cos(theta) + dphi = 2 * np.pi / tor_period + return ( + dd[0] / da**2 + + (R + r * np.cos(theta)) / (r * R) * d[0] / da + + dd[1] / (2 * np.pi * r) ** 2 + - np.sin(theta) / (r * R) * d[1] / (2 * np.pi) + + dd[2] / (dphi * R) ** 2 + ) + + +def rhs_fun(e1, e2, e3): + return -laplacian(e1, e2, e3) + + +rhs_perturbation = GenericPerturbation(rhs_fun, given_in_basis="0") + +model.em_fields.source.add_perturbation(rhs_perturbation) + + +if __name__ == "__main__": + # switches for debugging and testing + estimate_mem = False + run = True + pproc = True + save_figs = True + + if estimate_mem: + sim.estimate_mem(print_report=True) + exit() + + if run: + out = sim.run(profiling_activated=True, one_time_step=True) + else: + # use this for working with existing sim data from a previous run + out = sim.output + + if pproc: + out.pproc(create_vtk=True, parallel=True) + + from matplotlib import pyplot as plt + + def slice_axes(k, X, Y, Z, E1, E2, E3): + """Axes (horizontal, vertical, labels) for plotting the slice eta_k = 0.5: the poloidal plane (R, z) + at eta3 = 0.5 and the logical planes at eta1 = 0.5 and eta2 = 0.5.""" + if k == 0: + return E2, E3, "eta2", "eta3" + elif k == 1: + return E1, E3, "eta1", "eta3" + else: + return np.sqrt(X**2 + Y**2), Z, "R", "z" + + def plot_slices(num, exact, name): + """Numerical solution, exact solution and error on the slices eta_k = 0.5 (k = 1, 2, 3).""" + fig = plt.figure(figsize=(16, 12)) + for k in range(3): + idx = tuple(num.shape[k] // 2 if i == k else slice(None) for i in range(3)) + h, v, hlabel, vlabel = slice_axes(k, X[idx], Y[idx], Z[idx], E1[idx], E2[idx], E3[idx]) + eta_k = (E1, E2, E3)[k][idx].flat[0] + for row, (data, what) in enumerate( + [(num[idx], "from struphy"), (exact[idx], "exact"), (np.abs(num[idx] - exact[idx]), "error")] + ): + plt.subplot(3, 3, 3 * row + k + 1) + plt.pcolormesh(h, v, data, shading="gouraud") + plt.colorbar() + plt.xlabel(hlabel) + plt.ylabel(vlabel) + plt.title(f"{name} {what}, slice at eta{k + 1} = {eta_k:.2f}") + fig.tight_layout() + return fig + + if sim.rank == 0: + # Raw FEEC fields are evaluated directly from the saved spline coefficients at the + # cell centres of the simulation grid (serial Derham, safe on rank 0 only). + rhs = out.evaluate("em_fields/source", t=0).values + + phi_data = out.evaluate("em_fields/phi", t=-1) + phi = phi_data.values + + # logical coordinates of the evaluation points and their images x = F(eta) + E1, E2, E3 = np.meshgrid( + *(phi_data.coords[c].values for c in ("eta1", "eta2", "eta3")), + indexing="ij", + ) + X, Y, Z = (phi_data.coords[c].values for c in ("X", "Y", "Z")) + + # pushed-forward exact solutions: the 0-forms at x = F(eta) are their values at eta + rhs_exact = rhs_fun(E1, E2, E3) + phi_exact = exact_solution(E1, E2, E3) + + fig_rhs = plot_slices(rhs, rhs_exact, "RHS") + fig_phi = plot_slices(phi, phi_exact, "Phi") + + rel_err_rhs = np.max(np.abs(rhs - rhs_exact)) / np.max(np.abs(rhs_exact)) + rel_err_phi = np.max(np.abs(phi - phi_exact)) / np.max(np.abs(phi_exact)) + + print(f"Max relative error in RHS: {rel_err_rhs:.2e}") + print(f"Max relative error in Phi: {rel_err_phi:.2e}") + + assert rel_err_rhs < 2.2e-2, f"The computed RHS does not match the exact RHS, max rel error = {rel_err_rhs}." + assert rel_err_phi < 1.3e-2, ( + f"The computed solution does not match the exact solution, max rel error = {rel_err_phi}." + ) + + import os + + # `out.path_out` is the out's (absolute) output folder; `sim_folder` alone is a bare + # name resolved against the CWD. The profiling packaging picks these files up from + # here and uploads them as `results-out`. + results_dir = os.path.join(out.path_out, "results") + os.makedirs(results_dir, exist_ok=True) + + np.save(os.path.join(results_dir, "rel_err_rhs.npy"), rel_err_rhs) + np.save(os.path.join(results_dir, "rel_err_phi.npy"), rel_err_phi) + np.save(os.path.join(results_dir, "resolution.npy"), sim.grid.num_elements) + np.save(os.path.join(results_dir, "spline_degree.npy"), sim.derham_opts.degree) + + if save_figs: + fig_rhs.savefig(os.path.join(results_dir, "rhs_slices.png")) + fig_phi.savefig(os.path.join(results_dir, "phi_slices.png")) + else: + plt.show() diff --git a/profiling/submit_poisson_strong_scaling.py b/profiling/examples/Poisson/submit_poisson_cube_strong_scaling.py similarity index 65% rename from profiling/submit_poisson_strong_scaling.py rename to profiling/examples/Poisson/submit_poisson_cube_strong_scaling.py index d8387199a..b41fc652a 100644 --- a/profiling/submit_poisson_strong_scaling.py +++ b/profiling/examples/Poisson/submit_poisson_cube_strong_scaling.py @@ -1,4 +1,4 @@ -"""Poisson strong scaling profiling case. +"""Poisson strong scaling profiling case on 3D cube, no preconditioner. This file defines the Poisson strong scaling profiling case (the `ProfilingCase`) and submits it: for each rank count, `ProfilingCase.launch` builds and submits a @@ -10,9 +10,13 @@ """ import argparse +import sys from pathlib import Path -from profiling_job import ProfilingCase +# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/, two levels up +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + +from profiling_job import ProfilingCase # noqa: E402 def main() -> None: @@ -26,21 +30,32 @@ def main() -> None: action="store_true", help="Upload the packaged profiling results to the profiling-data repo.", ) + parser.add_argument( + "--language", + default="fortran", + choices=["fortran", "c"], + help="Pyccel language to compile the Struphy kernels with.", + ) + parser.add_argument( + "--compiler", + default="GNU", + help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").', + ) args = parser.parse_args() # Paths relative to this script's location, so it can be run from anywhere. script_dir = Path(__file__).resolve().parent - params_dir = script_dir / "examples" / "Poisson" / "cube_strong_scaling" + params_dir = script_dir / "cube_strong_scaling" profiling_case = ProfilingCase( label="poisson_cube_strong_scaling", - name="Poisson on cube strong scaling test, no preconditioner", - description="Strong scaling of the Poisson model with manufactured solution on 3D cube, solved with unpreconditioned CG.", + name="Poisson on 3D cube strong scaling test, no preconditioner", + description="Strong scaling of the Poisson model with manufactured solution on a 3D cube (Cuboid), solved with unpreconditioned CG.", physics_problem="Occurs in many plasma applications.", struphy_model_used="Poisson", params_source=params_dir / "params_poisson.py", - language="fortran", - compiler="GNU", + language=args.language, + compiler=args.compiler, upload=args.upload, ) diff --git a/profiling/examples/Poisson/submit_poisson_cube_strong_scaling_massmatrix.py b/profiling/examples/Poisson/submit_poisson_cube_strong_scaling_massmatrix.py new file mode 100644 index 000000000..67e49e47a --- /dev/null +++ b/profiling/examples/Poisson/submit_poisson_cube_strong_scaling_massmatrix.py @@ -0,0 +1,71 @@ +"""Poisson strong scaling profiling case on 3D cube, mass-matrix preconditioner. + +This file defines the Poisson strong scaling profiling case (the `ProfilingCase`) +and submits it: for each rank count, `ProfilingCase.launch` builds and submits a +SLURM script (using `clusters.SLURM_PRESETS` by default), or, without a batch +system, runs directly on this machine. `finalize_run` then packages and uploads +each run as soon as its own job finishes. +Each generated script runs the simulation itself by invoking `params_poisson.py` +directly (its `__main__` block is the worker). +""" + +import argparse +import sys +from pathlib import Path + +# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/, two levels up +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + +from profiling_job import ProfilingCase # noqa: E402 + + +def main() -> None: + + # Parse arguments, do not remove --upload + parser = argparse.ArgumentParser( + description=("Submit profiling jobs to a SLURM cluster and package the results for upload."), + ) + parser.add_argument( + "--upload", + action="store_true", + help="Upload the packaged profiling results to the profiling-data repo.", + ) + parser.add_argument( + "--language", + default="fortran", + choices=["fortran", "c"], + help="Pyccel language to compile the Struphy kernels with.", + ) + parser.add_argument( + "--compiler", + default="GNU", + help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").', + ) + args = parser.parse_args() + + # Paths relative to this script's location, so it can be run from anywhere. + script_dir = Path(__file__).resolve().parent + params_dir = script_dir / "cube_strong_scaling_massmatrix" + + profiling_case = ProfilingCase( + label="poisson_cube_strong_scaling_massmatrix", + name="Poisson on 3D cube strong scaling test, mass-matrix preconditioner", + description="Strong scaling of the Poisson model with manufactured solution on a 3D cube (Cuboid), solved with CG preconditioned by the Kronecker approximation of the inverse 0-form mass matrix (MassMatrixPreconditioner).", + physics_problem="Occurs in many plasma applications.", + struphy_model_used="Poisson", + params_source=params_dir / "params_poisson.py", + language=args.language, + compiler=args.compiler, + upload=args.upload, + ) + + # Launch one run per rank count + for num_tasks in (1, 2, 4, 8, 16, 32, 64, 128, 256): + profiling_case.launch(num_tasks) + + # Package and push each run as its own job finishes. + profiling_case.finalize_run() + + +if __name__ == "__main__": + main() diff --git a/profiling/submit_poisson_multigrid_strong_scaling.py b/profiling/examples/Poisson/submit_poisson_cube_strong_scaling_multigrid.py similarity index 64% rename from profiling/submit_poisson_multigrid_strong_scaling.py rename to profiling/examples/Poisson/submit_poisson_cube_strong_scaling_multigrid.py index 4cc991078..684ba2c41 100644 --- a/profiling/submit_poisson_multigrid_strong_scaling.py +++ b/profiling/examples/Poisson/submit_poisson_cube_strong_scaling_multigrid.py @@ -1,4 +1,4 @@ -"""Poisson strong scaling profiling case with multigrid preconditioner. +"""Poisson strong scaling profiling case on 3D cube, multigrid preconditioner. This file defines the Poisson strong scaling profiling case (the `ProfilingCase`) and submits it: for each rank count, `ProfilingCase.launch` builds and submits a @@ -10,9 +10,13 @@ """ import argparse +import sys from pathlib import Path -from profiling_job import ProfilingCase +# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/, two levels up +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + +from profiling_job import ProfilingCase # noqa: E402 def main() -> None: @@ -26,21 +30,32 @@ def main() -> None: action="store_true", help="Upload the packaged profiling results to the profiling-data repo.", ) + parser.add_argument( + "--language", + default="fortran", + choices=["fortran", "c"], + help="Pyccel language to compile the Struphy kernels with.", + ) + parser.add_argument( + "--compiler", + default="GNU", + help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").', + ) args = parser.parse_args() # Paths relative to this script's location, so it can be run from anywhere. script_dir = Path(__file__).resolve().parent - params_dir = script_dir / "examples" / "Poisson" / "cube_strong_scaling_multigrid" + params_dir = script_dir / "cube_strong_scaling_multigrid" profiling_case = ProfilingCase( label="poisson_cube_strong_scaling_multigrid", - name="Poisson on cube strong scaling test, multigrid preconditioner", - description="Strong scaling of the Poisson model with manufactured solution on 3D cube, solved with multigrid-preconditioned CG.", + name="Poisson on 3D cube strong scaling test, multigrid preconditioner", + description="Strong scaling of the Poisson model with manufactured solution on a 3D cube (Cuboid), solved with CG preconditioned by geometric multigrid.", physics_problem="Occurs in many plasma applications.", struphy_model_used="Poisson", params_source=params_dir / "params_poisson.py", - language="fortran", - compiler="GNU", + language=args.language, + compiler=args.compiler, upload=args.upload, ) diff --git a/profiling/examples/Poisson/submit_poisson_cube_strong_scaling_stiffness.py b/profiling/examples/Poisson/submit_poisson_cube_strong_scaling_stiffness.py new file mode 100644 index 000000000..82909ad33 --- /dev/null +++ b/profiling/examples/Poisson/submit_poisson_cube_strong_scaling_stiffness.py @@ -0,0 +1,71 @@ +"""Poisson strong scaling profiling case on 3D cube, stiffness preconditioner. + +This file defines the Poisson strong scaling profiling case (the `ProfilingCase`) +and submits it: for each rank count, `ProfilingCase.launch` builds and submits a +SLURM script (using `clusters.SLURM_PRESETS` by default), or, without a batch +system, runs directly on this machine. `finalize_run` then packages and uploads +each run as soon as its own job finishes. +Each generated script runs the simulation itself by invoking `params_poisson.py` +directly (its `__main__` block is the worker). +""" + +import argparse +import sys +from pathlib import Path + +# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/, two levels up +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + +from profiling_job import ProfilingCase # noqa: E402 + + +def main() -> None: + + # Parse arguments, do not remove --upload + parser = argparse.ArgumentParser( + description=("Submit profiling jobs to a SLURM cluster and package the results for upload."), + ) + parser.add_argument( + "--upload", + action="store_true", + help="Upload the packaged profiling results to the profiling-data repo.", + ) + parser.add_argument( + "--language", + default="fortran", + choices=["fortran", "c"], + help="Pyccel language to compile the Struphy kernels with.", + ) + parser.add_argument( + "--compiler", + default="GNU", + help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").', + ) + args = parser.parse_args() + + # Paths relative to this script's location, so it can be run from anywhere. + script_dir = Path(__file__).resolve().parent + params_dir = script_dir / "cube_strong_scaling_stiffness" + + profiling_case = ProfilingCase( + label="poisson_cube_strong_scaling_stiffness", + name="Poisson on 3D cube strong scaling test, stiffness preconditioner", + description="Strong scaling of the Poisson model with manufactured solution on a 3D cube (Cuboid), solved with CG preconditioned by the Kronecker approximation of the inverse Laplacian (StiffnessPreconditioner).", + physics_problem="Occurs in many plasma applications.", + struphy_model_used="Poisson", + params_source=params_dir / "params_poisson.py", + language=args.language, + compiler=args.compiler, + upload=args.upload, + ) + + # Launch one run per rank count + for num_tasks in (1, 2, 4, 8, 16, 32, 64, 128, 256): + profiling_case.launch(num_tasks) + + # Package and push each run as its own job finishes. + profiling_case.finalize_run() + + +if __name__ == "__main__": + main() diff --git a/profiling/examples/Poisson/submit_poisson_hollow_cylinder_strong_scaling.py b/profiling/examples/Poisson/submit_poisson_hollow_cylinder_strong_scaling.py new file mode 100644 index 000000000..d98ad53a3 --- /dev/null +++ b/profiling/examples/Poisson/submit_poisson_hollow_cylinder_strong_scaling.py @@ -0,0 +1,71 @@ +"""Poisson strong scaling profiling case on hollow cylinder, no preconditioner. + +This file defines the Poisson strong scaling profiling case (the `ProfilingCase`) +and submits it: for each rank count, `ProfilingCase.launch` builds and submits a +SLURM script (using `clusters.SLURM_PRESETS` by default), or, without a batch +system, runs directly on this machine. `finalize_run` then packages and uploads +each run as soon as its own job finishes. +Each generated script runs the simulation itself by invoking `params_poisson.py` +directly (its `__main__` block is the worker). +""" + +import argparse +import sys +from pathlib import Path + +# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/, two levels up +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + +from profiling_job import ProfilingCase # noqa: E402 + + +def main() -> None: + + # Parse arguments, do not remove --upload + parser = argparse.ArgumentParser( + description=("Submit profiling jobs to a SLURM cluster and package the results for upload."), + ) + parser.add_argument( + "--upload", + action="store_true", + help="Upload the packaged profiling results to the profiling-data repo.", + ) + parser.add_argument( + "--language", + default="fortran", + choices=["fortran", "c"], + help="Pyccel language to compile the Struphy kernels with.", + ) + parser.add_argument( + "--compiler", + default="GNU", + help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").', + ) + args = parser.parse_args() + + # Paths relative to this script's location, so it can be run from anywhere. + script_dir = Path(__file__).resolve().parent + params_dir = script_dir / "hollow_cylinder_strong_scaling" + + profiling_case = ProfilingCase( + label="poisson_hollow_cylinder_strong_scaling", + name="Poisson on hollow cylinder strong scaling test, no preconditioner", + description="Strong scaling of the Poisson model with manufactured solution on a hollow cylinder (HollowCylinder, a1=1e-2, a2=1, Lz=4), solved with unpreconditioned CG.", + physics_problem="Occurs in many plasma applications.", + struphy_model_used="Poisson", + params_source=params_dir / "params_poisson.py", + language=args.language, + compiler=args.compiler, + upload=args.upload, + ) + + # Launch one run per rank count + for num_tasks in (1, 2, 4, 8, 16, 32, 64, 128, 256): + profiling_case.launch(num_tasks) + + # Package and push each run as its own job finishes. + profiling_case.finalize_run() + + +if __name__ == "__main__": + main() diff --git a/profiling/examples/Poisson/submit_poisson_hollow_cylinder_strong_scaling_massmatrix.py b/profiling/examples/Poisson/submit_poisson_hollow_cylinder_strong_scaling_massmatrix.py new file mode 100644 index 000000000..7952e3b60 --- /dev/null +++ b/profiling/examples/Poisson/submit_poisson_hollow_cylinder_strong_scaling_massmatrix.py @@ -0,0 +1,71 @@ +"""Poisson strong scaling profiling case on hollow cylinder, mass-matrix preconditioner. + +This file defines the Poisson strong scaling profiling case (the `ProfilingCase`) +and submits it: for each rank count, `ProfilingCase.launch` builds and submits a +SLURM script (using `clusters.SLURM_PRESETS` by default), or, without a batch +system, runs directly on this machine. `finalize_run` then packages and uploads +each run as soon as its own job finishes. +Each generated script runs the simulation itself by invoking `params_poisson.py` +directly (its `__main__` block is the worker). +""" + +import argparse +import sys +from pathlib import Path + +# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/, two levels up +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + +from profiling_job import ProfilingCase # noqa: E402 + + +def main() -> None: + + # Parse arguments, do not remove --upload + parser = argparse.ArgumentParser( + description=("Submit profiling jobs to a SLURM cluster and package the results for upload."), + ) + parser.add_argument( + "--upload", + action="store_true", + help="Upload the packaged profiling results to the profiling-data repo.", + ) + parser.add_argument( + "--language", + default="fortran", + choices=["fortran", "c"], + help="Pyccel language to compile the Struphy kernels with.", + ) + parser.add_argument( + "--compiler", + default="GNU", + help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").', + ) + args = parser.parse_args() + + # Paths relative to this script's location, so it can be run from anywhere. + script_dir = Path(__file__).resolve().parent + params_dir = script_dir / "hollow_cylinder_strong_scaling_massmatrix" + + profiling_case = ProfilingCase( + label="poisson_hollow_cylinder_strong_scaling_massmatrix", + name="Poisson on hollow cylinder strong scaling test, mass-matrix preconditioner", + description="Strong scaling of the Poisson model with manufactured solution on a hollow cylinder (HollowCylinder, a1=1e-2, a2=1, Lz=4), solved with CG preconditioned by the Kronecker approximation of the inverse 0-form mass matrix (MassMatrixPreconditioner).", + physics_problem="Occurs in many plasma applications.", + struphy_model_used="Poisson", + params_source=params_dir / "params_poisson.py", + language=args.language, + compiler=args.compiler, + upload=args.upload, + ) + + # Launch one run per rank count + for num_tasks in (1, 2, 4, 8, 16, 32, 64, 128, 256): + profiling_case.launch(num_tasks) + + # Package and push each run as its own job finishes. + profiling_case.finalize_run() + + +if __name__ == "__main__": + main() diff --git a/profiling/examples/Poisson/submit_poisson_hollow_cylinder_strong_scaling_multigrid.py b/profiling/examples/Poisson/submit_poisson_hollow_cylinder_strong_scaling_multigrid.py new file mode 100644 index 000000000..85b895335 --- /dev/null +++ b/profiling/examples/Poisson/submit_poisson_hollow_cylinder_strong_scaling_multigrid.py @@ -0,0 +1,71 @@ +"""Poisson strong scaling profiling case on hollow cylinder, multigrid preconditioner. + +This file defines the Poisson strong scaling profiling case (the `ProfilingCase`) +and submits it: for each rank count, `ProfilingCase.launch` builds and submits a +SLURM script (using `clusters.SLURM_PRESETS` by default), or, without a batch +system, runs directly on this machine. `finalize_run` then packages and uploads +each run as soon as its own job finishes. +Each generated script runs the simulation itself by invoking `params_poisson.py` +directly (its `__main__` block is the worker). +""" + +import argparse +import sys +from pathlib import Path + +# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/, two levels up +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + +from profiling_job import ProfilingCase # noqa: E402 + + +def main() -> None: + + # Parse arguments, do not remove --upload + parser = argparse.ArgumentParser( + description=("Submit profiling jobs to a SLURM cluster and package the results for upload."), + ) + parser.add_argument( + "--upload", + action="store_true", + help="Upload the packaged profiling results to the profiling-data repo.", + ) + parser.add_argument( + "--language", + default="fortran", + choices=["fortran", "c"], + help="Pyccel language to compile the Struphy kernels with.", + ) + parser.add_argument( + "--compiler", + default="GNU", + help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").', + ) + args = parser.parse_args() + + # Paths relative to this script's location, so it can be run from anywhere. + script_dir = Path(__file__).resolve().parent + params_dir = script_dir / "hollow_cylinder_strong_scaling_multigrid" + + profiling_case = ProfilingCase( + label="poisson_hollow_cylinder_strong_scaling_multigrid", + name="Poisson on hollow cylinder strong scaling test, multigrid preconditioner", + description="Strong scaling of the Poisson model with manufactured solution on a hollow cylinder (HollowCylinder, a1=1e-2, a2=1, Lz=4), solved with CG preconditioned by geometric multigrid.", + physics_problem="Occurs in many plasma applications.", + struphy_model_used="Poisson", + params_source=params_dir / "params_poisson.py", + language=args.language, + compiler=args.compiler, + upload=args.upload, + ) + + # Launch one run per rank count + for num_tasks in (1, 2, 4, 8, 16, 32, 64, 128, 256): + profiling_case.launch(num_tasks) + + # Package and push each run as its own job finishes. + profiling_case.finalize_run() + + +if __name__ == "__main__": + main() diff --git a/profiling/examples/Poisson/submit_poisson_hollow_cylinder_strong_scaling_stiffness.py b/profiling/examples/Poisson/submit_poisson_hollow_cylinder_strong_scaling_stiffness.py new file mode 100644 index 000000000..d9ce5c9f3 --- /dev/null +++ b/profiling/examples/Poisson/submit_poisson_hollow_cylinder_strong_scaling_stiffness.py @@ -0,0 +1,71 @@ +"""Poisson strong scaling profiling case on hollow cylinder, stiffness preconditioner. + +This file defines the Poisson strong scaling profiling case (the `ProfilingCase`) +and submits it: for each rank count, `ProfilingCase.launch` builds and submits a +SLURM script (using `clusters.SLURM_PRESETS` by default), or, without a batch +system, runs directly on this machine. `finalize_run` then packages and uploads +each run as soon as its own job finishes. +Each generated script runs the simulation itself by invoking `params_poisson.py` +directly (its `__main__` block is the worker). +""" + +import argparse +import sys +from pathlib import Path + +# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/, two levels up +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + +from profiling_job import ProfilingCase # noqa: E402 + + +def main() -> None: + + # Parse arguments, do not remove --upload + parser = argparse.ArgumentParser( + description=("Submit profiling jobs to a SLURM cluster and package the results for upload."), + ) + parser.add_argument( + "--upload", + action="store_true", + help="Upload the packaged profiling results to the profiling-data repo.", + ) + parser.add_argument( + "--language", + default="fortran", + choices=["fortran", "c"], + help="Pyccel language to compile the Struphy kernels with.", + ) + parser.add_argument( + "--compiler", + default="GNU", + help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").', + ) + args = parser.parse_args() + + # Paths relative to this script's location, so it can be run from anywhere. + script_dir = Path(__file__).resolve().parent + params_dir = script_dir / "hollow_cylinder_strong_scaling_stiffness" + + profiling_case = ProfilingCase( + label="poisson_hollow_cylinder_strong_scaling_stiffness", + name="Poisson on hollow cylinder strong scaling test, stiffness preconditioner", + description="Strong scaling of the Poisson model with manufactured solution on a hollow cylinder (HollowCylinder, a1=1e-2, a2=1, Lz=4), solved with CG preconditioned by the Kronecker approximation of the inverse Laplacian (StiffnessPreconditioner).", + physics_problem="Occurs in many plasma applications.", + struphy_model_used="Poisson", + params_source=params_dir / "params_poisson.py", + language=args.language, + compiler=args.compiler, + upload=args.upload, + ) + + # Launch one run per rank count + for num_tasks in (1, 2, 4, 8, 16, 32, 64, 128, 256): + profiling_case.launch(num_tasks) + + # Package and push each run as its own job finishes. + profiling_case.finalize_run() + + +if __name__ == "__main__": + main() diff --git a/profiling/examples/Poisson/submit_poisson_hollow_torus_strong_scaling.py b/profiling/examples/Poisson/submit_poisson_hollow_torus_strong_scaling.py new file mode 100644 index 000000000..3d1b1c605 --- /dev/null +++ b/profiling/examples/Poisson/submit_poisson_hollow_torus_strong_scaling.py @@ -0,0 +1,71 @@ +"""Poisson strong scaling profiling case on hollow torus, no preconditioner. + +This file defines the Poisson strong scaling profiling case (the `ProfilingCase`) +and submits it: for each rank count, `ProfilingCase.launch` builds and submits a +SLURM script (using `clusters.SLURM_PRESETS` by default), or, without a batch +system, runs directly on this machine. `finalize_run` then packages and uploads +each run as soon as its own job finishes. +Each generated script runs the simulation itself by invoking `params_poisson.py` +directly (its `__main__` block is the worker). +""" + +import argparse +import sys +from pathlib import Path + +# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/, two levels up +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + +from profiling_job import ProfilingCase # noqa: E402 + + +def main() -> None: + + # Parse arguments, do not remove --upload + parser = argparse.ArgumentParser( + description=("Submit profiling jobs to a SLURM cluster and package the results for upload."), + ) + parser.add_argument( + "--upload", + action="store_true", + help="Upload the packaged profiling results to the profiling-data repo.", + ) + parser.add_argument( + "--language", + default="fortran", + choices=["fortran", "c"], + help="Pyccel language to compile the Struphy kernels with.", + ) + parser.add_argument( + "--compiler", + default="GNU", + help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").', + ) + args = parser.parse_args() + + # Paths relative to this script's location, so it can be run from anywhere. + script_dir = Path(__file__).resolve().parent + params_dir = script_dir / "hollow_torus_strong_scaling" + + profiling_case = ProfilingCase( + label="poisson_hollow_torus_strong_scaling", + name="Poisson on hollow torus strong scaling test, no preconditioner", + description="Strong scaling of the Poisson model with manufactured solution on a full hollow torus (HollowTorus, a1=1e-2, a2=1, R0=4), solved with unpreconditioned CG.", + physics_problem="Occurs in many plasma applications.", + struphy_model_used="Poisson", + params_source=params_dir / "params_poisson.py", + language=args.language, + compiler=args.compiler, + upload=args.upload, + ) + + # Launch one run per rank count + for num_tasks in (1, 2, 4, 8, 16, 32, 64, 128, 256): + profiling_case.launch(num_tasks) + + # Package and push each run as its own job finishes. + profiling_case.finalize_run() + + +if __name__ == "__main__": + main() diff --git a/profiling/examples/Poisson/submit_poisson_hollow_torus_strong_scaling_massmatrix.py b/profiling/examples/Poisson/submit_poisson_hollow_torus_strong_scaling_massmatrix.py new file mode 100644 index 000000000..cab233f62 --- /dev/null +++ b/profiling/examples/Poisson/submit_poisson_hollow_torus_strong_scaling_massmatrix.py @@ -0,0 +1,71 @@ +"""Poisson strong scaling profiling case on hollow torus, mass-matrix preconditioner. + +This file defines the Poisson strong scaling profiling case (the `ProfilingCase`) +and submits it: for each rank count, `ProfilingCase.launch` builds and submits a +SLURM script (using `clusters.SLURM_PRESETS` by default), or, without a batch +system, runs directly on this machine. `finalize_run` then packages and uploads +each run as soon as its own job finishes. +Each generated script runs the simulation itself by invoking `params_poisson.py` +directly (its `__main__` block is the worker). +""" + +import argparse +import sys +from pathlib import Path + +# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/, two levels up +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + +from profiling_job import ProfilingCase # noqa: E402 + + +def main() -> None: + + # Parse arguments, do not remove --upload + parser = argparse.ArgumentParser( + description=("Submit profiling jobs to a SLURM cluster and package the results for upload."), + ) + parser.add_argument( + "--upload", + action="store_true", + help="Upload the packaged profiling results to the profiling-data repo.", + ) + parser.add_argument( + "--language", + default="fortran", + choices=["fortran", "c"], + help="Pyccel language to compile the Struphy kernels with.", + ) + parser.add_argument( + "--compiler", + default="GNU", + help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").', + ) + args = parser.parse_args() + + # Paths relative to this script's location, so it can be run from anywhere. + script_dir = Path(__file__).resolve().parent + params_dir = script_dir / "hollow_torus_strong_scaling_massmatrix" + + profiling_case = ProfilingCase( + label="poisson_hollow_torus_strong_scaling_massmatrix", + name="Poisson on hollow torus strong scaling test, mass-matrix preconditioner", + description="Strong scaling of the Poisson model with manufactured solution on a full hollow torus (HollowTorus, a1=1e-2, a2=1, R0=4), solved with CG preconditioned by the Kronecker approximation of the inverse 0-form mass matrix (MassMatrixPreconditioner).", + physics_problem="Occurs in many plasma applications.", + struphy_model_used="Poisson", + params_source=params_dir / "params_poisson.py", + language=args.language, + compiler=args.compiler, + upload=args.upload, + ) + + # Launch one run per rank count + for num_tasks in (1, 2, 4, 8, 16, 32, 64, 128, 256): + profiling_case.launch(num_tasks) + + # Package and push each run as its own job finishes. + profiling_case.finalize_run() + + +if __name__ == "__main__": + main() diff --git a/profiling/examples/Poisson/submit_poisson_hollow_torus_strong_scaling_multigrid.py b/profiling/examples/Poisson/submit_poisson_hollow_torus_strong_scaling_multigrid.py new file mode 100644 index 000000000..2ebde34dd --- /dev/null +++ b/profiling/examples/Poisson/submit_poisson_hollow_torus_strong_scaling_multigrid.py @@ -0,0 +1,71 @@ +"""Poisson strong scaling profiling case on hollow torus, multigrid preconditioner. + +This file defines the Poisson strong scaling profiling case (the `ProfilingCase`) +and submits it: for each rank count, `ProfilingCase.launch` builds and submits a +SLURM script (using `clusters.SLURM_PRESETS` by default), or, without a batch +system, runs directly on this machine. `finalize_run` then packages and uploads +each run as soon as its own job finishes. +Each generated script runs the simulation itself by invoking `params_poisson.py` +directly (its `__main__` block is the worker). +""" + +import argparse +import sys +from pathlib import Path + +# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/, two levels up +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + +from profiling_job import ProfilingCase # noqa: E402 + + +def main() -> None: + + # Parse arguments, do not remove --upload + parser = argparse.ArgumentParser( + description=("Submit profiling jobs to a SLURM cluster and package the results for upload."), + ) + parser.add_argument( + "--upload", + action="store_true", + help="Upload the packaged profiling results to the profiling-data repo.", + ) + parser.add_argument( + "--language", + default="fortran", + choices=["fortran", "c"], + help="Pyccel language to compile the Struphy kernels with.", + ) + parser.add_argument( + "--compiler", + default="GNU", + help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").', + ) + args = parser.parse_args() + + # Paths relative to this script's location, so it can be run from anywhere. + script_dir = Path(__file__).resolve().parent + params_dir = script_dir / "hollow_torus_strong_scaling_multigrid" + + profiling_case = ProfilingCase( + label="poisson_hollow_torus_strong_scaling_multigrid", + name="Poisson on hollow torus strong scaling test, multigrid preconditioner", + description="Strong scaling of the Poisson model with manufactured solution on a full hollow torus (HollowTorus, a1=1e-2, a2=1, R0=4), solved with CG preconditioned by geometric multigrid.", + physics_problem="Occurs in many plasma applications.", + struphy_model_used="Poisson", + params_source=params_dir / "params_poisson.py", + language=args.language, + compiler=args.compiler, + upload=args.upload, + ) + + # Launch one run per rank count + for num_tasks in (1, 2, 4, 8, 16, 32, 64, 128, 256): + profiling_case.launch(num_tasks) + + # Package and push each run as its own job finishes. + profiling_case.finalize_run() + + +if __name__ == "__main__": + main() diff --git a/profiling/examples/Poisson/submit_poisson_hollow_torus_strong_scaling_stiffness.py b/profiling/examples/Poisson/submit_poisson_hollow_torus_strong_scaling_stiffness.py new file mode 100644 index 000000000..288630e20 --- /dev/null +++ b/profiling/examples/Poisson/submit_poisson_hollow_torus_strong_scaling_stiffness.py @@ -0,0 +1,71 @@ +"""Poisson strong scaling profiling case on hollow torus, stiffness preconditioner. + +This file defines the Poisson strong scaling profiling case (the `ProfilingCase`) +and submits it: for each rank count, `ProfilingCase.launch` builds and submits a +SLURM script (using `clusters.SLURM_PRESETS` by default), or, without a batch +system, runs directly on this machine. `finalize_run` then packages and uploads +each run as soon as its own job finishes. +Each generated script runs the simulation itself by invoking `params_poisson.py` +directly (its `__main__` block is the worker). +""" + +import argparse +import sys +from pathlib import Path + +# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/, two levels up +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + +from profiling_job import ProfilingCase # noqa: E402 + + +def main() -> None: + + # Parse arguments, do not remove --upload + parser = argparse.ArgumentParser( + description=("Submit profiling jobs to a SLURM cluster and package the results for upload."), + ) + parser.add_argument( + "--upload", + action="store_true", + help="Upload the packaged profiling results to the profiling-data repo.", + ) + parser.add_argument( + "--language", + default="fortran", + choices=["fortran", "c"], + help="Pyccel language to compile the Struphy kernels with.", + ) + parser.add_argument( + "--compiler", + default="GNU", + help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").', + ) + args = parser.parse_args() + + # Paths relative to this script's location, so it can be run from anywhere. + script_dir = Path(__file__).resolve().parent + params_dir = script_dir / "hollow_torus_strong_scaling_stiffness" + + profiling_case = ProfilingCase( + label="poisson_hollow_torus_strong_scaling_stiffness", + name="Poisson on hollow torus strong scaling test, stiffness preconditioner", + description="Strong scaling of the Poisson model with manufactured solution on a full hollow torus (HollowTorus, a1=1e-2, a2=1, R0=4), solved with CG preconditioned by the Kronecker approximation of the inverse Laplacian (StiffnessPreconditioner).", + physics_problem="Occurs in many plasma applications.", + struphy_model_used="Poisson", + params_source=params_dir / "params_poisson.py", + language=args.language, + compiler=args.compiler, + upload=args.upload, + ) + + # Launch one run per rank count + for num_tasks in (1, 2, 4, 8, 16, 32, 64, 128, 256): + profiling_case.launch(num_tasks) + + # Package and push each run as its own job finishes. + profiling_case.finalize_run() + + +if __name__ == "__main__": + main() diff --git a/profiling/submit_diocotron_strong_scaling.py b/profiling/examples/ToyGyrokinetic/submit_diocotron_strong_scaling.py similarity index 73% rename from profiling/submit_diocotron_strong_scaling.py rename to profiling/examples/ToyGyrokinetic/submit_diocotron_strong_scaling.py index 449f76fd0..cd0f930ff 100644 --- a/profiling/submit_diocotron_strong_scaling.py +++ b/profiling/examples/ToyGyrokinetic/submit_diocotron_strong_scaling.py @@ -12,9 +12,13 @@ """ import argparse +import sys from pathlib import Path -from profiling_job import ProfilingCase +# the profiling driver (profiling_job.py, clusters.py, ...) lives in profiling/, two levels up +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) + +from profiling_job import ProfilingCase # noqa: E402 def main() -> None: @@ -28,11 +32,22 @@ def main() -> None: action="store_true", help="Upload the packaged profiling results to the profiling-data repo.", ) + parser.add_argument( + "--language", + default="fortran", + choices=["fortran", "c"], + help="Pyccel language to compile the Struphy kernels with.", + ) + parser.add_argument( + "--compiler", + default="GNU", + help='Pyccel compiler family ("GNU", "intel", "PGI", "nvidia" or "LLVM").', + ) args = parser.parse_args() # Paths relative to this script's location, so it can be run from anywhere. script_dir = Path(__file__).resolve().parent - params_dir = script_dir / "examples" / "ToyGyrokinetic" / "diocotron_instability" + params_dir = script_dir / "diocotron_instability" profiling_case = ProfilingCase( label="diocotron_instability", @@ -41,8 +56,8 @@ def main() -> None: physics_problem="Diocotron instability in a non-neutral plasma.", struphy_model_used="ToyDrift", params_source=params_dir / "params_diocotron.py", - language="fortran", - compiler="GNU", + language=args.language, + compiler=args.compiler, upload=args.upload, ) diff --git a/src/struphy/feec/preconditioner.py b/src/struphy/feec/preconditioner.py index f671c39c2..1d3481b25 100644 --- a/src/struphy/feec/preconditioner.py +++ b/src/struphy/feec/preconditioner.py @@ -262,7 +262,11 @@ class MassMatrixPreconditioner(KroneckerPreconditioner): Whether to include boundary operators. dim_reduce : int | None - Axis along which the weight is kept; None for unit weights in all directions. + Axis along which the weight is kept (default 0); None for unit weights in all directions. + With an axis, the approximation is exact for weights (and mappings) that vary only along + that axis, e.g. radial profiles in eta1. With None, the Kronecker factors do not depend + on the weights, so ``update_mass_operator`` only recomputes the diagonal scaling: use it + for weights that change often (e.g. a density updated in every time step). weight_reduction : str How the weight is reduced in the other axes: ``"midpoint"`` (value at 0.5, default) @@ -271,7 +275,13 @@ class MassMatrixPreconditioner(KroneckerPreconditioner): diagonal_scaling : bool Whether to correct the approximation with the diagonals of the mass matrix and of - its approximation (see :class:`KroneckerPreconditioner`). + its approximation (see :class:`KroneckerPreconditioner`; default True). + + Notes + ----- + The diagonal scaling follows `G. Loli, G. Sangalli, M. Tani, "Easy and efficient preconditioning + of the isogeometric mass matrix", Comp. Math. Appl., Vol. 116, 2022 + `_. """ def __init__( @@ -280,7 +290,7 @@ def __init__( apply_bc: bool = True, dim_reduce: int | None = 0, weight_reduction: str = "midpoint", - diagonal_scaling: bool = False, + diagonal_scaling: bool = True, ): assert dim_reduce is None or dim_reduce < 3 assert weight_reduction in ("midpoint", "average"), f"Unknown weight_reduction {weight_reduction!r}." @@ -350,43 +360,6 @@ def transpose(self, conjugate: bool = False) -> MassMatrixPreconditioner: ) -class MassMatrixDiagonalPreconditioner(MassMatrixPreconditioner): - r""" - Preconditioner for inverting 3d weighted mass matrices. The mass matrix is approximated by - - .. math:: - D^{1/2} * \hat D^{-1/2} * \hat M * \hat D^{-1/2} * D^{1/2} - - Where $D$ is the diagonal of the matrix to invert, :math:`\hat M` is the mass matrix on the logical domain - that is a Kronecker product (fastly inverted) and :math:`\hat D^{-1/2}` is the diagonal of :math:`\hat M`. - - This is ``MassMatrixPreconditioner(mass_operator, apply_bc, dim_reduce=None, diagonal_scaling=True)``; - the class is kept for its name (e.g. in solver options). - - Notes - ----- - - Reference: `G. Loli, G. Sangalli, M. Tani, "Easy and efficient preconditioning of the isogeometric mass matrix", Comp. Math. Appl., Vol. 116, 2022 `_ - - Parameters - ---------- - mass_operator : WeightedMassOperator - The weighted mass operator for which the approximate inverse is needed. - - apply_bc : bool - Whether to include boundary operators. - """ - - def __init__(self, mass_operator: WeightedMassOperator, apply_bc: bool = True): - super().__init__(mass_operator, apply_bc=apply_bc, dim_reduce=None, diagonal_scaling=True) - - def transpose(self, conjugate: bool = False) -> MassMatrixDiagonalPreconditioner: - """ - Returns the transposed operator. - """ - return MassMatrixDiagonalPreconditioner(self._mass_operator.transpose(), self._apply_bc) - - class StiffnessPreconditioner(KroneckerPreconditioner): r""" Preconditioner for the stabilized stiffness operators @@ -587,7 +560,8 @@ def _build_approximation(self, bc): if d in stiff_dirs[c]: assert basis == "B", "Only B-spline directions are differentiated." M_D, _ = _dense_mass_1d(self._mass_operator, "M", d, shapes[d]) - D_d = _difference_matrix_1d(M_d.shape[0], M_D.shape[0]) + D_d = _dense_derivative_1d(self._d, stiff_dirs[c][d], c) + assert D_d.shape == (M_D.shape[0], M_d.shape[0]) S_d = mid_means[stiff_dirs[c][d]] * (D_d.T @ M_D @ D_d) if bc is not None and basis == "B": @@ -1149,17 +1123,18 @@ def _dense_mass_1d( return np.asarray(xp.to_numpy(M.toarray()), dtype=float), domain_decomposition -def _difference_matrix_1d(n_B: int, n_M: int) -> np.ndarray: +def _dense_derivative_1d(derivative: LinearOperator, r: int, c: int) -> np.ndarray: """ - 1d derivative matrix from B-spline to M-spline coefficients, :math:`(\\mathbb d c)_r = c_{r+1} - c_r` - (periodic if ``n_M == n_B``, else ``n_M == n_B - 1``). + Global 1d derivative matrix (B- to M-spline coefficients) of the block ``(r, c)`` of a + tensor-product derivative, as dense host array, taken from the 1d factor of the block itself + (:meth:`DirectionalDerivativeOperator.tokronstencil`), so that it matches the 3d operator + (e.g. it vanishes in a periodic direction with a single element). """ - assert n_M in (n_B, n_B - 1) - D = np.zeros((n_M, n_B)) - for r in range(n_M): - D[r, r] = -1.0 - D[r, (r + 1) % n_B] = 1.0 - return D + op = derivative[r, c] if isinstance(derivative, BlockLinearOperator) else derivative + assert isinstance(op, DirectionalDerivativeOperator) and not op.transposed + D = op.tokronstencil().mats[op.diffdir].toarray() + # the sign of a negated block is already included in its 1d factor + return np.asarray(xp.to_numpy(D), dtype=float) def _apply_bc_dense(A: np.ndarray, bc_d: tuple[bool, bool]) -> np.ndarray: diff --git a/src/struphy/feec/tests/test_preconditioner_transpose.py b/src/struphy/feec/tests/test_preconditioner_transpose.py index f0b49eaf6..5796df6bd 100644 --- a/src/struphy/feec/tests/test_preconditioner_transpose.py +++ b/src/struphy/feec/tests/test_preconditioner_transpose.py @@ -1,15 +1,15 @@ import pytest -@pytest.mark.parametrize("dim_reduce", [0, 2]) +@pytest.mark.parametrize("dim_reduce", [0, 2, None]) @pytest.mark.parametrize("weight_reduction", ["midpoint", "average"]) @pytest.mark.parametrize("diagonal_scaling", [False, True]) def test_mass_preconditioner_transpose(dim_reduce, weight_reduction, diagonal_scaling): - """The transposes of the mass-matrix preconditioners keep their class and options.""" + """The transpose of the mass-matrix preconditioner keeps its class and options.""" from struphy import domains from struphy.feec.mass import WeightedMassOperators - from struphy.feec.preconditioner import MassMatrixDiagonalPreconditioner, MassMatrixPreconditioner + from struphy.feec.preconditioner import MassMatrixPreconditioner from struphy.feec.psydac_derham import Derham from struphy.io.options import DerhamOptions from struphy.topology.grids import TensorProductGrid @@ -32,18 +32,14 @@ def test_mass_preconditioner_transpose(dim_reduce, weight_reduction, diagonal_sc assert pc_T.weight_reduction == weight_reduction assert pc_T.diagonal_scaling == diagonal_scaling - pc_diag_T = MassMatrixDiagonalPreconditioner(mass_ops.M1, apply_bc=False).transpose() - assert type(pc_diag_T) is MassMatrixDiagonalPreconditioner - assert pc_diag_T._apply_bc is False - -def test_mass_diagonal_preconditioner_assembles_1d_matrices_only(monkeypatch): - """MassMatrixDiagonalPreconditioner assembles only 1d mass matrices (no 3d logical mass matrix); +def test_mass_preconditioner_unweighted_assembles_1d_matrices_only(monkeypatch): + """MassMatrixPreconditioner with dim_reduce=None assembles only 1d mass matrices (no 3d logical mass matrix); the diagonal of the logical mass matrix is the one of its Kronecker approximation.""" from struphy import domains from struphy.feec.mass import WeightedMassOperator, WeightedMassOperators - from struphy.feec.preconditioner import MassMatrixDiagonalPreconditioner, MassMatrixPreconditioner + from struphy.feec.preconditioner import MassMatrixPreconditioner from struphy.feec.psydac_derham import Derham from struphy.io.options import DerhamOptions from struphy.topology.grids import TensorProductGrid @@ -61,14 +57,13 @@ def recording_init(self, derham, V, W, *args, **kwargs): init(self, derham, V, W, *args, **kwargs) monkeypatch.setattr(WeightedMassOperator, "__init__", recording_init) - pc = MassMatrixDiagonalPreconditioner(M1) + pc = MassMatrixPreconditioner(M1, dim_reduce=None) - assert isinstance(pc, MassMatrixPreconditioner) assert pc.dim_reduce is None and pc.diagonal_scaling assert created and all(ldim == 1 for ldim in created) -def test_mass_diagonal_preconditioner_zero_diagonal_block(): +def test_mass_preconditioner_unweighted_zero_diagonal_block(): """A mass operator with a zero (None) diagonal block, e.g. after reassembly with a zero weight: construction and update_mass_operator must work; the scaling of that block is 1.""" @@ -76,7 +71,7 @@ def test_mass_diagonal_preconditioner_zero_diagonal_block(): from struphy import domains from struphy.feec.mass import WeightedMassOperator, WeightedMassOperators - from struphy.feec.preconditioner import MassMatrixDiagonalPreconditioner + from struphy.feec.preconditioner import MassMatrixPreconditioner from struphy.feec.psydac_derham import Derham from struphy.feec.utilities import create_equal_random_arrays from struphy.io.options import DerhamOptions @@ -94,7 +89,7 @@ def one(e1, e2, e3): M.assemble() assert M.matrix[2, 2] is None - pc = MassMatrixDiagonalPreconditioner(mass_ops.M1) + pc = MassMatrixPreconditioner(mass_ops.M1, dim_reduce=None) pc.update_mass_operator(M) scaling = pc._scaling[0] assert xp.allclose(scaling[2, 2]._data, 1.0) @@ -103,7 +98,7 @@ def one(e1, e2, e3): assert all(xp.all(xp.isfinite(b.toarray())) for b in pc.dot(x).blocks) -def test_mass_diagonal_preconditioner_small_periodic_direction(): +def test_mass_preconditioner_unweighted_small_periodic_direction(): """Periodic direction with fewer points than the stencil width (1 element, degree 1): the main diagonal of the 1d factors must be at offset 0, otherwise the diagonal scaling vanishes.""" @@ -112,7 +107,7 @@ def test_mass_diagonal_preconditioner_small_periodic_direction(): from struphy import domains from struphy.feec.mass import WeightedMassOperators - from struphy.feec.preconditioner import MassMatrixDiagonalPreconditioner, _local_diagonal + from struphy.feec.preconditioner import MassMatrixPreconditioner, _local_diagonal from struphy.feec.psydac_derham import Derham from struphy.feec.utilities import create_equal_random_arrays from struphy.io.options import DerhamOptions @@ -121,7 +116,7 @@ def test_mass_diagonal_preconditioner_small_periodic_direction(): derham = Derham(grid=TensorProductGrid(num_elements=(4, 4, 1)), options=DerhamOptions(degree=(2, 2, 1))) mass_ops = WeightedMassOperators(derham=derham, domain=domains.Cuboid()) - pc = MassMatrixDiagonalPreconditioner(mass_ops.Mv) + pc = MassMatrixPreconditioner(mass_ops.Mv, dim_reduce=None) assert all(float(xp.min(d)) > 0.0 for d in _local_diagonal(pc.matrix)) _, x = create_equal_random_arrays(derham.Vvfem, seed=1) diff --git a/src/struphy/feec/tests/test_stiffness_preconditioner.py b/src/struphy/feec/tests/test_stiffness_preconditioner.py index 6305551ca..5364dffb9 100644 --- a/src/struphy/feec/tests/test_stiffness_preconditioner.py +++ b/src/struphy/feec/tests/test_stiffness_preconditioner.py @@ -19,7 +19,8 @@ def _max_abs(arrays, comm): @pytest.mark.parametrize("derivative", ["grad", "curl", "div"]) @pytest.mark.parametrize("sigma", [0.0, 0.7]) -def test_stiffness_approximation_unit_cube(derivative, sigma): +@pytest.mark.parametrize("num_elements", [[6, 5, 4], [6, 5, 1]]) +def test_stiffness_approximation_unit_cube(derivative, sigma, num_elements): """On the unit cube (unit weights), the diagonal blocks of the Kronecker approximation equal those of the stiffness operator, the solver inverts the approximation, and probing gives its exact diagonal.""" @@ -35,7 +36,9 @@ def test_stiffness_approximation_unit_cube(derivative, sigma): comm = MPI.COMM_WORLD derham = Derham( - TensorProductGrid(num_elements=[6, 5, 4]), DerhamOptions(degree=[2, 3, 2], bcs=(None, None, None)), comm=comm + TensorProductGrid(num_elements=num_elements), + DerhamOptions(degree=[2, 3, 2 if num_elements[2] > 1 else 1], bcs=(None, None, None)), + comm=comm, ) mass_ops = WeightedMassOperators(derham, domains.Cuboid()) @@ -59,7 +62,8 @@ def test_stiffness_approximation_unit_cube(derivative, sigma): b._data[:] = 0.0 r = A.dot(xc) - P.matrix.dot(xc) ref = A.dot(xc) - assert _max_abs([_blocks(r)[c].toarray()], comm) < 1e-12 * _max_abs([_blocks(ref)[c].toarray()], comm) + # <= since a block vanishes for a derivative in a direction with a single periodic element + assert _max_abs([_blocks(r)[c].toarray()], comm) <= 1e-12 * _max_abs([_blocks(ref)[c].toarray()], comm) if sigma > 0: r = P.solver.dot(P.matrix.dot(x)) - x diff --git a/src/struphy/feec/variational_utilities.py b/src/struphy/feec/variational_utilities.py index c90d1a03e..c8dc0fa72 100644 --- a/src/struphy/feec/variational_utilities.py +++ b/src/struphy/feec/variational_utilities.py @@ -1376,7 +1376,7 @@ def update_weight(self, coeffs): def _create_inv( self, type="pcg", - pc_type="MassMatrixDiagonalPreconditioner", + pc_type="MassMatrixPreconditioner", tol=1e-16, maxiter=500, verbose=False, @@ -1389,7 +1389,8 @@ def _create_inv( preconditioner, pc_type, ) - self._pc = pc_class(self.massop) + # dim_reduce=None: the Kronecker factors do not depend on the weights, cheap to update + self._pc = pc_class(self.massop, dim_reduce=None) self._inv = inverse( self.massop, diff --git a/src/struphy/io/options.py b/src/struphy/io/options.py index 289404c52..0035d0ba6 100644 --- a/src/struphy/io/options.py +++ b/src/struphy/io/options.py @@ -75,8 +75,13 @@ class LiteralOptions: # solvers OptsSymmSolver = Literal["pcg", "cg"] OptsGenSolver = Literal["pbicgstab", "bicgstab", "gmres"] - OptsMassPrecond = Literal["MassMatrixPreconditioner", "MassMatrixDiagonalPreconditioner", None] - OptsDiffusionPrecond = Literal["MultiGrid", "MassMatrixPreconditioner", "MassMatrixDiagonalPreconditioner", None] + OptsMassPrecond = Literal["MassMatrixPreconditioner", None] + OptsDiffusionPrecond = Literal[ + "MultiGrid", + "StiffnessPreconditioner", + "MassMatrixPreconditioner", + None, + ] OptsSaddlePointSolver = Literal["uzawa"] OptsDirectSolver = Literal["SparseSolver", "ScipySparse", "InexactNPInverse", "DirectNPInverse"] OptsNonlinearSolver = Literal["Picard", "Newton"] diff --git a/src/struphy/propagators/implicit_diffusion.py b/src/struphy/propagators/implicit_diffusion.py index 2af8fcd03..8d0c6ce1a 100644 --- a/src/struphy/propagators/implicit_diffusion.py +++ b/src/struphy/propagators/implicit_diffusion.py @@ -9,6 +9,7 @@ from line_profiler import profile from scope_profiler import ProfileManager +from struphy.feec import preconditioner from struphy.feec.mass import L2Projector, WeightedMassOperator from struphy.io.options import LiteralOptions, OptionsBase from struphy.linear_algebra.multigrid.preconditioner import MultiGridOptions, MultiGridPreconditioner @@ -183,11 +184,26 @@ class Options(OptionsBase): Name of the symmetric iterative solver passed to :func:`psydac.linalg.solvers.inverse`. - precond : LiteralOptions.OptsDiffusionPrecond, default="MassMatrixPreconditioner" - Name of the preconditioner configuration. - ``"MultiGrid"`` uses a geometric multigrid V-cycle - (:class:`~struphy.linear_algebra.multigrid.preconditioner.MultiGridPreconditioner`, - requires ``solver="pcg"``). The other values currently result in ``pc=None``. + precond : LiteralOptions.OptsDiffusionPrecond, default=None + Name of the preconditioner: + + - ``None``: no preconditioner. + - ``"MultiGrid"``: geometric multigrid V-cycle + (:class:`~struphy.linear_algebra.multigrid.preconditioner.MultiGridPreconditioner`, + requires ``solver="pcg"``). + - ``"MassMatrixPreconditioner"``: Kronecker + approximation of the inverse of the 0-form mass matrix ``M0`` + (see :mod:`struphy.feec.preconditioner`). + - ``"StiffnessPreconditioner"``: Kronecker approximation of the inverse of + ``grad.T @ M1 @ grad + sigma * M0`` + (:class:`~struphy.feec.preconditioner.StiffnessPreconditioner`), with + ``sigma = sigma_1`` (divided by dt if ``divide_by_dt=True``); any ``stab_mat`` is + approximated by ``M0`` (this keeps periodic or Neumann problems regular). + + precond_params : dict, default=None + Keyword arguments passed to the constructor of the mass-matrix or stiffness + preconditioner, e.g. ``{"weights": "unit"}`` for ``"StiffnessPreconditioner"`` + or ``{"dim_reduce": None}`` for ``"MassMatrixPreconditioner"``. multigrid : MultiGridOptions, default=None Options of the multigrid preconditioner (if ``precond="MultiGrid"``). @@ -223,7 +239,8 @@ class Options(OptionsBase): diffusion_mat: OptsDiffusionMat = "M1" x0: StencilVector = None solver: LiteralOptions.OptsSymmSolver = "pcg" - precond: LiteralOptions.OptsDiffusionPrecond = "MassMatrixPreconditioner" + precond: LiteralOptions.OptsDiffusionPrecond = None + precond_params: dict = None multigrid: MultiGridOptions = None solver_params: SolverParameters = None filter_params: dict[PICVariable | SPHVariable, FilterParameters] = None @@ -234,10 +251,12 @@ def __post_init__(self): check_option(self.diffusion_mat, self.OptsDiffusionMat) check_option(self.solver, LiteralOptions.OptsSymmSolver) check_option(self.precond, LiteralOptions.OptsDiffusionPrecond) - if self.precond == "MultiGrid": - assert self.solver == "pcg", "precond='MultiGrid' requires solver='pcg'." + if self.precond is not None: + assert self.solver == "pcg", f"precond={self.precond!r} requires solver='pcg'." # defaults + if self.precond_params is None: + self.precond_params = {} if self.solver_params is None: self.solver_params = SolverParameters() if self.multigrid is None: @@ -357,10 +376,15 @@ def verify_rhs(rho) -> StencilVector | FEECVariable | AccumulatorVector: self._diffusion_op = self.derham.grad.T @ diffusion_mat @ self.derham.grad # preconditioner and solver for Ax=b + # the multigrid and stiffness preconditioners depend on sigma_1, they are updated in __call__ + # if it changes (e.g. with dt); _pc_sig_1 is the sigma_1 they were built with self._mg = None - if self.options.precond == "MultiGrid": - # the operator is updated in __call__ if sigma_1 changes (e.g. with dt) - self._mg_sig_1 = self._sigma_1 + self._stiffness_pc = None + self._pc_sig_1 = self._sigma_1 + precond = self.options.precond + if precond is None: + pc = None + elif precond == "MultiGrid": self._mg = MultiGridPreconditioner( self._sigma_1 * stab_mat + self._diffusion_op, self.derham, @@ -369,15 +393,20 @@ def verify_rhs(rho) -> StencilVector | FEECVariable | AccumulatorVector: mass_ops=self.mass_ops, ) pc = self._mg + elif precond == "StiffnessPreconditioner": + self._stiffness_pc = self._build_stiffness_pc(self._sigma_1) + pc = self._stiffness_pc else: - # TODO: mass-matrix preconditioners are not effective for this operator - pc = None + pc_class = getattr(preconditioner, precond) + pc = pc_class(self.mass_ops.M0, **self.options.precond_params) # solver just with A_2, but will be set during call with dt + # (only "pcg" takes a preconditioner) + pc_kwargs = {"pc": pc} if self.options.solver == "pcg" else {} self._solver = inverse( self._diffusion_op, self.options.solver, - pc=pc, + **pc_kwargs, x0=self.x0, tol=self.options.solver_params.tol, maxiter=self.options.solver_params.maxiter, @@ -391,6 +420,18 @@ def verify_rhs(rho) -> StencilVector | FEECVariable | AccumulatorVector: self._rhs2 = phi.space.zeros() self._tmp_src = phi.space.zeros() + def _build_stiffness_pc(self, sig_1: float) -> preconditioner.StiffnessPreconditioner: + """StiffnessPreconditioner for ``sig_1 * stab_mat + grad.T @ M1 @ grad``; the stabilization + is approximated by ``sig_1 * M0`` for any ``stab_mat`` (with ``sigma = 0``, the pseudo-inverse + would drop the constants, on which the operator is not singular for periodic or Neumann + boundary conditions).""" + return preconditioner.StiffnessPreconditioner( + self.mass_ops, + "grad", + sigma=sig_1, + **self.options.precond_params, + ) + @property def sources(self) -> list[StencilVector | FEECVariable | AccumulatorVector]: """ @@ -481,9 +522,13 @@ def __call__(self, dt): # compute lhs self._solver.linop = sig_1 * self._stab_mat + self._diffusion_op - if self._mg is not None and sig_1 != self._mg_sig_1: - self._mg.update(self._solver.linop) - self._mg_sig_1 = sig_1 + if sig_1 != self._pc_sig_1: + if self._mg is not None: + self._mg.update(self._solver.linop) + if self._stiffness_pc is not None: + self._stiffness_pc = self._build_stiffness_pc(sig_1) + self._solver.set_options(pc=self._stiffness_pc) + self._pc_sig_1 = sig_1 # solve with ProfileManager.profile_region(self._solve_region, functions=[self._solver.solve]): diff --git a/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py b/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py index 7cacf5434..fbf6dab48 100644 --- a/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py +++ b/src/struphy/propagators/poisson_adiabatic_gyrokinetic.py @@ -131,10 +131,15 @@ class Options: Name of the symmetric iterative solver passed to :func:`psydac.linalg.solvers.inverse`. - precond : LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner" - Name of the preconditioner configuration. - Currently this class inherits the same behavior as - :class:`ImplicitDiffusion`, where ``pc=None`` is used internally. + precond : LiteralOptions.OptsDiffusionPrecond, default=None + Name of the preconditioner, see :class:`ImplicitDiffusion`. Requires ``solver="pcg"``. + ``"MultiGrid"`` is not supported: multigrid would be built for the stabilization without + the average, which is applied after ``ImplicitDiffusion.allocate()``. + + precond_params : dict, default=None + Keyword arguments passed to the constructor of the mass-matrix or stiffness + preconditioner, see :class:`ImplicitDiffusion`. + solver_params : SolverParameters, default=None Iterative-solver controls (for example ``tol``, ``maxiter``, @@ -162,7 +167,8 @@ class Options: which_geometry: OptsGeometry = "cylindrical" x0: StencilVector = None solver: LiteralOptions.OptsSymmSolver = "pcg" - precond: LiteralOptions.OptsMassPrecond = "MassMatrixPreconditioner" + precond: LiteralOptions.OptsDiffusionPrecond = None + precond_params: dict = None solver_params: SolverParameters = None filter_params: dict[PICVariable | SPHVariable, FilterParameters] = None @@ -170,9 +176,17 @@ def __post_init__(self): # checks check_option(self.stab_mat, self.OptsStabMat) check_option(self.solver, LiteralOptions.OptsSymmSolver) - check_option(self.precond, LiteralOptions.OptsMassPrecond) + check_option(self.precond, LiteralOptions.OptsDiffusionPrecond) + if self.precond is not None: + assert self.solver == "pcg", f"precond={self.precond!r} requires solver='pcg'." + if self.precond == "MultiGrid": + raise ValueError( + "precond='MultiGrid' is not supported (it would not see the average in the stabilization)." + ) # defaults + if self.precond_params is None: + self.precond_params = {} if self.solver_params is None: self.solver_params = SolverParameters() diff --git a/src/struphy/propagators/poisson_solve.py b/src/struphy/propagators/poisson_solve.py index 9e71cb73c..9dd8998cc 100644 --- a/src/struphy/propagators/poisson_solve.py +++ b/src/struphy/propagators/poisson_solve.py @@ -45,7 +45,7 @@ class Options(OptionsBase): Internally mapped to ``sigma_1 = stab_eps`` in the parent :class:`ImplicitDiffusion` formulation. - stab_mat : {"M0", "M0ad", "Id"}, default="Id" + stab_mat : {"M0", "M0ad", "Id"}, default="M0" Stabilization matrix multiplied by ``stab_eps``. - ``"M0"``: standard weighted 0-form mass operator. @@ -67,9 +67,15 @@ class Options(OptionsBase): Name of the symmetric iterative solver passed to :func:`psydac.linalg.solvers.inverse`. - precond : LiteralOptions.OptsDiffusionPrecond, default="MassMatrixPreconditioner" - Name of the preconditioner configuration, see :class:`ImplicitDiffusion` - (``"MultiGrid"`` for geometric multigrid). + precond : LiteralOptions.OptsDiffusionPrecond, default=None + Name of the preconditioner, see :class:`ImplicitDiffusion` + (``"MultiGrid"`` for geometric multigrid, ``"StiffnessPreconditioner"`` for the + Kronecker approximation of the Laplacian, ``"MassMatrixPreconditioner"`` for the + Kronecker approximation of ``M0``). Requires ``solver="pcg"``. + + precond_params : dict, default=None + Keyword arguments passed to the constructor of the mass-matrix or stiffness + preconditioner, see :class:`ImplicitDiffusion`. multigrid : MultiGridOptions, default=None Options of the multigrid preconditioner (if ``precond="MultiGrid"``). @@ -97,11 +103,12 @@ class Options(OptionsBase): OptsDiffusionMat = Literal["M1", "M1perp", "M1para", "M1gyro"] # propagator options stab_eps: float = 0.0 - stab_mat: OptsStabMat = "Id" + stab_mat: OptsStabMat = "M0" diffusion_mat: OptsDiffusionMat = "M1" x0: StencilVector = None solver: LiteralOptions.OptsSymmSolver = "pcg" - precond: LiteralOptions.OptsDiffusionPrecond = "MassMatrixPreconditioner" + precond: LiteralOptions.OptsDiffusionPrecond = None + precond_params: dict = None multigrid: MultiGridOptions = None solver_params: SolverParameters = None filter_params: dict[PICVariable | SPHVariable, FilterParameters] = None @@ -112,10 +119,12 @@ def __post_init__(self): check_option(self.diffusion_mat, self.OptsDiffusionMat) check_option(self.solver, LiteralOptions.OptsSymmSolver) check_option(self.precond, LiteralOptions.OptsDiffusionPrecond) - if self.precond == "MultiGrid": - assert self.solver == "pcg", "precond='MultiGrid' requires solver='pcg'." + if self.precond is not None: + assert self.solver == "pcg", f"precond={self.precond!r} requires solver='pcg'." # defaults + if self.precond_params is None: + self.precond_params = {} if self.solver_params is None: self.solver_params = SolverParameters() if self.multigrid is None: diff --git a/src/struphy/propagators/shear_alfven_current_coupling_5d.py b/src/struphy/propagators/shear_alfven_current_coupling_5d.py index 0ee45bcb4..117cf94f3 100644 --- a/src/struphy/propagators/shear_alfven_current_coupling_5d.py +++ b/src/struphy/propagators/shear_alfven_current_coupling_5d.py @@ -122,7 +122,7 @@ class Options(OptionsBase): Time stepping algorithm. solver : LiteralOptions.OptsSymmSolver, default="pcg" Symmetric iterative solver. - precond : LiteralOptions.OptsMassPrecond, default="MassMatrixDiagonalPreconditioner" + precond : LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner" Preconditioner for mass-matrix block. solver_params : SolverParameters, default=None Solver controls. @@ -141,7 +141,7 @@ class Options(OptionsBase): u_space: LiteralOptions.OptsVecSpace = "Hdiv" algo: OptsAlgo = "implicit" solver: LiteralOptions.OptsSymmSolver = "pcg" - precond: LiteralOptions.OptsMassPrecond = "MassMatrixDiagonalPreconditioner" + precond: LiteralOptions.OptsMassPrecond = "MassMatrixPreconditioner" solver_params: SolverParameters = None filter_params: FilterParameters = None butcher: ButcherTableau = None diff --git a/src/struphy/propagators/tests/test_gyrokinetic_poisson.py b/src/struphy/propagators/tests/test_gyrokinetic_poisson.py index 08d1f5f4a..3e7d7414f 100644 --- a/src/struphy/propagators/tests/test_gyrokinetic_poisson.py +++ b/src/struphy/propagators/tests/test_gyrokinetic_poisson.py @@ -170,7 +170,7 @@ def rho_pulled(e1, e2, e3): divide_by_dt=True, diffusion_mat="M1perp", solver="pcg", - precond="MassMatrixPreconditioner", + precond=None, solver_params=solver_params, ) @@ -377,7 +377,7 @@ def rho2_pulled(e1, e2, e3): divide_by_dt=True, diffusion_mat="M1perp", solver="pcg", - precond="MassMatrixPreconditioner", + precond=None, solver_params=solver_params, ) @@ -396,7 +396,7 @@ def rho2_pulled(e1, e2, e3): divide_by_dt=True, diffusion_mat="M1perp", solver="pcg", - precond="MassMatrixPreconditioner", + precond=None, solver_params=solver_params, ) @@ -519,7 +519,7 @@ def rho(e1, e2, e3): stab_eps=1e-8, diffusion_mat="M1", solver="pcg", - precond="MassMatrixPreconditioner", + precond=None, solver_params=solver_params, ) @@ -538,7 +538,7 @@ def rho(e1, e2, e3): stab_eps=1e-8, diffusion_mat="M1perp", solver="pcg", - precond="MassMatrixPreconditioner", + precond=None, solver_params=solver_params, ) @@ -694,7 +694,7 @@ def rho_physical(x, y, z): stab_eps=1e-8, diffusion_mat="M1perp", solver="pcg", - precond="MassMatrixPreconditioner", + precond=None, solver_params=solver_params, ) @@ -748,7 +748,7 @@ def rho_physical(x, y, z): stab_eps=1e-8, diffusion_mat="M1", solver="pcg", - precond="MassMatrixPreconditioner", + precond=None, solver_params=solver_params, ) poisson_solver_2p5d.allocate() diff --git a/src/struphy/propagators/tests/test_poisson.py b/src/struphy/propagators/tests/test_poisson.py index 84524ee77..9c15ee437 100644 --- a/src/struphy/propagators/tests/test_poisson.py +++ b/src/struphy/propagators/tests/test_poisson.py @@ -203,7 +203,7 @@ def rho_pulled(e1, e2, e3): # sigma_2=0.0, # sigma_3=1.0, solver="pcg", - precond="MassMatrixPreconditioner", + precond=None, solver_params=solver_params, ) @@ -368,7 +368,7 @@ def test_poisson_accum_1d(mapping, do_plot=False): # sigma_2=0.0, # sigma_3=1.0, solver="pcg", - precond="MassMatrixPreconditioner", + precond=None, solver_params=solver_params, ) @@ -470,7 +470,7 @@ def test_poisson_accum_full_f_background_1d(): poisson_solver.options = poisson_solver.Options( stab_eps=1e-6, solver="pcg", - precond="MassMatrixPreconditioner", + precond=None, solver_params=SolverParameters(tol=1.0e-12, maxiter=3000), ) poisson_solver.allocate() @@ -695,7 +695,7 @@ def rho2_pulled(e1, e2, e3): # sigma_2=0.0, # sigma_3=1.0, solver="pcg", - precond="MassMatrixPreconditioner", + precond=None, solver_params=solver_params, ) @@ -723,7 +723,7 @@ def rho2_pulled(e1, e2, e3): # sigma_2=0.0, # sigma_3=1.0, solver="pcg", - precond="MassMatrixPreconditioner", + precond=None, solver_params=solver_params, ) @@ -808,7 +808,7 @@ def rho(e1, e2, e3): phis = [] infos = [] - for precond in ["MassMatrixPreconditioner", "MultiGrid"]: + for precond in [None, "MultiGrid"]: phi = FEECVariable(space="H1") phi.allocate(derham=derham, domain=domain) solver = PoissonSolve(rho=rho) @@ -880,6 +880,105 @@ def test_implicit_diffusion_multigrid_dt(): assert prop._solver._info["niter"] <= 15 +@pytest.mark.parametrize("stab_mat", ["M0", "Id"]) +@pytest.mark.parametrize("bc_type", ["periodic", "dirichlet", "neumann"]) +def test_poisson_2d_kronecker_preconds(bc_type, stab_mat): + """PoissonSolve with the Kronecker preconditioners agrees with the unpreconditioned solve (pcg and cg); + the stiffness preconditioner needs fewer iterations than no preconditioner (the mass-matrix + preconditioner does not approximate the Laplacian and is only checked for the solution).""" + domain = domains.Colella(Lx=4.0, Ly=2.0, alpha=0.1, Lz=1.0) + bcs = { + "periodic": (None, None, None), + "dirichlet": (("dirichlet", "dirichlet"), None, None), + "neumann": (("free", "free"), None, None), + }[bc_type] + derham = Derham(TensorProductGrid(num_elements=[32, 32, 1]), DerhamOptions(degree=[2, 2, 1], bcs=bcs), comm=comm) + mass_ops = WeightedMassOperators(derham, domain) + Propagator.derham = derham + Propagator.domain = domain + Propagator.mass_ops = mass_ops + + def rho(e1, e2, e3): + return xp.cos(2 * xp.pi * e1) * xp.sin(2 * xp.pi * e2) + 0.3 * xp.sin(4 * xp.pi * e2) + + phis = [] + infos = [] + for solver_name, precond in [ + ("pcg", None), + ("pcg", "MassMatrixPreconditioner"), + ("pcg", "StiffnessPreconditioner"), + ("cg", None), + ]: + phi = FEECVariable(space="H1") + phi.allocate(derham=derham, domain=domain) + solver = PoissonSolve(rho=rho) + solver.variables.phi = phi + solver.options = solver.Options( + stab_eps=1e-8, + stab_mat=stab_mat, + solver=solver_name, + precond=precond, + solver_params=SolverParameters(tol=1e-11, maxiter=3000, recycle=False), + ) + solver.allocate() + solver(1.0) + phis.append(phi.spline.vector.toarray()) + infos.append(solver._solver._info) + + # global coefficient arrays (toarray only fills the local part) + if comm.Get_size() > 1: + phis = [comm.allreduce(p, op=MPI.SUM) for p in phis] + if bc_type != "dirichlet": + # solutions are defined up to a constant (the stabilization is tiny) + phis = [p - xp.mean(p) for p in phis] + for p in phis[1:]: + assert xp.max(xp.abs(p - phis[0])) < 1e-6 * xp.max(xp.abs(phis[0])) + assert all(info["success"] for info in infos) + assert infos[2]["niter"] < infos[0]["niter"] + + +def test_implicit_diffusion_stiffness_dt(): + """With divide_by_dt, the stiffness preconditioner follows changes of dt.""" + from struphy.propagators.implicit_diffusion import ImplicitDiffusion + + domain = domains.Cuboid(l1=0.0, r1=2.0, l2=0.0, r2=1.0, l3=0.0, r3=1.0) + derham = Derham( + TensorProductGrid(num_elements=[32, 16, 1]), + DerhamOptions(degree=[2, 2, 1], bcs=(("dirichlet", "dirichlet"), None, None)), + comm=comm, + ) + mass_ops = WeightedMassOperators(derham, domain) + Propagator.derham = derham + Propagator.domain = domain + Propagator.mass_ops = mass_ops + + phi = FEECVariable(space="H1") + phi.allocate(derham=derham, domain=domain) + phi.spline.vector = derham.P0(lambda e1, e2, e3: xp.sin(xp.pi * e1) * xp.cos(2 * xp.pi * e2)) + + prop = ImplicitDiffusion() + prop.variables.phi = phi + prop.options = prop.Options( + sigma_1=1.0, + sigma_2=1.0, + sigma_3=0.0, + divide_by_dt=True, + precond="StiffnessPreconditioner", + solver_params=SolverParameters(tol=1e-12, maxiter=100, recycle=False), + ) + prop.allocate() + + for dt in [0.1, 0.1, 0.01]: + rhs = (1.0 / dt) * mass_ops.M0.dot(phi.spline.vector) + prop(dt) + assert prop._stiffness_pc.sigma == 1.0 / dt + A = (1.0 / dt) * mass_ops.M0 + derham.grad.T @ mass_ops.M1 @ derham.grad + r = rhs - A.dot(phi.spline.vector) + assert xp.sqrt(r.inner(r)) < 1e-9 * xp.sqrt(rhs.inner(rhs)) + # exact inverse on the Cuboid (up to round-off) + assert prop._solver._info["niter"] <= 3 + + if __name__ == "__main__": # direction = 0 # bc_type = "dirichlet" diff --git a/src/struphy/propagators/variational_density_evolve.py b/src/struphy/propagators/variational_density_evolve.py index f63d80a8e..2b54296f7 100644 --- a/src/struphy/propagators/variational_density_evolve.py +++ b/src/struphy/propagators/variational_density_evolve.py @@ -11,7 +11,7 @@ from maybempi import MPI from struphy.feec import preconditioner -from struphy.feec.preconditioner import MassMatrixDiagonalPreconditioner +from struphy.feec.preconditioner import MassMatrixPreconditioner from struphy.feec.variational_utilities import ( InternalEnergyEvaluator, KineticEnergyEvaluator, @@ -202,7 +202,7 @@ def allocate(self): self._info = self.options.nonlin_solver.info and (MPI.COMM_WORLD.Get_rank() == 0) self._Mrho = self.mass_ops.WMMnew - pc = MassMatrixDiagonalPreconditioner(self._Mrho) + pc = MassMatrixPreconditioner(self._Mrho, dim_reduce=None) self._Mrho_inv = inverse( self._Mrho, self.options.solver, @@ -408,8 +408,9 @@ def _initialize_projectors_and_mass(self): self.divPirhoT = self.divPirho.T # Inverse mass matrix needed to compute the error - self.pc_Mv = preconditioner.MassMatrixDiagonalPreconditioner( + self.pc_Mv = preconditioner.MassMatrixPreconditioner( self.mass_ops.Mv, + dim_reduce=None, ) self._inv_Mv = inverse( self.mass_ops.Mv, @@ -532,9 +533,7 @@ def _update_weighted_MM(self, rho): self._Mrho.assemble() logger.debug(f"In VariationalDensityEvolve: {self._Mrho_inv._options.get('pc') = }") - if hasattr(self, "_Mrho_inv") and isinstance( - self._Mrho_inv._options.get("pc"), MassMatrixDiagonalPreconditioner - ): + if hasattr(self, "_Mrho_inv") and isinstance(self._Mrho_inv._options.get("pc"), MassMatrixPreconditioner): self._Mrho_inv._options["pc"].update_mass_operator(self._Mrho) def _update_linear_form_dl_drho(self, rhon, rhon1, un, un1, sn): diff --git a/src/struphy/propagators/variational_entropy_evolve.py b/src/struphy/propagators/variational_entropy_evolve.py index 4e7cbdce4..d59b8308a 100644 --- a/src/struphy/propagators/variational_entropy_evolve.py +++ b/src/struphy/propagators/variational_entropy_evolve.py @@ -11,7 +11,7 @@ from maybempi import MPI from struphy.feec import preconditioner -from struphy.feec.preconditioner import MassMatrixDiagonalPreconditioner +from struphy.feec.preconditioner import MassMatrixPreconditioner from struphy.feec.variational_utilities import InternalEnergyEvaluator from struphy.io.options import LiteralOptions, OptionsBase from struphy.linear_algebra.solver import NonlinearSolverParameters, SolverParameters @@ -182,7 +182,7 @@ def allocate(self): # assembly of WMMnew happens in VariationalDensityEvolve self._Mrho = self.mass_ops.WMMnew - pc = MassMatrixDiagonalPreconditioner(self._Mrho) + pc = MassMatrixPreconditioner(self._Mrho, dim_reduce=None) self._Mrho_inv = inverse( self._Mrho, "pcg", @@ -320,8 +320,9 @@ def _initialize_projectors_and_mass(self): self.divPisT = self.divPis.T # Inverse mass matrix needed to compute the error - self.pc_Mv = preconditioner.MassMatrixDiagonalPreconditioner( + self.pc_Mv = preconditioner.MassMatrixPreconditioner( self.mass_ops.Mv, + dim_reduce=None, ) self._inv_Mv = inverse( self.mass_ops.Mv, diff --git a/src/struphy/propagators/variational_mag_field_evolve.py b/src/struphy/propagators/variational_mag_field_evolve.py index 1a7c98604..1e501ad65 100644 --- a/src/struphy/propagators/variational_mag_field_evolve.py +++ b/src/struphy/propagators/variational_mag_field_evolve.py @@ -10,7 +10,7 @@ from maybempi import MPI from struphy.feec import preconditioner -from struphy.feec.preconditioner import MassMatrixDiagonalPreconditioner +from struphy.feec.preconditioner import MassMatrixPreconditioner from struphy.feec.variational_utilities import Hdiv0_transport_operator from struphy.io.options import LiteralOptions, OptionsBase from struphy.linear_algebra.schur_solver import SchurSolverFull @@ -165,7 +165,7 @@ def allocate(self): # assembly of WMMnew happens in VariationalDensityEvolve self._Mrho = self.mass_ops.WMMnew - pc = MassMatrixDiagonalPreconditioner(self._Mrho) + pc = MassMatrixPreconditioner(self._Mrho, dim_reduce=None) self._Mrho_inv = inverse( self._Mrho, "pcg", @@ -325,8 +325,9 @@ def _initialize_projectors_and_mass(self): self.curlPibT = self.curlPib.T # Inverse mass matrix needed to compute the error - self.pc_Mv = preconditioner.MassMatrixDiagonalPreconditioner( + self.pc_Mv = preconditioner.MassMatrixPreconditioner( self.mass_ops.Mv, + dim_reduce=None, ) self._inv_Mv = inverse( self.mass_ops.Mv, diff --git a/src/struphy/propagators/variational_momentum_advection.py b/src/struphy/propagators/variational_momentum_advection.py index 47af8f5ad..7fc6e64ce 100644 --- a/src/struphy/propagators/variational_momentum_advection.py +++ b/src/struphy/propagators/variational_momentum_advection.py @@ -8,7 +8,7 @@ from scope_profiler import ProfileManager from struphy.feec import preconditioner -from struphy.feec.preconditioner import MassMatrixDiagonalPreconditioner +from struphy.feec.preconditioner import MassMatrixPreconditioner from struphy.feec.variational_utilities import BracketOperator from struphy.io.options import LiteralOptions, OptionsBase from struphy.linear_algebra.solver import NonlinearSolverParameters, SolverParameters @@ -130,7 +130,7 @@ def allocate(self): # assembly of WMMnew happens in VariationalDensityEvolve self._Mrho = self.mass_ops.WMMnew - pc = MassMatrixDiagonalPreconditioner(self._Mrho) + pc = MassMatrixPreconditioner(self._Mrho, dim_reduce=None) self._Mrho_inv = inverse( self._Mrho, self.options.solver, @@ -276,8 +276,9 @@ def __call_picard(self, dt): def _initialize_mass(self): """Initialization of the mass matrix solver""" # Inverse mass matrix needed to compute the error - self.pc_Mv = preconditioner.MassMatrixDiagonalPreconditioner( + self.pc_Mv = preconditioner.MassMatrixPreconditioner( self.mass_ops.Mv, + dim_reduce=None, ) self._inv_Mv = inverse( self.mass_ops.Mv, diff --git a/src/struphy/propagators/variational_pb_evolve.py b/src/struphy/propagators/variational_pb_evolve.py index a28f879c9..5d5da96e6 100644 --- a/src/struphy/propagators/variational_pb_evolve.py +++ b/src/struphy/propagators/variational_pb_evolve.py @@ -12,7 +12,7 @@ from struphy.feec import preconditioner from struphy.feec.mass import L2Projector -from struphy.feec.preconditioner import MassMatrixDiagonalPreconditioner +from struphy.feec.preconditioner import MassMatrixPreconditioner from struphy.feec.variational_utilities import ( Hdiv0_transport_operator, Pressure_transport_operator, @@ -247,7 +247,7 @@ def allocate(self): # assembly of WMMnew happens in VariationalDensityEvolve self._Mrho = self.mass_ops.WMMnew - pc = MassMatrixDiagonalPreconditioner(self._Mrho) + pc = MassMatrixPreconditioner(self._Mrho, dim_reduce=None) self._Mrho_inv = inverse( self._Mrho, "pcg", @@ -524,8 +524,9 @@ def _initialize_projectors_and_mass(self): self._tmp_int_grid = xp.zeros(grid_shape, dtype=float) # Inverse mass matrix needed to compute the error - self.pc_Mv = preconditioner.MassMatrixDiagonalPreconditioner( + self.pc_Mv = preconditioner.MassMatrixPreconditioner( self.mass_ops.Mv, + dim_reduce=None, ) self._inv_Mv = inverse( self.mass_ops.Mv, diff --git a/src/struphy/propagators/variational_qb_evolve.py b/src/struphy/propagators/variational_qb_evolve.py index 2d6c9f53f..868771f0d 100644 --- a/src/struphy/propagators/variational_qb_evolve.py +++ b/src/struphy/propagators/variational_qb_evolve.py @@ -11,7 +11,7 @@ from struphy.feec import preconditioner from struphy.feec.mass import L2Projector -from struphy.feec.preconditioner import MassMatrixDiagonalPreconditioner +from struphy.feec.preconditioner import MassMatrixPreconditioner from struphy.feec.variational_utilities import ( Hdiv0_transport_operator, Pressure_transport_operator, @@ -248,7 +248,7 @@ def allocate(self): # assembly of WMMnew happens in VariationalDensityEvolve self._Mrho = self.mass_ops.WMMnew - pc = MassMatrixDiagonalPreconditioner(self._Mrho) + pc = MassMatrixPreconditioner(self._Mrho, dim_reduce=None) self._Mrho_inv = inverse( self._Mrho, "pcg", @@ -518,8 +518,9 @@ def _initialize_projectors_and_mass(self): self._tmp_int_grid = xp.zeros(grid_shape, dtype=float) # Inverse mass matrix needed to compute the error - self.pc_Mv = preconditioner.MassMatrixDiagonalPreconditioner( + self.pc_Mv = preconditioner.MassMatrixPreconditioner( self.mass_ops.Mv, + dim_reduce=None, ) self._inv_Mv = inverse( self.mass_ops.Mv, diff --git a/src/struphy/propagators/variational_resistivity.py b/src/struphy/propagators/variational_resistivity.py index 758a654be..b7830e086 100644 --- a/src/struphy/propagators/variational_resistivity.py +++ b/src/struphy/propagators/variational_resistivity.py @@ -124,7 +124,7 @@ class Options(OptionsBase): Adiabatic index. solver : LiteralOptions.OptsSymmSolver, default="pcg" Linear solver for implicit subproblems. - precond : LiteralOptions.OptsMassPrecond, default="MassMatrixDiagonalPreconditioner" + precond : LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner" Preconditioner used in linear solves. solver_params : SolverParameters, default=None Linear-solver controls. @@ -146,7 +146,7 @@ class Options(OptionsBase): model: OptsModel = "full" gamma: float = 5.0 / 3.0 solver: LiteralOptions.OptsSymmSolver = "pcg" - precond: LiteralOptions.OptsMassPrecond = "MassMatrixDiagonalPreconditioner" + precond: LiteralOptions.OptsMassPrecond = "MassMatrixPreconditioner" solver_params: SolverParameters = None nonlin_solver: NonlinearSolverParameters = None linearize_current: bool = False @@ -521,8 +521,9 @@ def __call_newton(self, dt): def _initialize_projectors_and_mass(self): """Initialization of all the `BasisProjectionOperator` and needed to compute the bracket term""" - pc_M1 = preconditioner.MassMatrixDiagonalPreconditioner( + pc_M1 = preconditioner.MassMatrixPreconditioner( self.mass_ops.M1, + dim_reduce=None, ) inv_M1 = inverse( self.mass_ops.M1, @@ -533,8 +534,9 @@ def _initialize_projectors_and_mass(self): verbose=False, ) - pc_M3 = preconditioner.MassMatrixDiagonalPreconditioner( + pc_M3 = preconditioner.MassMatrixPreconditioner( self.mass_ops.M3, + dim_reduce=None, ) self._inv_M3 = inverse( self.mass_ops.M3, @@ -558,7 +560,8 @@ def _initialize_projectors_and_mass(self): preconditioner, self.options.precond, ) - self.pc_jac = pc_class(self.M_de_ds) + # dim_reduce=None: the Kronecker factors do not depend on the weights, cheap to update + self.pc_jac = pc_class(self.M_de_ds, dim_reduce=None) self.inv_jac = inverse( self.M_de_ds, @@ -588,7 +591,7 @@ def _initialize_projectors_and_mass(self): preconditioner, self.options.precond, ) - self.pc = pc_class(M2) + self.pc = pc_class(M2, dim_reduce=None) self.inv_lop = inverse( self.l_op, diff --git a/src/struphy/propagators/variational_viscosity.py b/src/struphy/propagators/variational_viscosity.py index a9b8f11c7..7d8a4ad5c 100644 --- a/src/struphy/propagators/variational_viscosity.py +++ b/src/struphy/propagators/variational_viscosity.py @@ -11,7 +11,7 @@ from struphy.feec import preconditioner from struphy.feec.basis_projection_ops import CoordinateProjector from struphy.feec.mass import L2Projector -from struphy.feec.preconditioner import MassMatrixDiagonalPreconditioner +from struphy.feec.preconditioner import MassMatrixPreconditioner from struphy.feec.variational_utilities import InternalEnergyEvaluator from struphy.io.options import LiteralOptions, OptionsBase from struphy.linear_algebra.solver import NonlinearSolverParameters, SolverParameters @@ -130,7 +130,7 @@ class Options(OptionsBase): Adiabatic index. solver : LiteralOptions.OptsSymmSolver, default="pcg" Linear solver for implicit subproblems. - precond : LiteralOptions.OptsMassPrecond, default="MassMatrixDiagonalPreconditioner" + precond : LiteralOptions.OptsMassPrecond, default="MassMatrixPreconditioner" Preconditioner used in linear solves. solver_params : SolverParameters, default=None Linear-solver controls. @@ -152,7 +152,7 @@ class Options(OptionsBase): model: OptsModel = "full" gamma: float = 5.0 / 3.0 solver: LiteralOptions.OptsSymmSolver = "pcg" - precond: LiteralOptions.OptsMassPrecond = "MassMatrixDiagonalPreconditioner" + precond: LiteralOptions.OptsMassPrecond = "MassMatrixPreconditioner" solver_params: SolverParameters = None nonlin_solver: NonlinearSolverParameters = None mu: float = 0.0 @@ -200,7 +200,7 @@ def allocate(self): # assembly of WMMnew happens in VariationalDensityEvolve self._Mrho = self.mass_ops.WMMnew - pc = MassMatrixDiagonalPreconditioner(self._Mrho) + pc = MassMatrixPreconditioner(self._Mrho, dim_reduce=None) self._Mrho_inv = inverse( self._Mrho, "pcg", @@ -464,8 +464,9 @@ def _initialize_projectors_and_mass(self): M1 = self.mass_ops.M1 self.M1_du = self.mass_ops.create_weighted_mass("Hcurl", "Hcurl") - self.pc_M3 = preconditioner.MassMatrixDiagonalPreconditioner( + self.pc_M3 = preconditioner.MassMatrixPreconditioner( self.mass_ops.M3, + dim_reduce=None, ) self._inv_M3 = inverse( self.mass_ops.M3, @@ -485,7 +486,8 @@ def _initialize_projectors_and_mass(self): preconditioner, self.options.precond, ) - self.pc_jac = pc_class(self.M_de_ds) + # dim_reduce=None: the Kronecker factors do not depend on the weights, cheap to update + self.pc_jac = pc_class(self.M_de_ds, dim_reduce=None) self.inv_jac = inverse( self.M_de_ds,