From bfe20734d997658ea9f2d4525e96e9824dd352bb Mon Sep 17 00:00:00 2001 From: Syun'ichi Shiraiwa Date: Wed, 19 Nov 2025 23:37:52 -0800 Subject: [PATCH 1/6] (WIP) serial version --- _build_system/build_utils.py | 17 +- mfem/_ser/bilininteg_hdg.i | 25 ++ mfem/_ser/darcyform.i | 30 +++ mfem/_ser/darcyhybridization.i | 25 ++ mfem/_ser/darcyreduction.i | 25 ++ mfem/_ser/estimators_hdg.i | 25 ++ mfem/_ser/nonlininteg.i | 5 +- mfem/_ser/setup.py | 4 +- mfem/common/bilininteg_ext.i | 317 ------------------------- mfem/common/generate_bilininteg_ext.py | 2 + mfem/common/lininteg_ext.i | 78 ------ 11 files changed, 150 insertions(+), 403 deletions(-) create mode 100644 mfem/_ser/bilininteg_hdg.i create mode 100644 mfem/_ser/darcyform.i create mode 100644 mfem/_ser/darcyhybridization.i create mode 100644 mfem/_ser/darcyreduction.i create mode 100644 mfem/_ser/estimators_hdg.i delete mode 100644 mfem/common/bilininteg_ext.i delete mode 100644 mfem/common/lininteg_ext.i diff --git a/_build_system/build_utils.py b/_build_system/build_utils.py index 0cff67ef..6f4ebd30 100644 --- a/_build_system/build_utils.py +++ b/_build_system/build_utils.py @@ -102,19 +102,26 @@ def make_call(command, target='', force_verbose=False, env=None): env.update(os.environ) myverbose = bglb.verbose or force_verbose - if not myverbose: - kwargs['stdout'] = subprocess.DEVNULL - kwargs['stderr'] = subprocess.DEVNULL - + + kwargs['stdout'] = subprocess.PIPE + kwargs['stderr'] = subprocess.PIPE + p = subprocess.Popen(command, **kwargs) - p.communicate() + stdout, stderr= p.communicate() if p.returncode != 0: if target == '': target = " ".join(command) + + print("!!!!!!!!!!!!!!!", p) + print(stdout) + print(stderr) + print("Failed when calling command: " + target) raise subprocess.CalledProcessError(p.returncode, " ".join(command)) + if not myverbose: + print(stdout) def chdir(path): ''' diff --git a/mfem/_ser/bilininteg_hdg.i b/mfem/_ser/bilininteg_hdg.i new file mode 100644 index 00000000..fed4ce1c --- /dev/null +++ b/mfem/_ser/bilininteg_hdg.i @@ -0,0 +1,25 @@ +// +// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved. +// +%module(package="mfem._ser") bilininteg_hdg + +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "../common/pyoperator.hpp" +#include "../common/pyintrules.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_FEM_DARCY_BILININTEG_HDG + +%init %{ +import_array(); +%} +%include "exception.i" +%include "../common/exception.i" + +%include "fem/darcy/bilininteg_hdg.hpp" + +#endif + diff --git a/mfem/_ser/darcyform.i b/mfem/_ser/darcyform.i new file mode 100644 index 00000000..c1c4bf7d --- /dev/null +++ b/mfem/_ser/darcyform.i @@ -0,0 +1,30 @@ +// +// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved. +// +%module(package="mfem._ser") darcyform + +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "../common/pyoperator.hpp" +#include "../common/pyintrules.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_FEM_DARCY_DARCYFORM + +%init %{ +import_array(); +%} +%include "exception.i" +%include "../common/exception.i" + +%import "fespace.i" +%import "blockvector.i" +%import "bilinearform.i" +%import "nonlinearform.i" + +%include "fem/darcy/darcyform.hpp" + +#endif + diff --git a/mfem/_ser/darcyhybridization.i b/mfem/_ser/darcyhybridization.i new file mode 100644 index 00000000..d95ad474 --- /dev/null +++ b/mfem/_ser/darcyhybridization.i @@ -0,0 +1,25 @@ +// +// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved. +// +%module(package="mfem._ser") darcyhybridization + +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "../common/pyoperator.hpp" +#include "../common/pyintrules.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_FEM_DARCY_DARCYHYBRIDIZATION + +%init %{ +import_array(); +%} +%include "exception.i" +%include "../common/exception.i" + +%include "fem/darcy/darcyhybridization.hpp" + +#endif + diff --git a/mfem/_ser/darcyreduction.i b/mfem/_ser/darcyreduction.i new file mode 100644 index 00000000..a3a1c4c2 --- /dev/null +++ b/mfem/_ser/darcyreduction.i @@ -0,0 +1,25 @@ +// +// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved. +// +%module(package="mfem._ser") darcyreduction + +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "../common/pyoperator.hpp" +#include "../common/pyintrules.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_FEM_DARCY_DARCYREDUCTION + +%init %{ +import_array(); +%} +%include "exception.i" +%include "../common/exception.i" + +%include "fem/darcy/darcyreduction.hpp" + +#endif + diff --git a/mfem/_ser/estimators_hdg.i b/mfem/_ser/estimators_hdg.i new file mode 100644 index 00000000..476c3826 --- /dev/null +++ b/mfem/_ser/estimators_hdg.i @@ -0,0 +1,25 @@ +// +// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved. +// +%module(package="mfem._ser") estimators_hdg + +%{ +#include "mfem.hpp" +#include "numpy/arrayobject.h" +#include "../common/pyoperator.hpp" +#include "../common/pyintrules.hpp" +%} + +%include "../common/existing_mfem_headers.i" +#ifdef FILE_EXISTS_FEM_DARCY_ESTIMATORS_HDG + +%init %{ +import_array(); +%} +%include "exception.i" +%include "../common/exception.i" + +%include "fem/darcy/estimators_hdg.hpp" + +#endif + diff --git a/mfem/_ser/nonlininteg.i b/mfem/_ser/nonlininteg.i index 0b746101..118322e3 100644 --- a/mfem/_ser/nonlininteg.i +++ b/mfem/_ser/nonlininteg.i @@ -3,6 +3,7 @@ // %module(package="mfem._ser", directors="1") nonlininteg %{ +#include "numpy/arrayobject.h" #include "mfem.hpp" #include "../common/pycoefficient.hpp" #include "../common/pyoperator.hpp" @@ -10,11 +11,11 @@ #include "../common/pynonlininteg.hpp" using namespace mfem; %} -/* + %init %{ import_array(); %} -*/ + /* %exception { try { $action } diff --git a/mfem/_ser/setup.py b/mfem/_ser/setup.py index 86b0c81a..0a062e21 100644 --- a/mfem/_ser/setup.py +++ b/mfem/_ser/setup.py @@ -110,7 +110,9 @@ def get_extensions(): "attribute_sets", "arrays_by_name", "hyperbolic", "complex_densemat", "complexstaticcond", "complexweakform", - "bounds", "integrator"] + "bounds", "integrator", + "darcyform","estimators_hdg", "bilininteg_hdg", + "darcyreductoin", "darcyhybridization"] if add_cuda == '1': from setup_local import cudainc diff --git a/mfem/common/bilininteg_ext.i b/mfem/common/bilininteg_ext.i deleted file mode 100644 index 0929cdc1..00000000 --- a/mfem/common/bilininteg_ext.i +++ /dev/null @@ -1,317 +0,0 @@ -namespace mfem { -%pythonappend BilinearFormIntegrator::BilinearFormIntegrator %{ - self._coeff = args -%} - -%pythonappend TransposeIntegrator::TransposeIntegrator %{ - if own_bfi_ == 1: bfi_.thisown = 0 -%} - -%pythonappend LumpedIntegrator::LumpedIntegrator %{ - if own_bfi_ == 1: bfi_.thisown = 0 -%} - -%pythonappend InverseIntegrator::InverseIntegrator %{ - if own_integ == 1: integ.thisown = 0 -%} - -%pythonappend SumIntegrator::SumIntegrator %{ - self.own_integs = own_integs -%} - -%pythonappend MixedScalarIntegrator::MixedScalarIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedVectorIntegrator::MixedVectorIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedScalarVectorIntegrator::MixedScalarVectorIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedScalarMassIntegrator::MixedScalarMassIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedVectorProductIntegrator::MixedVectorProductIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedScalarDerivativeIntegrator::MixedScalarDerivativeIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedScalarWeakDerivativeIntegrator::MixedScalarWeakDerivativeIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedScalarDivergenceIntegrator::MixedScalarDivergenceIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedVectorDivergenceIntegrator::MixedVectorDivergenceIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedScalarWeakGradientIntegrator::MixedScalarWeakGradientIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedScalarCurlIntegrator::MixedScalarCurlIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedScalarWeakCurlIntegrator::MixedScalarWeakCurlIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedVectorMassIntegrator::MixedVectorMassIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedCrossProductIntegrator::MixedCrossProductIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedDotProductIntegrator::MixedDotProductIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedWeakGradDotIntegrator::MixedWeakGradDotIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedWeakDivCrossIntegrator::MixedWeakDivCrossIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedGradGradIntegrator::MixedGradGradIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedCrossGradGradIntegrator::MixedCrossGradGradIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedCurlCurlIntegrator::MixedCurlCurlIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedCrossCurlCurlIntegrator::MixedCrossCurlCurlIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedCrossCurlGradIntegrator::MixedCrossCurlGradIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedCrossGradCurlIntegrator::MixedCrossGradCurlIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedWeakCurlCrossIntegrator::MixedWeakCurlCrossIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedScalarWeakCurlCrossIntegrator::MixedScalarWeakCurlCrossIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedCrossGradIntegrator::MixedCrossGradIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedCrossCurlIntegrator::MixedCrossCurlIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedScalarCrossCurlIntegrator::MixedScalarCrossCurlIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedScalarCrossGradIntegrator::MixedScalarCrossGradIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedScalarCrossProductIntegrator::MixedScalarCrossProductIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedScalarWeakCrossProductIntegrator::MixedScalarWeakCrossProductIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedDirectionalDerivativeIntegrator::MixedDirectionalDerivativeIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedGradDivIntegrator::MixedGradDivIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedDivGradIntegrator::MixedDivGradIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedScalarWeakDivergenceIntegrator::MixedScalarWeakDivergenceIntegrator %{ - self._coeff = vq -%} - -%pythonappend MixedVectorGradientIntegrator::MixedVectorGradientIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedVectorCurlIntegrator::MixedVectorCurlIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedVectorWeakCurlIntegrator::MixedVectorWeakCurlIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedVectorWeakDivergenceIntegrator::MixedVectorWeakDivergenceIntegrator %{ - self._coeff = args -%} - -%pythonappend GradientIntegrator::GradientIntegrator %{ - self._coeff = args -%} - -%pythonappend DiffusionIntegrator::DiffusionIntegrator %{ - self._coeff = args -%} - -%pythonappend MassIntegrator::MassIntegrator %{ - self._coeff = args -%} - -%pythonappend BoundaryMassIntegrator::BoundaryMassIntegrator %{ - self._coeff = q -%} - -%pythonappend ConvectionIntegrator::ConvectionIntegrator %{ - self._coeff = q -%} - -%pythonappend ConservativeConvectionIntegrator::ConservativeConvectionIntegrator %{ - self._coeff = q -%} - -%pythonappend GroupConvectionIntegrator::GroupConvectionIntegrator %{ - self._coeff = q -%} - -%pythonappend VectorMassIntegrator::VectorMassIntegrator %{ - self._coeff = args -%} - -%pythonappend VectorFEDivergenceIntegrator::VectorFEDivergenceIntegrator %{ - self._coeff = args -%} - -%pythonappend VectorFEWeakDivergenceIntegrator::VectorFEWeakDivergenceIntegrator %{ - self._coeff = args -%} - -%pythonappend VectorFECurlIntegrator::VectorFECurlIntegrator %{ - self._coeff = args -%} - -%pythonappend VectorFEBoundaryFluxIntegrator::VectorFEBoundaryFluxIntegrator %{ - self._coeff = args -%} - -%pythonappend DerivativeIntegrator::DerivativeIntegrator %{ - self._coeff = q -%} - -%pythonappend CurlCurlIntegrator::CurlCurlIntegrator %{ - self._coeff = args -%} - -%pythonappend VectorCurlCurlIntegrator::VectorCurlCurlIntegrator %{ - self._coeff = args -%} - -%pythonappend MixedCurlIntegrator::MixedCurlIntegrator %{ - self._coeff = args -%} - -%pythonappend VectorFEMassIntegrator::VectorFEMassIntegrator %{ - self._coeff = args -%} - -%pythonappend VectorDivergenceIntegrator::VectorDivergenceIntegrator %{ - self._coeff = args -%} - -%pythonappend DivDivIntegrator::DivDivIntegrator %{ - self._coeff = args -%} - -%pythonappend VectorDiffusionIntegrator::VectorDiffusionIntegrator %{ - self._coeff = args -%} - -%pythonappend ElasticityIntegrator::ElasticityIntegrator %{ - self._coeff = args -%} - -%pythonappend ElasticityComponentIntegrator::ElasticityComponentIntegrator %{ - self._coeff = parent_ -%} - -%pythonappend DGTraceIntegrator::DGTraceIntegrator %{ - self._coeff = args -%} - -%pythonappend NonconservativeDGTraceIntegrator::NonconservativeDGTraceIntegrator %{ - self._coeff = args -%} - -%pythonappend DGDiffusionIntegrator::DGDiffusionIntegrator %{ - self._coeff = args -%} - -%pythonappend DGDiffusionBR2Integrator::DGDiffusionBR2Integrator %{ - self._coeff = args -%} - -%pythonappend DGElasticityIntegrator::DGElasticityIntegrator %{ - self._coeff = args -%} - -%pythonappend ScalarProductInterpolator::ScalarProductInterpolator %{ - self._coeff = sc -%} - -%pythonappend ScalarVectorProductInterpolator::ScalarVectorProductInterpolator %{ - self._coeff = sc -%} - -%pythonappend VectorScalarProductInterpolator::VectorScalarProductInterpolator %{ - self._coeff = vc -%} - -%pythonappend ScalarCrossProductInterpolator::ScalarCrossProductInterpolator %{ - self._coeff = vc -%} - -%pythonappend VectorCrossProductInterpolator::VectorCrossProductInterpolator %{ - self._coeff = vc -%} - -%pythonappend VectorInnerProductInterpolator::VectorInnerProductInterpolator %{ - self._coeff = vc -%} - -%pythonappend PyBilinearFormIntegrator::PyBilinearFormIntegrator %{ - self._ir=ir -%} - -%pythonappend SumIntegrator::AddIntegrator %{ - if self.own_integs == 1: integ.thisown = 0 -%} -} diff --git a/mfem/common/generate_bilininteg_ext.py b/mfem/common/generate_bilininteg_ext.py index 0d02825d..7267f841 100644 --- a/mfem/common/generate_bilininteg_ext.py +++ b/mfem/common/generate_bilininteg_ext.py @@ -41,6 +41,8 @@ pass elif line.find("(self, vdim_)") != -1: pass + elif line.find("(self, sign_=1.)") != -1: + pass else: print(cname) print(line) diff --git a/mfem/common/lininteg_ext.i b/mfem/common/lininteg_ext.i deleted file mode 100644 index 7e02196b..00000000 --- a/mfem/common/lininteg_ext.i +++ /dev/null @@ -1,78 +0,0 @@ -namespace mfem { -%pythonappend LinearFormIntegrator::LinearFormIntegrator %{ - self._coeff = args -%} -%pythonappend DeltaLFIntegrator::DeltaLFIntegrator %{ - self._coeff = args -%} -%pythonappend DomainLFIntegrator::DomainLFIntegrator %{ - self._coeff = args -%} -%pythonappend DomainLFGradIntegrator::DomainLFGradIntegrator %{ - self._coeff = QF -%} -%pythonappend BoundaryLFIntegrator::BoundaryLFIntegrator %{ - self._coeff = QG -%} -%pythonappend BoundaryNormalLFIntegrator::BoundaryNormalLFIntegrator %{ - self._coeff = QG -%} -%pythonappend BoundaryTangentialLFIntegrator::BoundaryTangentialLFIntegrator %{ - self._coeff = QG -%} -%pythonappend VectorDomainLFIntegrator::VectorDomainLFIntegrator %{ - self._ir=ir - self._coeff = QF -%} -%pythonappend VectorDomainLFGradIntegrator::VectorDomainLFGradIntegrator %{ - self._coeff = QF -%} -%pythonappend VectorBoundaryLFIntegrator::VectorBoundaryLFIntegrator %{ - self._coeff = QG -%} -%pythonappend VectorFEDomainLFIntegrator::VectorFEDomainLFIntegrator %{ - self._coeff = F -%} -%pythonappend VectorFEDomainLFCurlIntegrator::VectorFEDomainLFCurlIntegrator %{ - self._coeff = F -%} -%pythonappend VectorFEDomainLFDivIntegrator::VectorFEDomainLFDivIntegrator %{ - self._coeff = QF -%} -%pythonappend VectorBoundaryFluxLFIntegrator::VectorBoundaryFluxLFIntegrator %{ - self._ir=ir - self._coeff = f -%} -%pythonappend VectorFEBoundaryFluxLFIntegrator::VectorFEBoundaryFluxLFIntegrator %{ - self._coeff = args -%} -%pythonappend VectorFEBoundaryNormalLFIntegrator::VectorFEBoundaryNormalLFIntegrator %{ - self._coeff = f -%} -%pythonappend VectorFEBoundaryTangentLFIntegrator::VectorFEBoundaryTangentLFIntegrator %{ - self._coeff = QG -%} -%pythonappend BoundaryFlowIntegrator::BoundaryFlowIntegrator %{ - self._coeff = args -%} -%pythonappend DGDirichletLFIntegrator::DGDirichletLFIntegrator %{ - self._coeff = args -%} -%pythonappend DGElasticityDirichletLFIntegrator::DGElasticityDirichletLFIntegrator %{ - self._coeff = uD_ -%} -%pythonappend WhiteGaussianNoiseDomainLFIntegrator::WhiteGaussianNoiseDomainLFIntegrator %{ - self._coeff = QG -%} -%pythonappend VectorQuadratureLFIntegrator::VectorQuadratureLFIntegrator %{ - self._ir=ir - self._coeff = vqfc -%} -%pythonappend QuadratureLFIntegrator::QuadratureLFIntegrator %{ - self._ir=ir - self._coeff = qfc -%} -%pythonappend PyLinearFormIntegrator::PyLinearFormIntegrator %{ - self._ir=ir -%} -} \ No newline at end of file From 903df0a0f4e65e06a8a8623757aef109ccd5ab61 Mon Sep 17 00:00:00 2001 From: Syun'ichi Shiraiwa Date: Fri, 21 Nov 2025 00:31:09 -0500 Subject: [PATCH 2/6] (WIP) serial version built --- _build_system/build_pymfem.py | 12 ++++++++++++ _build_system/build_utils.py | 22 +++++++++------------- mfem/_ser/bilininteg_hdg.i | 6 +++++- mfem/_ser/darcyform.i | 4 ++++ mfem/_ser/darcyhybridization.i | 6 +++++- mfem/_ser/darcyreduction.i | 10 +++++++++- mfem/_ser/estimators_hdg.i | 5 ++++- mfem/_ser/setup.py | 2 +- 8 files changed, 49 insertions(+), 18 deletions(-) diff --git a/_build_system/build_pymfem.py b/_build_system/build_pymfem.py index 140a37be..809303f5 100644 --- a/_build_system/build_pymfem.py +++ b/_build_system/build_pymfem.py @@ -192,6 +192,14 @@ def update_header_exists(mfem_source): serflag.append('-I' + os.path.join(bglb.suitesparse_prefix, 'include', 'suitesparse')) + # generate xxinteg_ext.i + # create empty files + # run swig with *.i files + # call update scripts + for filename in ('../common/bilininteg_ext.i', + '../common/lininteg_ext.i'): + fid = open(filename) + fid.close() for filename in ['lininteg.i', 'bilininteg.i']: command = [swig_command] + swigflag + serflag + [filename] make_call(command) @@ -208,6 +216,10 @@ def update_header_exists(mfem_source): with mp_pool: mp_pool.map(subprocess.run, commands) + #for c in commands: + # print(" ".join(c)) + # subprocess.run(c) + if not do_parallel: os.chdir(pwd) return diff --git a/_build_system/build_utils.py b/_build_system/build_utils.py index 6f4ebd30..2bb98b37 100644 --- a/_build_system/build_utils.py +++ b/_build_system/build_utils.py @@ -19,7 +19,7 @@ __all__ = ["read_mfem_tplflags", "abspath", "external_install_prefix", "make_call", "chdir", "remove_files", "make", "make_install", "download", "gitclone", - "record_mfem_sha", "cmake", + "record_mfem_sha", "cmake", "get_numpy_inc", "get_mpi4py_inc", "find_libpath_from_prefix", "clean_so", ] @@ -102,27 +102,26 @@ def make_call(command, target='', force_verbose=False, env=None): env.update(os.environ) myverbose = bglb.verbose or force_verbose - + kwargs['stdout'] = subprocess.PIPE kwargs['stderr'] = subprocess.PIPE - + p = subprocess.Popen(command, **kwargs) - stdout, stderr= p.communicate() + stdout, stderr = p.communicate() if p.returncode != 0: if target == '': target = " ".join(command) - print("!!!!!!!!!!!!!!!", p) - print(stdout) - print(stderr) - + print("!!!!!!!!!!!!!!! make call faield !!!!!!!!!!!!!!!!!") print("Failed when calling command: " + target) + print(stdout) + print(stderr) raise subprocess.CalledProcessError(p.returncode, " ".join(command)) - if not myverbose: print(stdout) + def chdir(path): ''' change directory to `path`; returns the previous directory @@ -253,7 +252,6 @@ def cmake(path, **kwargs): make_call(command) - def get_numpy_inc(): python = sys.executable @@ -314,6 +312,7 @@ def find_libpath_from_prefix(lib, prefix0): return '' + def clean_so(all=None): python = sys.executable command = [python, "setup.py", "clean"] @@ -337,6 +336,3 @@ def clean_so(all=None): make_call(command) chdir(pwd) - - - diff --git a/mfem/_ser/bilininteg_hdg.i b/mfem/_ser/bilininteg_hdg.i index fed4ce1c..708de0d2 100644 --- a/mfem/_ser/bilininteg_hdg.i +++ b/mfem/_ser/bilininteg_hdg.i @@ -8,6 +8,9 @@ #include "numpy/arrayobject.h" #include "../common/pyoperator.hpp" #include "../common/pyintrules.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pynonlininteg.hpp" %} %include "../common/existing_mfem_headers.i" @@ -19,7 +22,8 @@ import_array(); %include "exception.i" %include "../common/exception.i" +%import "bilininteg.i" + %include "fem/darcy/bilininteg_hdg.hpp" #endif - diff --git a/mfem/_ser/darcyform.i b/mfem/_ser/darcyform.i index c1c4bf7d..b0e33c2c 100644 --- a/mfem/_ser/darcyform.i +++ b/mfem/_ser/darcyform.i @@ -8,6 +8,8 @@ #include "numpy/arrayobject.h" #include "../common/pyoperator.hpp" #include "../common/pyintrules.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pybilininteg.hpp" %} %include "../common/existing_mfem_headers.i" @@ -23,6 +25,8 @@ import_array(); %import "blockvector.i" %import "bilinearform.i" %import "nonlinearform.i" +%import "darcyreduction.i" +%import "darcyhybridization.i" %include "fem/darcy/darcyform.hpp" diff --git a/mfem/_ser/darcyhybridization.i b/mfem/_ser/darcyhybridization.i index d95ad474..8f628c9f 100644 --- a/mfem/_ser/darcyhybridization.i +++ b/mfem/_ser/darcyhybridization.i @@ -8,6 +8,7 @@ #include "numpy/arrayobject.h" #include "../common/pyoperator.hpp" #include "../common/pyintrules.hpp" +#include "../common/pycoefficient.hpp" %} %include "../common/existing_mfem_headers.i" @@ -19,7 +20,10 @@ import_array(); %include "exception.i" %include "../common/exception.i" + //%ignore mfem::DarcyHybridization::m_nlfi_u; +%ignore mfem::DarcyHybridization; +%import "nonlininteg.i" + %include "fem/darcy/darcyhybridization.hpp" #endif - diff --git a/mfem/_ser/darcyreduction.i b/mfem/_ser/darcyreduction.i index a3a1c4c2..52d77761 100644 --- a/mfem/_ser/darcyreduction.i +++ b/mfem/_ser/darcyreduction.i @@ -8,6 +8,9 @@ #include "numpy/arrayobject.h" #include "../common/pyoperator.hpp" #include "../common/pyintrules.hpp" +#include "../common/pycoefficient.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pynonlininteg.hpp" %} %include "../common/existing_mfem_headers.i" @@ -19,7 +22,12 @@ import_array(); %include "exception.i" %include "../common/exception.i" +%import "fespace.i" +%import "nonlinearform.i" +%import "sparsemat.i" +%import "blockmatrix.i" +%import "blockvector.i" + %include "fem/darcy/darcyreduction.hpp" #endif - diff --git a/mfem/_ser/estimators_hdg.i b/mfem/_ser/estimators_hdg.i index 476c3826..f910245f 100644 --- a/mfem/_ser/estimators_hdg.i +++ b/mfem/_ser/estimators_hdg.i @@ -8,6 +8,8 @@ #include "numpy/arrayobject.h" #include "../common/pyoperator.hpp" #include "../common/pyintrules.hpp" +#include "../common/pybilininteg.hpp" +#include "../common/pycoefficient.hpp" %} %include "../common/existing_mfem_headers.i" @@ -19,7 +21,8 @@ import_array(); %include "exception.i" %include "../common/exception.i" +%import "bilininteg.i" + %include "fem/darcy/estimators_hdg.hpp" #endif - diff --git a/mfem/_ser/setup.py b/mfem/_ser/setup.py index 0a062e21..b3d0d23a 100644 --- a/mfem/_ser/setup.py +++ b/mfem/_ser/setup.py @@ -112,7 +112,7 @@ def get_extensions(): "complex_densemat", "complexstaticcond", "complexweakform", "bounds", "integrator", "darcyform","estimators_hdg", "bilininteg_hdg", - "darcyreductoin", "darcyhybridization"] + "darcyreduction", "darcyhybridization"] if add_cuda == '1': from setup_local import cudainc From 84e6e63bc143663479ec596e3e8b760b26b00e7b Mon Sep 17 00:00:00 2001 From: Syun'ichi Shiraiwa Date: Fri, 21 Nov 2025 11:50:26 -0500 Subject: [PATCH 3/6] fixed typo in build system --- _build_system/build_pymfem.py | 21 +++++++++++++-------- _build_system/build_utils.py | 2 +- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/_build_system/build_pymfem.py b/_build_system/build_pymfem.py index 809303f5..41406326 100644 --- a/_build_system/build_pymfem.py +++ b/_build_system/build_pymfem.py @@ -147,7 +147,19 @@ def check_new(ifile): return os.path.getmtime(ifile) > os.path.getmtime(wfile) def update_integrator_exts(): - pwd = chdir(os.path.join(rootdir, 'mfem', 'common')) + commdir = os.path.join(rootdir, 'mfem', 'common') + pwd = chdir(commdir) + + for filename in ('bilininteg_ext.i', 'lininteg_ext.i'): + fid = open(filename, "w") + fid.close() + + os.chdir(pwd) + for filename in ['lininteg.i', 'bilininteg.i']: + command = [swig_command] + swigflag + serflag + [filename] + make_call(command) + + os.chdir(commdir) command1 = [sys.executable, "generate_lininteg_ext.py"] command2 = [sys.executable, "generate_bilininteg_ext.py"] make_call(command1) @@ -196,13 +208,6 @@ def update_header_exists(mfem_source): # create empty files # run swig with *.i files # call update scripts - for filename in ('../common/bilininteg_ext.i', - '../common/lininteg_ext.i'): - fid = open(filename) - fid.close() - for filename in ['lininteg.i', 'bilininteg.i']: - command = [swig_command] + swigflag + serflag + [filename] - make_call(command) update_integrator_exts() commands = [] diff --git a/_build_system/build_utils.py b/_build_system/build_utils.py index 2bb98b37..30821667 100644 --- a/_build_system/build_utils.py +++ b/_build_system/build_utils.py @@ -118,7 +118,7 @@ def make_call(command, target='', force_verbose=False, env=None): print(stderr) raise subprocess.CalledProcessError(p.returncode, " ".join(command)) - if not myverbose: + if myverbose: print(stdout) From a0ab96fea5b5317840aa0df36fa0afd638b42700 Mon Sep 17 00:00:00 2001 From: Syun'ichi Shiraiwa Date: Fri, 21 Nov 2025 12:17:00 -0500 Subject: [PATCH 4/6] added error checking of command-line options --- _build_system/build_config.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/_build_system/build_config.py b/_build_system/build_config.py index 17595562..21e6d321 100644 --- a/_build_system/build_config.py +++ b/_build_system/build_config.py @@ -230,6 +230,9 @@ def process_cmd_options(command_obj, cfs): else: setattr(command_obj, attr, False) + if len(cfs) != 0: + assert False, "Following build options are undefined: " + ", ".join(list(cfs)) + def process_setup_options(command_obj, args): for item in args: From 6de722f6d90962a2915912b415d98d3d26a5846e Mon Sep 17 00:00:00 2001 From: Syun'ichi Shiraiwa Date: Fri, 21 Nov 2025 13:31:43 -0500 Subject: [PATCH 5/6] (WIP) working on ex5-hdg.py --- examples/ex5-hdg.py | 445 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 445 insertions(+) create mode 100644 examples/ex5-hdg.py diff --git a/examples/ex5-hdg.py b/examples/ex5-hdg.py new file mode 100644 index 00000000..3e01143a --- /dev/null +++ b/examples/ex5-hdg.py @@ -0,0 +1,445 @@ +''' + MFEM example 5 + + See c++ version in the MFEM library for more detail +''' +import os +import time +import mfem.ser as mfem +from mfem.ser import intArray +from os.path import expanduser, join, dirname +import numpy as np +from numpy import sin, cos, exp + +def run(order=1, + meshfile='', + nx=0, + ny=0, + dg=False, + brt=False, + td=False, + hb=False, + rd=False, + sn=False, + visualization=False, + device='cpu', + numba=False, + pa=False): + + def pFunc_exact(x): + xi = float(x[0]) + yi = float(x[1]) + zi = 0.0 + if len(x) == 3: + zi = x[2] + from numpy import sin, cos, exp + return exp(xi)*sin(yi)*cos(zi) + + class uFunc_ex(mfem.VectorPyCoefficient): + def EvalValue(self, x): + xi = float(x[0]) + yi = float(x[1]) + zi = 0.0 + if len(x) == 3: + zi = x[2] + ret = [- exp(xi)*sin(yi)*cos(zi), + - exp(xi)*cos(yi)*cos(zi)] + if len(x) == 3: + ret.append(exp(xi)*sin(yi)*sin(zi)) + return ret + + class pFunc_ex(mfem.PyCoefficient): + def EvalValue(self, x): + return pFunc_exact(x) + + class fFunc(mfem.VectorPyCoefficient): + def EvalValue(self, x): + if len(x) == 3: + return [0., 0., 0.] + else: + return [0., 0.] + + class gFunc(mfem.PyCoefficient): + def EvalValue(self, x): + if len(x) == 3: + return -pFunc_exact(x) + return 0. + + class f_natural(mfem.PyCoefficient): + def EvalValue(self, x): + return -pFunc_exact(x) + + device = mfem.Device(device) + device.Print() + + if meshfile != "": + mesh = mfem.Mesh(meshfile, 1, 1) + else: + if ny <= 0: + ny = nx + mesh = mfem.Mesh.MakeCartesian2D(nx, ny, mfem.Element.QUADRILATERL) + + dim = mesh.Dimension() + + # 4. Refine the mesh to increase the resolution. In this example we do + # 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the + # largest number that gives a final mesh with no more than 10,000 + # elements. + + ref_levels = int(np.floor(np.log(10000./mesh.GetNE())/np.log(2.)/dim)) + for x in range(ref_levels): + mesh.UniformRefinement() + + # 5. Define a finite element space on the mesh. Here we use the + # Raviart-Thomas finite elements of the specified order. + if dg: + R_coll = mfem.L2_FECollection(order, dim, mfem.BasisType.GaussLobatto) + elif brt: + R_coll = mfem.BrokenRT_FECollection(order, dim) + R_coll_dg = mfem.L2_FECollection(order+1, dim) + else: + R_coll = mfem.RT_FECollection(order, dim) + + W_coll = mfem.L2_FECollection(order, dim) + + R_space = mfem.FiniteElementSpace(mesh, R_coll, dim if dg else 1) + if brt: + R_space_hg = mfem.FiniteElementSpace(mesh, R_coll_dg, dim) + W_space = mfem.FiniteElementSpace(mesh, W_coll) + + darcy = mfem.DarcyForm(R_space, W_space) + + # 6. Define the BlockStructure of the problem, i.e. define the array of + # offsets for each variable. The last component of the Array is the sum + # of the dimensions of each block. + + block_offsets = darcy.GetOffsets().ToList() + + print("***********************************************************") + print("dim(R) = " + str(block_offsets[1] - block_offsets[0])) + print("dim(W) = " + str(block_offsets[2] - block_offsets[1])) + print("dim(R+W) = " + str(block_offsets[-1])) + print("***********************************************************") + + # 7. Define the coefficients, analytical solution, and rhs of the PDE. + + kcoeff = mfem.ConstantCoefficient(1.0) # acoustic resistance + ikcoeff = mfem.RatioCoefficient(1., kcoeff) # inverse acoustic resistance + + fcoeff = fFunc(dim) + fnatcoeff = f_natural() + gcoeff = gFunc() + ucoeff = uFunc_ex(dim) + pcoeff = pFunc_ex() + + # 8. Allocate memory (x, rhs) for the analytical solution and the right hand + # side. Define the GridFunction u,p for the finite element solution and + # linear forms fform and gform for the right hand side. The data + # allocated by x and rhs are passed as a reference to the grid functions + # (u,p) and the linear forms (fform, gform). + + mt = device.GetMemoryType(); + x = mfem.BlockVector(block_offsets, mt) + rhs = mfem.BlockVector(block_offsets, mt) + + fform = mfem.LinearForm() + fform.Update(R_space, rhs.GetBlock(0), 0) + + if dg: + fform.AddDomainIntegrator(mfem.VectorDomainLFIntegrator(fcoeff)) + fform.AddBdrFaceIntegrator(mfem.VectorBoundaryFluxLFIntegrator(fnatcoeff)) + else + { + fform.AddDomainIntegrator(mfem.VectorFEDomainLFIntegrator(fcoeff)) + if brt: + fform.AddBdrFaceIntegrator(mfem.VectorFEBoundaryFluxLFIntegrator(fnatcoeff)) + else: + fform.AddBoundaryIntegrator(mfem.VectorFEBoundaryFluxLFIntegrator(fnatcoeff)) + + + fform.Assemble() + fform.SyncAliasMemory(rhs) + + gform = mfem.LinearForm() + gform.Update(W_space, rhs.GetBlock(1), 0) + gform.AddDomainIntegrator(mfem.DomainLFIntegrator(gcoeff)) + gform.Assemble() + gform.SyncAliasMemory(rhs) + + # 9. Assemble the finite element matrices for the Darcy operator + # D = [ M B^T ] + # [ B 0 ] + # where: + # + # M = \int_\Omega k u_h \cdot v_h d\Omega u_h, v_h \in R_h + # B = -\int_\Omega \div u_h q_h d\Omega u_h \in R_h, q_h \in W_h + mVarf = darcy.GetFluxMassForm() + bVarf = darcy.GetFluxDivForm() + + if dg: + mtVarf = darcy.GetPotentialMassForm() + mVarf.AddDomainIntegrator(mfem.VectorMassIntegrator(kcoeff)) + bVarf.AddDomainIntegrator(mfem.VectorDivergenceIntegrator()) + bVarf.AddInteriorFaceIntegrator(mfem.TransposeIntegrator( + mfem.DGNormalTraceIntegrator(-1.))) + mtVarf.AddInteriorFaceIntegrator(mfem.HDGDiffusionIntegrator(ikcoeff, td)) + else: + mVarf.AddDomainIntegrator(mfem.VectorFEMassIntegrator(kcoeff)) + bVarf.AddDomainIntegrator(mfem.VectorFEDivergenceIntegrator) + if brt: + bVarf.AddInteriorFaceIntegrator(mfem.TransposeIntegrator( + mfem.DGNormalTraceIntegrator(-1.))) + + #set hybridization / assembly level + chrono = time.time() + + ess_flux_tdofs_list = mfem.intArray() + if hybridization: + trace_coll = mfem.DG_Interface_FECollection(order, dim); + trace_space = mfem.FiniteElementSpace(mesh, trace_coll); + darcy.EnableHybridization(trace_space, + mfem.NormalTraceJumpIntegrator(), + ess_flux_tdofs_list) + elif (reduction and (dg or brt)): + darcy.EnableFluxReduction() + + if pa: + mVarf.SetAssemblyLevel(mfem.AssemblyLevel_PARTIAL) + + darcy.Assemble() + + pDarcyOp = mfem.OperatorHandle() + X = mfem.Vector() + B = mfem.Vector() + x.Assign(0.0) + darcy.FormLinearSystem(ess_flux_tdofs_list, x, rhs, pDarcyOp, X, B) + + print("Assembly took " + str(time.time() - chrono)) + + maxIter = 1000 + rtol = 1e-6 + atol = 1e-10 + + if (hybridization or (reduction and (dg or brt))): + # 10. Construct the preconditioner + prec = mfem.GSSmoother(pDarcyOp.AsSparseMatrix()) + + # 11. Solve the linear system with GMRES. + # Check the norm of the unpreconditioned residual. + chrono = time.time() + solver = mfem.GMRESSolver() + solver.SetAbsTol(atol) + solver.SetRelTol(rtol) + solver.SetMaxIter(maxIter) + solver.SetOperator(pDarcyOp) + solver.SetPreconditioner(prec) + solver.SetPrintLevel(1) + + p + + + mVarf.AddDomainIntegrator(mfem.VectorFEMassIntegrator(k)) + mVarf.Assemble() + if not pa: + mVarf.Finalize() + + if pa: + bVarf.SetAssemblyLevel(mfem.AssemblyLevel_PARTIAL) + bVarf.AddDomainIntegrator(mfem.VectorFEDivergenceIntegrator()) + bVarf.Assemble() + if not pa: + bVarf.Finalize() + + darcyOp = mfem.BlockOperator(block_offsets) + + if pa: + Bt = mfem.TransposeOperator(bVarf) + + darcyOp.SetBlock(0, 0, mVarf) + darcyOp.SetBlock(0, 1, Bt, -1) + darcyOp.SetBlock(1, 0, bVarf, -1) + else: + M = mVarf.SpMat() + B = bVarf.SpMat() + B *= -1 + if mfem.Device.IsEnabled(): + B = mfem.Transpose(B) + Bt = mfem.TransposeOperator(B) + + darcyOp.SetBlock(0, 0, M) + darcyOp.SetBlock(0, 1, Bt) + darcyOp.SetBlock(1, 0, B) + + Md = mfem.Vector(mVarf.Height()) + darcyPrec = mfem.BlockDiagonalPreconditioner(block_offsets) + + if pa: + mVarf.AssembleDiagonal(Md) + Md_host = mfem.Vector(Md.HostRead(), mVarf.Height()) + invMd = mfem.Vector(mVarf.Height()) + + for i in range(invMd.Size()): + invMd[i] = 1.0/Md_host[i] + + BMBt_diag = mfem.Vector(bVarf.Height()) + bVarf.AssembleDiagonal_ADAt(invMd, BMBt_diag) + ess_tdof_list = mfem.intArray() + invM = mfem.OperatorJacobiSmoother(Md, ess_tdof_list) + invS = mfem.OperatorJacobiSmoother(BMBt_diag, ess_tdof_list) + else: + MinvBt = mfem.Transpose(B) + M.GetDiag(Md) + + for i in range(Md.Size()): + MinvBt.ScaleRow(i, 1/Md[i]) + S = mfem.Mult(B, MinvBt) + + invM = mfem.DSmoother(M) + invS = mfem.GSSmoother(S) + + invM.iterative_mode = False + invS.iterative_mode = False + + darcyPrec.SetDiagonalBlock(0, invM) + darcyPrec.SetDiagonalBlock(1, invS) + + maxIter = 1000 + rtol = 1e-6 + atol = 1e-10 + + solver = mfem.MINRESSolver() + solver.SetAbsTol(atol) + solver.SetRelTol(rtol) + solver.SetMaxIter(maxIter) + solver.SetOperator(darcyOp) + solver.SetPreconditioner(darcyPrec) + solver.SetPrintLevel(1) + x.Assign(0.0) + solver.Mult(rhs, x) + + solve_time = clock() - stime + + if solver.GetConverged(): + print("MINRES converged in " + str(solver.GetNumIterations()) + + " iterations with a residual norm of " + "{:g}".format(solver.GetFinalNorm())) + else: + print("MINRES did not converge in " + str(solver.GetNumIterations()) + + " iterations. Residual norm is " + "{:g}".format(solver.GetFinalNorm())) + print("MINRES solver took " + str(solve_time) + "s.") + + u = mfem.GridFunction() + p = mfem.GridFunction() + u.MakeRef(R_space, x.GetBlock(0), 0) + p.MakeRef(W_space, x.GetBlock(1), 0) + + order_quad = max(2, 2*order+1) + + irs = [mfem.IntRules.Get(i, order_quad) + for i in range(mfem.Geometry.NumGeom)] + + norm_p = mfem.ComputeLpNorm(2, pcoeff, mesh, irs) + norm_u = mfem.ComputeLpNorm(2, ucoeff, mesh, irs) + err_u = u.ComputeL2Error(ucoeff, irs) + err_p = p.ComputeL2Error(pcoeff, irs) + + print("|| u_h - u_ex || / || u_ex || = " + "{:g}".format(err_u / norm_u)) + print("|| p_h - p_ex || / || p_ex || = " + "{:g}".format(err_p / norm_p)) + + mesh.Print("ex5.mesh", 8) + u.Save("sol_u.gf", 8) + p.Save("sol_p.gf", 8) + + visit_dc = mfem.VisItDataCollection("Example5", mesh) + visit_dc.RegisterField("velocity", u) + visit_dc.RegisterField("pressure", p) + visit_dc.Save() + + paraview_dc = mfem.ParaViewDataCollection("Example5", mesh) + paraview_dc.SetPrefixPath("ParaView") + paraview_dc.SetLevelsOfDetail(order) + paraview_dc.SetCycle(0) + paraview_dc.SetDataFormat(mfem.VTKFormat_BINARY) + paraview_dc.SetHighOrderOutput(True) + paraview_dc.SetTime(0.0) + paraview_dc.RegisterField("velocity", u) + paraview_dc.RegisterField("pressure", p) + paraview_dc.Save() + + if visualization: + sout_u = mfem.socketstream("localhost", 19916) + sout_u.precision(8) + sout_u << "solution\n" << mesh << u << "window_title 'Velocity'\n" + sout_p = mfem.socketstream("localhost", 19916) + sout_p.precision(8) + sout_p << "solution\n" << mesh << p << "window_title 'Pressure'\n" + + +if __name__ == "__main__": + from mfem.common.arg_parser import ArgParser + + parser = ArgParser(description='Ex5 (Darcy Problem)') + parser.add_argument('-m', '--mesh', + default="", + action='store', type=str, + help='Mesh file to use.') + parser.add_argument("-nx", "--ncells-x", + action='store', type=int, default=0, + help="Number of cells in x.") + parser.add_argument("-ny", "--ncells-y", + action='store', type=int, default=0, + help="Number of cells in y.") + parser.add_argument('-o', '--order', + action='store', default=1, type=int, + help="Finite element order (polynomial degree)."); + parser.add_argument("-dg", "--discontinuous", + action='store_true', + help="Enable DG elements for fluxes."); + parser.add_argument("-brt", "--broken-RT", + action='store_true', + help="Enable broken RT elements for fluxes.") + parser.add_argument("-td", "--stab_diff", + action='store', type=float, default=0.5, + help="Diffusion stabilization factor (1/2=default)") + parser.add_argument("-hb", "--hybridization", + action='store_true', + help="Enable hybridization.") + parser.add_argument("-rd", "--reduction", + action='store_true', + help="Enable reduction of DG flux.") + parser.add_argument("-sn", "--solution-norm", + action='store_true', + help="Solution norm (0=native, 1=flux, 2=potential).") + 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().") + parser.add_argument('-vis', '--visualization', + action='store_true', + help='Enable GLVis visualization') + + + args = parser.parse_args() + parser.print_options(args) + + order = args.order + meshfile = expanduser( + join(dirname(__file__), '..', 'data', args.mesh)) + visualization = args.visualization + device = args.device + pa = args.partial_assembly + + run(order=order, + meshfile=meshfile, + nx=args.ncells_x, + ny=args.ncells_y, + dg=args.discontinous, + brt=args.broken_RF, + td=args.stab_diff, + hb=args.hybridizatoin, + rd=args.reduciton, + sn=args.solution_norm, + visualization=visualization, + device=device, + pa=pa) From 14e16787e4fcd91c02a04237d1de44e7a67fabd2 Mon Sep 17 00:00:00 2001 From: Syun'ichi Shiraiwa Date: Fri, 21 Nov 2025 15:59:27 -0500 Subject: [PATCH 6/6] ex5-hdg.py runs w/o crashing --- examples/ex5-hdg.py | 303 +++++++++++++++++++++++--------------------- mfem/ser.py | 2 + 2 files changed, 163 insertions(+), 142 deletions(-) diff --git a/examples/ex5-hdg.py b/examples/ex5-hdg.py index 3e01143a..4e7e3690 100644 --- a/examples/ex5-hdg.py +++ b/examples/ex5-hdg.py @@ -6,11 +6,15 @@ import os import time import mfem.ser as mfem + from mfem.ser import intArray from os.path import expanduser, join, dirname import numpy as np from numpy import sin, cos, exp +use_umfpack = hasattr(mfem, "UMFPackSolver") + + def run(order=1, meshfile='', nx=0, @@ -77,7 +81,8 @@ def EvalValue(self, x): else: if ny <= 0: ny = nx - mesh = mfem.Mesh.MakeCartesian2D(nx, ny, mfem.Element.QUADRILATERL) + print(nx, ny) + mesh = mfem.Mesh.MakeCartesian2D(nx, ny, mfem.Element.QUADRILATERAL) dim = mesh.Dimension() @@ -85,7 +90,7 @@ def EvalValue(self, x): # 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the # largest number that gives a final mesh with no more than 10,000 # elements. - + ref_levels = int(np.floor(np.log(10000./mesh.GetNE())/np.log(2.)/dim)) for x in range(ref_levels): mesh.UniformRefinement() @@ -99,9 +104,9 @@ def EvalValue(self, x): R_coll_dg = mfem.L2_FECollection(order+1, dim) else: R_coll = mfem.RT_FECollection(order, dim) - + W_coll = mfem.L2_FECollection(order, dim) - + R_space = mfem.FiniteElementSpace(mesh, R_coll, dim if dg else 1) if brt: R_space_hg = mfem.FiniteElementSpace(mesh, R_coll_dg, dim) @@ -113,18 +118,18 @@ def EvalValue(self, x): # offsets for each variable. The last component of the Array is the sum # of the dimensions of each block. - block_offsets = darcy.GetOffsets().ToList() - + block_offsets = darcy.GetOffsets() + offsets = block_offsets.ToList() print("***********************************************************") - print("dim(R) = " + str(block_offsets[1] - block_offsets[0])) - print("dim(W) = " + str(block_offsets[2] - block_offsets[1])) - print("dim(R+W) = " + str(block_offsets[-1])) + print("dim(R) = " + str(offsets[1] - offsets[0])) + print("dim(W) = " + str(offsets[2] - offsets[1])) + print("dim(R+W) = " + str(offsets[-1])) print("***********************************************************") # 7. Define the coefficients, analytical solution, and rhs of the PDE. - kcoeff = mfem.ConstantCoefficient(1.0) # acoustic resistance - ikcoeff = mfem.RatioCoefficient(1., kcoeff) # inverse acoustic resistance + kcoeff = mfem.ConstantCoefficient(1.0) # acoustic resistance + ikcoeff = mfem.RatioCoefficient(1., kcoeff) # inverse acoustic resistance fcoeff = fFunc(dim) fnatcoeff = f_natural() @@ -137,35 +142,36 @@ def EvalValue(self, x): # linear forms fform and gform for the right hand side. The data # allocated by x and rhs are passed as a reference to the grid functions # (u,p) and the linear forms (fform, gform). - - mt = device.GetMemoryType(); + + mt = mfem.MemoryType(device.GetDeviceMemoryType()) x = mfem.BlockVector(block_offsets, mt) rhs = mfem.BlockVector(block_offsets, mt) - + fform = mfem.LinearForm() fform.Update(R_space, rhs.GetBlock(0), 0) if dg: fform.AddDomainIntegrator(mfem.VectorDomainLFIntegrator(fcoeff)) - fform.AddBdrFaceIntegrator(mfem.VectorBoundaryFluxLFIntegrator(fnatcoeff)) - else - { - fform.AddDomainIntegrator(mfem.VectorFEDomainLFIntegrator(fcoeff)) - if brt: - fform.AddBdrFaceIntegrator(mfem.VectorFEBoundaryFluxLFIntegrator(fnatcoeff)) - else: - fform.AddBoundaryIntegrator(mfem.VectorFEBoundaryFluxLFIntegrator(fnatcoeff)) - + fform.AddBdrFaceIntegrator( + mfem.VectorBoundaryFluxLFIntegrator(fnatcoeff)) + else: + fform.AddDomainIntegrator(mfem.VectorFEDomainLFIntegrator(fcoeff)) + if brt: + fform.AddBdrFaceIntegrator( + mfem.VectorFEBoundaryFluxLFIntegrator(fnatcoeff)) + else: + fform.AddBoundaryIntegrator( + mfem.VectorFEBoundaryFluxLFIntegrator(fnatcoeff)) fform.Assemble() fform.SyncAliasMemory(rhs) - + gform = mfem.LinearForm() gform.Update(W_space, rhs.GetBlock(1), 0) gform.AddDomainIntegrator(mfem.DomainLFIntegrator(gcoeff)) gform.Assemble() gform.SyncAliasMemory(rhs) - + # 9. Assemble the finite element matrices for the Darcy operator # D = [ M B^T ] # [ B 0 ] @@ -177,32 +183,33 @@ def EvalValue(self, x): bVarf = darcy.GetFluxDivForm() if dg: - mtVarf = darcy.GetPotentialMassForm() - mVarf.AddDomainIntegrator(mfem.VectorMassIntegrator(kcoeff)) - bVarf.AddDomainIntegrator(mfem.VectorDivergenceIntegrator()) - bVarf.AddInteriorFaceIntegrator(mfem.TransposeIntegrator( - mfem.DGNormalTraceIntegrator(-1.))) - mtVarf.AddInteriorFaceIntegrator(mfem.HDGDiffusionIntegrator(ikcoeff, td)) + mtVarf = darcy.GetPotentialMassForm() + mVarf.AddDomainIntegrator(mfem.VectorMassIntegrator(kcoeff)) + bVarf.AddDomainIntegrator(mfem.VectorDivergenceIntegrator()) + bVarf.AddInteriorFaceIntegrator(mfem.TransposeIntegrator( + mfem.DGNormalTraceIntegrator(-1.))) + mtVarf.AddInteriorFaceIntegrator( + mfem.HDGDiffusionIntegrator(ikcoeff, td)) else: mVarf.AddDomainIntegrator(mfem.VectorFEMassIntegrator(kcoeff)) - bVarf.AddDomainIntegrator(mfem.VectorFEDivergenceIntegrator) - if brt: - bVarf.AddInteriorFaceIntegrator(mfem.TransposeIntegrator( - mfem.DGNormalTraceIntegrator(-1.))) + bVarf.AddDomainIntegrator(mfem.VectorFEDivergenceIntegrator()) + if brt: + bVarf.AddInteriorFaceIntegrator(mfem.TransposeIntegrator( + mfem.DGNormalTraceIntegrator(-1.))) - #set hybridization / assembly level + # set hybridization / assembly level chrono = time.time() - + ess_flux_tdofs_list = mfem.intArray() - if hybridization: - trace_coll = mfem.DG_Interface_FECollection(order, dim); - trace_space = mfem.FiniteElementSpace(mesh, trace_coll); - darcy.EnableHybridization(trace_space, - mfem.NormalTraceJumpIntegrator(), - ess_flux_tdofs_list) - elif (reduction and (dg or brt)): - darcy.EnableFluxReduction() - + if hb: + trace_coll = mfem.DG_Interface_FECollection(order, dim) + trace_space = mfem.FiniteElementSpace(mesh, trace_coll) + darcy.EnableHybridization(trace_space, + mfem.NormalTraceJumpIntegrator(), + ess_flux_tdofs_list) + elif (rd and (dg or brt)): + darcy.EnableFluxReduction() + if pa: mVarf.SetAssemblyLevel(mfem.AssemblyLevel_PARTIAL) @@ -210,92 +217,91 @@ def EvalValue(self, x): pDarcyOp = mfem.OperatorHandle() X = mfem.Vector() - B = mfem.Vector() + B = mfem.Vector() x.Assign(0.0) darcy.FormLinearSystem(ess_flux_tdofs_list, x, rhs, pDarcyOp, X, B) print("Assembly took " + str(time.time() - chrono)) - + maxIter = 1000 rtol = 1e-6 atol = 1e-10 - if (hybridization or (reduction and (dg or brt))): - # 10. Construct the preconditioner - prec = mfem.GSSmoother(pDarcyOp.AsSparseMatrix()) - - # 11. Solve the linear system with GMRES. - # Check the norm of the unpreconditioned residual. - chrono = time.time() - solver = mfem.GMRESSolver() - solver.SetAbsTol(atol) - solver.SetRelTol(rtol) - solver.SetMaxIter(maxIter) - solver.SetOperator(pDarcyOp) - solver.SetPreconditioner(prec) - solver.SetPrintLevel(1) - - p - - - mVarf.AddDomainIntegrator(mfem.VectorFEMassIntegrator(k)) - mVarf.Assemble() - if not pa: - mVarf.Finalize() - - if pa: - bVarf.SetAssemblyLevel(mfem.AssemblyLevel_PARTIAL) - bVarf.AddDomainIntegrator(mfem.VectorFEDivergenceIntegrator()) - bVarf.Assemble() - if not pa: - bVarf.Finalize() - - darcyOp = mfem.BlockOperator(block_offsets) - - if pa: - Bt = mfem.TransposeOperator(bVarf) - - darcyOp.SetBlock(0, 0, mVarf) - darcyOp.SetBlock(0, 1, Bt, -1) - darcyOp.SetBlock(1, 0, bVarf, -1) - else: - M = mVarf.SpMat() - B = bVarf.SpMat() - B *= -1 - if mfem.Device.IsEnabled(): - B = mfem.Transpose(B) - Bt = mfem.TransposeOperator(B) - - darcyOp.SetBlock(0, 0, M) - darcyOp.SetBlock(0, 1, Bt) - darcyOp.SetBlock(1, 0, B) - - Md = mfem.Vector(mVarf.Height()) - darcyPrec = mfem.BlockDiagonalPreconditioner(block_offsets) + if hb or (rd and (dg or brt)): + # 10. Construct the preconditioner + prec = mfem.GSSmoother(pDarcyOp.AsSparseMatrix()) + + # 11. Solve the linear system with GMRES. + # Check the norm of the unpreconditioned residual. + chrono = time.time() + solver = mfem.GMRESSolver() + solver.SetAbsTol(atol) + solver.SetRelTol(rtol) + solver.SetMaxIter(maxIter) + solver.SetOperator(pDarcyOp) + solver.SetPreconditioner(prec) + solver.SetPrintLevel(1) + + solver.Mult(B, X) + darcy.RecoverFEMSolution(X, rhs, x) + + if solver.GetConverged(): + print("GMRES converged in " + str(solver.GetNumIterations()) + + " iterations with a residual norm of " + str(solver.GetFinalNorm())) + else: + print("GMRES did not converge in " + str(solver.GetNumIterations()) + + " iterations with a residual norm of " + str(solver.GetFinalNorm())) + + print("GMRES solver took " + str(time.time() - chrono) + "s.\n") - if pa: - mVarf.AssembleDiagonal(Md) - Md_host = mfem.Vector(Md.HostRead(), mVarf.Height()) - invMd = mfem.Vector(mVarf.Height()) - - for i in range(invMd.Size()): - invMd[i] = 1.0/Md_host[i] - - BMBt_diag = mfem.Vector(bVarf.Height()) - bVarf.AssembleDiagonal_ADAt(invMd, BMBt_diag) - ess_tdof_list = mfem.intArray() - invM = mfem.OperatorJacobiSmoother(Md, ess_tdof_list) - invS = mfem.OperatorJacobiSmoother(BMBt_diag, ess_tdof_list) else: - MinvBt = mfem.Transpose(B) - M.GetDiag(Md) - - for i in range(Md.Size()): - MinvBt.ScaleRow(i, 1/Md[i]) - S = mfem.Mult(B, MinvBt) - - invM = mfem.DSmoother(M) - invS = mfem.GSSmoother(S) + # 10. Construct the operators for preconditioner + # + # P = [ diag(M) 0 ] + # [ 0 B diag(M)^-1 B^T ] + # + # Here we use Symmetric Gauss-Seidel to approximate the inverse of the + # pressure Schur Complement + MinvBt = mfem.SparseMatrix() + Md = mfem.Vector() + darcyPrec = mfem.BlockDiagonalPreconditioner(block_offsets) + + if pa: + mVarf.AssembleDiagonal(Md) + Md_host = Md.HostRead() + intMd = mfem.Vector(mVarf.Height()) + + intMd.GetDataArray()[:] = 1./Md_host.GetDataArray() + + BMBt_diag = mfem.Vector(bVarf.Height()) + bVarf.AssembleDiagonal_ADAt(invMd, BMBt_diag) + + ess_tdof_list = intArray() + + invM = mfem.OperatorJacobiSmoother(Md, ess_tdof_list) + invS = mfem.OperatorJacobiSmoother(BMBt_diag, ess_tdof_list) + else: + M = mfem.SparseMatrix(mVarf.SpMat()) + M.GetDiag(Md) + Md.HostReadWrite() + + B = mfem.SparseMatrix(bVarf.SpMat()) + MinvBt = mfem.Transpose(B) + + for i in range(Md.Size()): + MinvBt.ScaleRow(i, 1./Md[i]) + + S = mfem.Mult(B, MinvBt) + if dg: + Mtm = mfem.SparseMatrix(mtVarf.SpMat()) + Snew = mfem.Add(Mtm, S) + S = Snew + invM = mfem.DSmoother(M) + + if use_umfpack: + invS = mfem.UMFPackSolver(S) + else: + invS = mfem.GSSmoother(S) invM.iterative_mode = False invS.iterative_mode = False @@ -303,21 +309,18 @@ def EvalValue(self, x): darcyPrec.SetDiagonalBlock(0, invM) darcyPrec.SetDiagonalBlock(1, invS) - maxIter = 1000 - rtol = 1e-6 - atol = 1e-10 + chrono = time.time() solver = mfem.MINRESSolver() solver.SetAbsTol(atol) solver.SetRelTol(rtol) solver.SetMaxIter(maxIter) - solver.SetOperator(darcyOp) + solver.SetOperator(pDarcyOp) solver.SetPreconditioner(darcyPrec) solver.SetPrintLevel(1) - x.Assign(0.0) - solver.Mult(rhs, x) - solve_time = clock() - stime + solver.Mult(rhs, x) + solve_time = time.time() - chrono if solver.GetConverged(): print("MINRES converged in " + str(solver.GetNumIterations()) + @@ -325,11 +328,20 @@ def EvalValue(self, x): else: print("MINRES did not converge in " + str(solver.GetNumIterations()) + " iterations. Residual norm is " + "{:g}".format(solver.GetFinalNorm())) + print("MINRES solver took " + str(solve_time) + "s.") + # 12. Create the grid functions u and p. Compute the L2 error norms. + u = mfem.GridFunction() p = mfem.GridFunction() - u.MakeRef(R_space, x.GetBlock(0), 0) + if dg: + u_broken = mfem.GridFunction(R_space, x.GetBlock(0), 0) + coeff = mfem.VectorGridFunctionCoefficient(u_broken) + u.SetSpace(R_space_dg) + u.ProjectCoefficient(coeff) + else: + u.MakeRef(R_space, x.GetBlock(0), 0) p.MakeRef(W_space, x.GetBlock(1), 0) order_quad = max(2, 2*order+1) @@ -369,9 +381,14 @@ def EvalValue(self, x): sout_u = mfem.socketstream("localhost", 19916) sout_u.precision(8) sout_u << "solution\n" << mesh << u << "window_title 'Velocity'\n" + sout_u << "keys Rljvvvvvmmc" + sout_u.flush() + sout_p = mfem.socketstream("localhost", 19916) sout_p.precision(8) sout_p << "solution\n" << mesh << p << "window_title 'Pressure'\n" + sout_p << "keys Rljmmc" + sout_p.flush() if __name__ == "__main__": @@ -383,17 +400,17 @@ def EvalValue(self, x): action='store', type=str, help='Mesh file to use.') parser.add_argument("-nx", "--ncells-x", - action='store', type=int, default=0, + action='store', type=int, default=5, help="Number of cells in x.") parser.add_argument("-ny", "--ncells-y", - action='store', type=int, default=0, + action='store', type=int, default=5, help="Number of cells in y.") parser.add_argument('-o', '--order', action='store', default=1, type=int, - help="Finite element order (polynomial degree)."); + help="Finite element order (polynomial degree).") parser.add_argument("-dg", "--discontinuous", action='store_true', - help="Enable DG elements for fluxes."); + help="Enable DG elements for fluxes.") parser.add_argument("-brt", "--broken-RT", action='store_true', help="Enable broken RT elements for fluxes.") @@ -419,13 +436,15 @@ def EvalValue(self, x): action='store_true', help='Enable GLVis visualization') - args = parser.parse_args() parser.print_options(args) order = args.order - meshfile = expanduser( - join(dirname(__file__), '..', 'data', args.mesh)) + if args.mesh != "": + meshfile = expanduser( + join(dirname(__file__), '..', 'data', args.mesh)) + else: + meshfile = "" visualization = args.visualization device = args.device pa = args.partial_assembly @@ -433,12 +452,12 @@ def EvalValue(self, x): run(order=order, meshfile=meshfile, nx=args.ncells_x, - ny=args.ncells_y, - dg=args.discontinous, - brt=args.broken_RF, + ny=args.ncells_y, + dg=args.discontinuous, + brt=args.broken_RT, td=args.stab_diff, - hb=args.hybridizatoin, - rd=args.reduciton, + hb=args.hybridization, + rd=args.reduction, sn=args.solution_norm, visualization=visualization, device=device, diff --git a/mfem/ser.py b/mfem/ser.py index a66f05e0..5c8b6d9c 100644 --- a/mfem/ser.py +++ b/mfem/ser.py @@ -85,6 +85,7 @@ from mfem._ser.submesh import * from mfem._ser.transfermap import * from mfem._ser.hyperbolic import * +from mfem._ser.darcyform import * import mfem._ser.array as array import mfem._ser.blockoperator as blockoperator @@ -103,6 +104,7 @@ import mfem._ser.sparsemat as sparsemat import mfem._ser.tmop_modules as tmop + # # modules not a part of standard build #