diff --git a/docs/changelog.txt b/docs/changelog.txt index 52562f54..661cec33 100644 --- a/docs/changelog.txt +++ b/docs/changelog.txt @@ -39,6 +39,8 @@ * ex33, ex33p is added. * ex15 is updated to use Numba + * Initial support of PETSc + MFEM + 2022 03. * Added --cuda-arch option to specify the compute cuda archtecture diff --git a/docs/install.txt b/docs/install.txt index c3a7224f..5bbbf8bb 100644 --- a/docs/install.txt +++ b/docs/install.txt @@ -121,6 +121,10 @@ If you are sure, you could use --skip-swig option, so that it compiles the wrapp codes without re-generating it. $ python setup.py install --with-parallel --skip-ext --skip-swig --mfem-branch="master" +## PETSC (experimental) + +$ pip install petsc +$ python setup.py install --with-parallel --with-petsc ## Other options --unverifiedSSL : diff --git a/examples/petsc/ex1p.py b/examples/petsc/ex1p.py new file mode 100644 index 00000000..4befa291 --- /dev/null +++ b/examples/petsc/ex1p.py @@ -0,0 +1,177 @@ +''' + MFEM example 1p + + See c++ version in the MFEM library for more detail +''' +import sys +import os +from os.path import expanduser, join +import numpy as np + +from mpi4py import MPI + +from mfem.common.arg_parser import ArgParser +from mfem import path +import mfem.par as mfem + +def_meshfile = expanduser(join(path, 'data', 'star.mesh')) + +num_proc = MPI.COMM_WORLD.size +myid = MPI.COMM_WORLD.rank + + +def run(order=1, static_cond=False, + meshfile=def_meshfile, + visualization=False, + device='cpu', + pa=False): + + algebraic_ceed = False + + device = mfem.Device(device) + if myid == 0: + device.Print() + + mfem.MFEMInitializePetsc(None, None, "", None) + + mesh = mfem.Mesh(meshfile, 1, 1) + dim = mesh.Dimension() + + ref_levels = int(np.floor(np.log(10000./mesh.GetNE())/np.log(2.)/dim)) + for x in range(ref_levels): + mesh.UniformRefinement() + mesh.ReorientTetMesh() + pmesh = mfem.ParMesh(MPI.COMM_WORLD, mesh) + del mesh + + par_ref_levels = 2 + for l in range(par_ref_levels): + pmesh.UniformRefinement() + + if order > 0: + fec = mfem.H1_FECollection(order, dim) + elif mesh.GetNodes(): + fec = mesh.GetNodes().OwnFEC() + print("Using isoparametric FEs: " + str(fec.Name())) + else: + order = 1 + fec = mfem.H1_FECollection(order, dim) + + fespace = mfem.ParFiniteElementSpace(pmesh, fec) + fe_size = fespace.GlobalTrueVSize() + + if (myid == 0): + print('Number of finite element unknowns: ' + str(fe_size)) + + ess_tdof_list = mfem.intArray() + if pmesh.bdr_attributes.Size() > 0: + ess_bdr = mfem.intArray(pmesh.bdr_attributes.Max()) + ess_bdr.Assign(1) + fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list) + + # the basis functions in the finite element fespace. + b = mfem.ParLinearForm(fespace) + one = mfem.ConstantCoefficient(1.0) + b.AddDomainIntegrator(mfem.DomainLFIntegrator(one)) + b.Assemble() + + x = mfem.ParGridFunction(fespace) + x.Assign(0.0) + + a = mfem.ParBilinearForm(fespace) + if pa: + a.SetAssemblyLevel(mfem.AssemblyLevel_PARTIAL) + + a.AddDomainIntegrator(mfem.DiffusionIntegrator(one)) + + if static_cond: + a.EnableStaticCondensation() + + a.Assemble() + + A = mfem.OperatorPtr() + B = mfem.Vector() + X = mfem.Vector() + a.FormLinearSystem(ess_tdof_list, x, b, A, X, B) + + if pa: + if mfem.UsesTensorBasis(fespace): + if algebraic_ceed: + assert False, "not supported" + # prec = new ceed::AlgebraicSolver(a, ess_tdof_list); + else: + prec = mfem.OperatorJacobiSmoother(a, ess_tdof_list) + else: + prec = mfem.HypreBoomerAMG() + + ''' + cg = mfem.CGSolver(MPI.COMM_WORLD) + cg.SetRelTol(1e-12) + cg.SetMaxIter(200) + cg.SetPrintLevel(1) + cg.SetPreconditioner(prec) + cg.SetOperator(A.Ptr()) + cg.Mult(B, X) + ''' + hA = A.AsHypreParMatrix() + pcg = mfem.PetscPCGSolver(hA, False) + prec = mfem.HypreBoomerAMG(hA) + pcg.SetPreconditioner(prec) + pcg.iterative_mode = True + pcg.SetRelTol(1e-12) + pcg.SetAbsTol(1e-12) + pcg.SetMaxIter(200) + pcg.SetPrintLevel(1) + pcg.Mult(B, X) + + a.RecoverFEMSolution(X, b, x) + + smyid = '{:0>6d}'.format(myid) + mesh_name = "mesh."+smyid + sol_name = "sol."+smyid + + pmesh.Print(mesh_name, 8) + x.Save(sol_name, 8) + + +if __name__ == "__main__": + parser = ArgParser(description='Ex1 (Laplace Problem)') + parser.add_argument('-m', '--mesh', + default='star.mesh', + action='store', type=str, + help='Mesh file to use.') + parser.add_argument('-vis', '--visualization', + action='store_true', + help='Enable GLVis visualization') + parser.add_argument('-o', '--order', + action='store', default=1, type=int, + help="Finite element order (polynomial degree) or -1 for isoparametric space.") + parser.add_argument('-sc', '--static-condensation', + action='store_true', + help="Enable static condensation.") + parser.add_argument("-pa", "--partial-assembly", + action='store_true', + help="Enable Partial Assembly.") + parser.add_argument("-d", "--device", + default="cpu", type=str, + help="Device configuration string, see Device::Configure().") + + args = parser.parse_args() + + if myid == 0: + parser.print_options(args) + + order = args.order + static_cond = args.static_condensation + meshfile = expanduser( + join(os.path.dirname(__file__), '..', '..', 'data', args.mesh)) + visualization = args.visualization + device = args.device + pa = args.partial_assembly + + run(order=order, + static_cond=static_cond, + meshfile=meshfile, + visualization=visualization, + device=device, + pa=pa) diff --git a/examples/petsc/rc_ex1p b/examples/petsc/rc_ex1p new file mode 100644 index 00000000..3aa2c36a --- /dev/null +++ b/examples/petsc/rc_ex1p @@ -0,0 +1,17 @@ +# Options for the Krylov solver +-ksp_view +-ksp_converged_reason +# Options for the preconditioner +-pc_type gamg +-pc_gamg_type agg +-pc_gamg_agg_nsmooths 1 +-pc_gamg_coarse_eq_limit 100 +-pc_gamg_reuse_interpolation +-pc_gamg_square_graph 1 +-pc_gamg_threshold 0.0 +-mg_levels_ksp_max_it 2 +-mg_levels_ksp_type chebyshev +-mg_levels_esteig_ksp_type cg +-mg_levels_esteig_ksp_max_it 10 +-mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 +-mg_levels_pc_type sor diff --git a/mfem/_par/handle.i b/mfem/_par/handle.i index f83829b8..d1e8cd2c 100644 --- a/mfem/_par/handle.i +++ b/mfem/_par/handle.i @@ -12,6 +12,7 @@ #include "fem/pfespace.hpp" #include "numpy/arrayobject.h" #include "pyoperator.hpp" +#include "../common/pycoefficient.hpp" %} %include "../common/mfem_config.i" @@ -34,7 +35,7 @@ import_array(); %import "hypre.i" #endif #ifdef MFEM_USE_PETSC -%include "petsc.i" +%import "petsc.i" #endif %import "mem_manager.i" diff --git a/mfem/_par/petsc.i b/mfem/_par/petsc.i new file mode 100644 index 00000000..8cbc3dd0 --- /dev/null +++ b/mfem/_par/petsc.i @@ -0,0 +1,43 @@ +%module(package="mfem._par") petsc +%{ +#include "mfem.hpp" + +#include "../common/pycoefficient.hpp" +#include "numpy/arrayobject.h" +#include "pyoperator.hpp" +#include "../common/io_stream.hpp" +%} + +%include "../common/mfem_config.i" + +#ifdef MFEM_USE_MPI +%include mpi4py/mpi4py.i +%mpi4py_typemap(Comm, MPI_Comm); +#endif + +%init %{ +import_array(); +%} + +%include "exception.i" +%import "fe.i" +%import "fe_fixed_order.i" +%import "element.i" +%import "mesh.i" +%import "operators.i" +%import "coefficient.i" + +%include "../common/typemap_macros.i" +%include "../common/exception.i" + +%import "../common/io_stream_typemap.i" +OSTREAM_TYPEMAP(std::ostream&) +ISTREAM_TYPEMAP(std::istream&) + +%import "petscconf.h" +%include "linalg/petsc.hpp" + + + + + diff --git a/mfem/_par/pyoperator.hpp b/mfem/_par/pyoperator.hpp index d6e78549..1f65fb79 100644 --- a/mfem/_par/pyoperator.hpp +++ b/mfem/_par/pyoperator.hpp @@ -1,3 +1,6 @@ +#ifndef PYMFEM_PYOPERATOR +#define PYMFEM_PYOPERATOR + #include "linalg/operator.hpp" namespace mfem{ @@ -20,3 +23,5 @@ class PyOperatorBase : public Operator virtual ~PyTimeDependentOperatorBase() {} }; } + +#endif diff --git a/mfem/_par/setup.py b/mfem/_par/setup.py index aa430497..b84699e3 100644 --- a/mfem/_par/setup.py +++ b/mfem/_par/setup.py @@ -44,10 +44,11 @@ def get_extensions(): from setup_local import (mfembuilddir, mfemincdir, mfemsrcdir, mfemlnkdir, mfemptpl, build_mfem, hypreinc, metisinc, hyprelib, metis5lib, + petscinc, petsclib, cc_par, cxx_par, cxx11flag, add_pumi, add_cuda, add_libceed, add_strumpack, - add_suitesparse, add_gslibp) + add_suitesparse, add_gslibp, add_petsc) include_dirs = [mfembuilddir, mfemincdir, mfemsrcdir, numpy.get_include(), @@ -154,6 +155,11 @@ def get_extensions(): include_dirs.append(gslibpinc) modules.append("gslib") + if add_petsc == '1': + include_dirs.append(petscinc) + library_dirs.append(petsclib) + modules.append("petsc") + sources = {name: [name + "_wrap.cxx"] for name in modules} proxy_names = {name: '_'+name for name in modules} @@ -173,6 +179,7 @@ def get_extensions(): else: runtime_library_dirs = library_dirs + ext_modules = [Extension(proxy_names[modules[0]], sources=sources[modules[0]], extra_compile_args=extra_compile_args, diff --git a/mfem/common/io_stream.hpp b/mfem/common/io_stream.hpp index a479cae1..98347b5f 100644 --- a/mfem/common/io_stream.hpp +++ b/mfem/common/io_stream.hpp @@ -1,3 +1,6 @@ +#ifndef PYMFEM_IO_STREAM +#define PYMFEM_IO_STREAM + namespace PyMFEM { class wFILE @@ -45,3 +48,4 @@ class wFILE }; } +#endif diff --git a/mfem/common/pycoefficient.hpp b/mfem/common/pycoefficient.hpp index c4088db3..4ca86f4a 100644 --- a/mfem/common/pycoefficient.hpp +++ b/mfem/common/pycoefficient.hpp @@ -1,3 +1,6 @@ +#ifndef PYMFEM_PYCOEFFICIENT +#define PYMFEM_PYCOEFFICIENT + #include "fem/coefficient.hpp" #define assertm(exp, msg) assert(((void)msg, exp)) @@ -191,3 +194,4 @@ class MatrixNumbaCoefficient : public mfem::MatrixFunctionCoefficient, public N }; +#endif diff --git a/mfem/par.py b/mfem/par.py index 514aa344..46061fa0 100644 --- a/mfem/par.py +++ b/mfem/par.py @@ -97,6 +97,10 @@ from mfem._par.gslib import * except: pass +try: + from mfem._par.petsc import * +except: + pass import mfem._par.array as array import mfem._par.blockoperator as blockoperator diff --git a/setup.py b/setup.py index d568a708..99a81349 100644 --- a/setup.py +++ b/setup.py @@ -120,6 +120,9 @@ gslibp_prefix = '' gslib_only = False +enable_petsc = False +petsc_prefix = "" + enable_suitesparse = False suitesparse_prefix = "" @@ -810,8 +813,20 @@ def add_rpath(p): if suitesparse_prefix != '': cmake_opts['DSuiteSparse_DIR'] = suitesparse_prefix + if not serial and enable_petsc: + cmake_opts['DMFEM_USE_PETSC'] = '1' + cmake_opts['DPETSC_DIR'] = petsc_prefix + petsclibpath = os.path.dirname( + find_libpath_from_prefix( + 'petsc', petsc_prefix)) + cmake_opts['DPETSC_LIBRARIES'] = "-L" + petsclibpath + " -lpetsc" + cmake_opts['DPETSC_INCLUDES'] = os.path.join(petsc_prefix, 'include') + cmake_opts['DPETSC_ARCH'] = "" + cmake_opts["DPETSC_EXECUTABLE_RUNS"] = "YES" + if enable_lapack: cmake_opts['DMFEM_USE_LAPACK'] = '1' + if blas_libraries != "": cmake_opts['DBLAS_LIBRARIES'] = blas_libraries if lapack_libraries != "": @@ -901,6 +916,7 @@ def write_setup_local(): 'add_gslib': '', 'add_gslibp': '', 'add_gslibs': '', + 'add_petsc': '', 'libceedinc': os.path.join(libceed_prefix, 'include'), 'gslibsinc': os.path.join(gslibs_prefix, 'include'), 'gslibpinc': os.path.join(gslibp_prefix, 'include'), @@ -941,6 +957,8 @@ def add_extra(xxx, inc_sub=None): add_extra('gslibs') if enable_gslib: add_extra('gslibp') + if enable_petsc: + add_extra('petsc') pwd = chdir(rootdir) @@ -1062,6 +1080,8 @@ def update_header_exists(mfem_source): parflag.append('-I' + os.path.join(pumi_prefix, 'include')) if enable_strumpack: parflag.append('-I' + os.path.join(strumpack_prefix, 'include')) + if enable_petsc: + parflag.append('-I' + os.path.join(petsc_prefix, 'include')) if enable_suitesparse: parflag.append('-I' + os.path.join(suitesparse_prefix, 'include', 'suitesparse')) @@ -1194,6 +1214,8 @@ def print_config(): print(" hypre prefix", hypre_prefix) print(" metis prefix", metis_prefix) + if enable_petsc: + print(" petsc prefix", petsc_prefix) print(" c compiler : " + cc_command) print(" c++ compiler : " + cxx_command) print(" mpi-c compiler : " + mpicc_command) @@ -1230,6 +1252,7 @@ def configure_install(self): global enable_libceed, libceed_prefix, libceed_only global enable_gslib, gslibs_prefix, gslibp_prefix, gslib_only global enable_suitesparse, suitesparse_prefix + global enable_petsc, petsc_prefix global enable_lapack, blas_libraries, lapack_libraries verbose = bool(self.vv) if verbose == -1 else verbose @@ -1259,8 +1282,10 @@ def configure_install(self): enable_gslib = bool(self.with_gslib) gslib_only = bool(self.gslib_only) enable_suitesparse = bool(self.with_suitesparse) + enable_petsc = bool(self.with_petsc) enable_lapack = bool(self.with_lapack) + build_parallel = bool(self.with_parallel) # controlls PyMFEM parallel build_serial = not bool(self.no_serial) @@ -1358,6 +1383,11 @@ def configure_install(self): else: strumpack_prefix = mfem_prefix + if self.petsc_prefix != '': + petsc_prefix = abspath(self.petsc_prefix) + else: + assert not enable_petsc, "Petsc is enabled but prefix is not given" + if enable_cuda: nvcc = find_command('nvcc') cuda_prefix = os.path.dirname(os.path.dirname(nvcc)) @@ -1513,23 +1543,26 @@ class Install(_install): ('cuda-arch=', None, 'set cuda compute capability. Ex if A100, set to 80'), ('with-metis64', None, 'use 64bit int in metis'), ('with-pumi', None, 'enable pumi (parallel only)'), - ('pumi-prefix=', None, 'Specify locaiton of pumi'), + ('pumi-prefix=', None, 'locaiton of pumi'), ('with-suitesparse', None, 'build MFEM with suitesparse (MFEM_USE_SUITESPARSE=YES) (parallel only)'), ('suitesparse-prefix=', None, 'Specify locaiton of suitesparse (=SuiteSparse_DIR)'), ('with-libceed', None, 'enable libceed'), - ('libceed-prefix=', None, 'Specify locaiton of libceed'), + ('libceed-prefix=', None, 'locaiton of libceed'), ('libceed-only', None, 'Build libceed only'), - ('gslib-prefix=', None, 'Specify locaiton of gslib'), + ('gslib-prefix=', None, 'locaiton of gslib'), ('with-gslib', None, 'enable gslib'), ('gslib-only', None, 'Build gslib only'), ('with-strumpack', None, 'enable strumpack (parallel only)'), + ('strumpack-prefix=', None, 'Specify locaiton of strumpack'), ('with-lapack', None, 'build MFEM with lapack'), ('blas-libraries=', None, 'Specify locaiton of Blas library (used to build MFEM)'), ('lapack-libraries=', None, - 'Specify locaiton of Lapack library (used to build MFEM)'), + 'locaiton of Lapack library (used to build MFEM)'), + ('with-petsc', None, 'enable petsc (default = on if import petsc works)'), + ('petsc-prefilx=', None, 'location of petsc'), ] def initialize_options(self): @@ -1572,6 +1605,14 @@ def initialize_options(self): self.libceed_prefix = '' self.libceed_only = False + try: + import petsc + self.petsc_prefix = petsc.get_petsc_dir() + self.with_petsc = True + except ImportError: + self.with_petsc = False + self.petsc_prefix = "" + self.with_gslib = False self.gslib_prefix = '' self.gslib_only = False