From e218a41e87dc0b9853d925b70525b11010de19bd Mon Sep 17 00:00:00 2001 From: John Stephens Date: Tue, 15 Jul 2025 18:58:24 -0500 Subject: [PATCH 1/6] pipe bounds into Galerkin-in-time --- irksome/base_time_stepper.py | 2 +- irksome/discontinuous_galerkin_stepper.py | 8 ++++++-- irksome/galerkin_stepper.py | 8 ++++++-- 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/irksome/base_time_stepper.py b/irksome/base_time_stepper.py index cf0410c4..4ea4383d 100644 --- a/irksome/base_time_stepper.py +++ b/irksome/base_time_stepper.py @@ -167,7 +167,7 @@ def get_stage_bounds(self, bounds=None): if upper is None: sub = Function(Vbig).assign(PETSc.INFINITY) - if bounds_type == "stage": + if (bounds_type == "stage") or (bounds_type == "galerkin"): if lower is not None: dats = [lower.dat] * (self.num_stages) slb = Function(Vbig, val=flatten_dats(dats)) diff --git a/irksome/discontinuous_galerkin_stepper.py b/irksome/discontinuous_galerkin_stepper.py index d083cc42..68592542 100644 --- a/irksome/discontinuous_galerkin_stepper.py +++ b/irksome/discontinuous_galerkin_stepper.py @@ -177,9 +177,10 @@ class DiscontinuousGalerkinTimeStepper(StageCoupledTimeStepper): instance to be passed to a `firedrake.MixedVectorSpaceBasis` over the larger space associated with the Runge-Kutta method + :arg bounds: An optional kwarg used in certain bounds-constrained methods. """ def __init__(self, F, order, t, dt, u0, bcs=None, basis_type=None, - quadrature=None, **kwargs): + quadrature=None, bounds=None, **kwargs): assert order >= 0 self.order = order self.basis_type = basis_type @@ -198,7 +199,10 @@ def __init__(self, F, order, t, dt, u0, bcs=None, basis_type=None, self.update_b = vecconst(self.el.tabulate(0, (1.0,))[(0,)]) - super().__init__(F, t, dt, u0, num_stages, bcs=bcs, **kwargs) + if bounds is not None: + assert (basis_type is None) or (basis_type == "Lagrange") or (basis_type == "Bernstein") + + super().__init__(F, t, dt, u0, num_stages, bcs=bcs, bounds=bounds, **kwargs) def get_form_and_bcs(self, stages, basis_type=None, order=None, quadrature=None, F=None): if basis_type is None: diff --git a/irksome/galerkin_stepper.py b/irksome/galerkin_stepper.py index 6e919b66..d6830416 100644 --- a/irksome/galerkin_stepper.py +++ b/irksome/galerkin_stepper.py @@ -169,9 +169,10 @@ class GalerkinTimeStepper(StageCoupledTimeStepper): instance to be passed to a `firedrake.MixedVectorSpaceBasis` over the larger space associated with the Runge-Kutta method + :arg bounds: An optional kwarg used in certain bounds-constrained methods. """ def __init__(self, F, order, t, dt, u0, bcs=None, basis_type=None, - quadrature=None, **kwargs): + quadrature=None, bounds=None, **kwargs): assert order >= 1 self.order = order self.basis_type = basis_type @@ -186,7 +187,10 @@ def __init__(self, F, order, t, dt, u0, bcs=None, basis_type=None, self.quadrature = quadrature assert np.size(quadrature.get_points()) >= order - super().__init__(F, t, dt, u0, order, bcs=bcs, **kwargs) + if bounds is not None: + assert (basis_type is None) or (basis_type == "Lagrange") or (basis_type == "Bernstein") + + super().__init__(F, t, dt, u0, order, bcs=bcs, bounds=bounds, **kwargs) def get_form_and_bcs(self, stages, basis_type=None, order=None, quadrature=None, F=None): if basis_type is None: From 713debd4088d38fe975ee8dd7f1bb77801b7f5fd Mon Sep 17 00:00:00 2001 From: John Stephens Date: Tue, 15 Jul 2025 21:19:31 -0500 Subject: [PATCH 2/6] some (failing) tests --- tests/test_galerkinBounds.py | 328 +++++++++++++++++++++++++++++++++++ 1 file changed, 328 insertions(+) create mode 100644 tests/test_galerkinBounds.py diff --git a/tests/test_galerkinBounds.py b/tests/test_galerkinBounds.py new file mode 100644 index 00000000..ef3d1088 --- /dev/null +++ b/tests/test_galerkinBounds.py @@ -0,0 +1,328 @@ +import numpy as np +import pytest +from firedrake import (DirichletBC, Function, FunctionSpace, SpatialCoordinate, TestFunction, + TestFunctions, UnitSquareMesh, diff, div, dx, exp, grad, inner, + pi, sin, split, tanh, sqrt, NonlinearVariationalProblem, + NonlinearVariationalSolver) +from irksome import (Dt, MeshConstant, BoundsConstrainedDirichletBC, + GalerkinTimeStepper, DiscontinuousGalerkinTimeStepper) +from ufl.algorithms import expand_derivatives + +from FIAT.quadrature import (RadauQuadratureLineRule, GaussLobattoLegendreQuadratureLineRule) +from FIAT import ufc_simplex + +lu_params = { + "snes_type": "ksponly", + "ksp_type": "preonly", + "mat_type": "aij", + "pc_type": "lu" +} + +vi_params = { + "snes_type": "vinewtonrsls", + "snes_max_it": 300, + "snes_atol": 1.e-8, + "ksp_type": "preonly", + "pc_type": "lu" +} + + +def heat_CG(quad_rule, order, basis_type, bounds_type): + + N = 16 + msh = UnitSquareMesh(N, N) + V = FunctionSpace(msh, "Lagrange", 1) + + MC = MeshConstant(msh) + dt = MC.Constant(2 / N) + t = MC.Constant(0.0) + + x, y = SpatialCoordinate(msh) + + uexact = 0.5 * exp(-t) * (1 + (tanh((0.1 - sqrt((x - 0.5) ** 2 + (y - 0.5) ** 2)) / 0.015))) + + rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact)) + + v = TestFunction(V) + u_init = Function(V) + + G = inner(u_init - uexact, v) * dx + + nlvp = NonlinearVariationalProblem(G, u_init) + nlvs = NonlinearVariationalSolver(nlvp, solver_parameters=vi_params) + + lb = Function(V) + ub = Function(V) + + ub.assign(np.inf) + lb.assign(0.0) + + nlvs.solve(bounds=(lb, ub)) + u_c = Function(V) + u_c.assign(u_init) + + v_c = TestFunction(V) + + F_c = (inner(Dt(u_c), v_c) * dx + inner(grad(u_c), grad(v_c)) * dx - inner(rhs, v_c) * dx) + + bounds = (bounds_type, lb, ub) + + bc = BoundsConstrainedDirichletBC(V, uexact, "on_boundary", (lb, ub), solver_parameters=vi_params) + + if quad_rule is not None: + quad_rule = quad_rule(ufc_simplex(1), order+1) + + stepper = GalerkinTimeStepper(F_c, order, t, dt, u_c, quadrature=quad_rule, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) + + violations_for_constrained_method = [] + + for _ in range(5): + stepper.advance() + t += dt + min_value_c = min(u_c.dat.data) + print(min_value_c) + if min_value_c < 0: + violations_for_constrained_method.append(min_value_c) + + return violations_for_constrained_method + + +def heat_DG(quad_rule, order, basis_type, bounds_type): + + N = 16 + msh = UnitSquareMesh(N, N) + V = FunctionSpace(msh, "Lagrange", 1) + + MC = MeshConstant(msh) + dt = MC.Constant(2 / N) + t = MC.Constant(0.0) + + x, y = SpatialCoordinate(msh) + + uexact = 0.5 * exp(-t) * (1 + (tanh((0.1 - sqrt((x - 0.5) ** 2 + (y - 0.5) ** 2)) / 0.015))) + + rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact)) + + v = TestFunction(V) + u_init = Function(V) + + G = inner(u_init - uexact, v) * dx + + nlvp = NonlinearVariationalProblem(G, u_init) + nlvs = NonlinearVariationalSolver(nlvp, solver_parameters=vi_params) + + lb = Function(V) + ub = Function(V) + + ub.assign(np.inf) + lb.assign(0.0) + + nlvs.solve(bounds=(lb, ub)) + u_c = Function(V) + u_c.assign(u_init) + + v_c = TestFunction(V) + + F_c = (inner(Dt(u_c), v_c) * dx + inner(grad(u_c), grad(v_c)) * dx - inner(rhs, v_c) * dx) + + bounds = (bounds_type, lb, ub) + + bc = BoundsConstrainedDirichletBC(V, uexact, "on_boundary", (lb, ub), solver_parameters=vi_params) + + if quad_rule is not None: + quad_rule = quad_rule(ufc_simplex(1), order+1) + + stepper = DiscontinuousGalerkinTimeStepper(F_c, order, t, dt, u_c, quadrature=quad_rule, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) + + violations_for_constrained_method = [] + + for _ in range(5): + stepper.advance() + t += dt + min_value_c = min(u_c.dat.data) + print(min_value_c) + if min_value_c < 0: + violations_for_constrained_method.append(min_value_c) + + return violations_for_constrained_method + + +def wave_H1_CG(quad_rule, order, basis_type, bounds_type): + + msh = UnitSquareMesh(16, 16, name='mesh') + W = FunctionSpace(msh, "Lagrange", 1) + Z = W * W + + x, y = SpatialCoordinate(msh) + + MC = MeshConstant(msh) + t = MC.Constant(0.0) + dt = MC.Constant(np.sqrt(2) / 160.0) + Tf = MC.Constant(np.sqrt(2) / 10.0) + + u_init = sin(5*pi*x)*sin(5*pi*y) + + uv = Function(Z, name='uv') + + phi0 = TestFunction(W) + + bc = DirichletBC(W, 0, "on_boundary") + + projection_problem = NonlinearVariationalProblem( + inner(uv.subfunctions[0] - u_init, phi0) * dx, + uv.subfunctions[0], bcs=bc) + + projection_solver = NonlinearVariationalSolver( + projection_problem, solver_parameters=vi_params) + + upper = Function(Z) + lower = Function(Z) + + upper.subfunctions[1].assign(np.inf) + lower.subfunctions[1].assign(-np.inf) + + upper.subfunctions[0].assign(1.0) + lower.subfunctions[0].assign(-1.0) + + projection_solver.solve(bounds=(lower.subfunctions[0], upper.subfunctions[0])) + u, v = split(uv) + + phi, psi = TestFunctions(Z) + F = (inner(Dt(u), phi) * dx - inner(v, phi) * dx + + inner(Dt(v), psi) * dx + inner(grad(u), grad(psi)) * dx) + + bc = [DirichletBC(Z.sub(0), 0, "on_boundary"), DirichletBC(Z.sub(1), 0, "on_boundary")] + + bounds = (bounds_type, lower, upper) + + if quad_rule is not None: + quad_rule = quad_rule(ufc_simplex(1), order+1) + + stepper = GalerkinTimeStepper(F, order, t, dt, uv, quadrature=quad_rule, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) + + bounds_violations = [] + + while (float(t) < float(Tf)): + if float(t) + float(dt) > float(Tf): + dt.assign(float(Tf) - float(t)) + + stepper.advance() + t += dt + + min_val = min(uv.subfunctions[0].dat.data) + max_val = max(uv.subfunctions[0].dat.data) + + if min_val < -1.0: + bounds_violations.append(min_val) + if max_val > 1.0: + bounds_violations.append(max_val) + + return bounds_violations + + +def wave_H1_DG(quad_rule, order, basis_type, bounds_type): + msh = UnitSquareMesh(16, 16, name='mesh') + W = FunctionSpace(msh, "Lagrange", 1) + Z = W * W + + x, y = SpatialCoordinate(msh) + + MC = MeshConstant(msh) + t = MC.Constant(0.0) + dt = MC.Constant(np.sqrt(2) / 160.0) + Tf = MC.Constant(np.sqrt(2) / 10.0) + + u_init = sin(5*pi*x)*sin(5*pi*y) + + uv = Function(Z, name='uv') + + phi0 = TestFunction(W) + + bc = DirichletBC(W, 0, "on_boundary") + + projection_problem = NonlinearVariationalProblem( + inner(uv.subfunctions[0] - u_init, phi0) * dx, + uv.subfunctions[0], bcs=bc) + + projection_solver = NonlinearVariationalSolver( + projection_problem, solver_parameters=vi_params) + + upper = Function(Z) + lower = Function(Z) + + upper.subfunctions[1].assign(np.inf) + lower.subfunctions[1].assign(-np.inf) + + upper.subfunctions[0].assign(1.0) + lower.subfunctions[0].assign(-1.0) + + projection_solver.solve(bounds=(lower.subfunctions[0], upper.subfunctions[0])) + u, v = split(uv) + + phi, psi = TestFunctions(Z) + F = (inner(Dt(u), phi) * dx - inner(v, phi) * dx + + inner(Dt(v), psi) * dx + inner(grad(u), grad(psi)) * dx) + + bc = [DirichletBC(Z.sub(0), 0, "on_boundary"), DirichletBC(Z.sub(1), 0, "on_boundary")] + + bounds = (bounds_type, lower, upper) + + if quad_rule is not None: + quad_rule = quad_rule(ufc_simplex(1), order+1) + + stepper = DiscontinuousGalerkinTimeStepper(F, order, t, dt, uv, quadrature=quad_rule, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) + + bounds_violations = [] + + while (float(t) < float(Tf)): + if float(t) + float(dt) > float(Tf): + dt.assign(float(Tf) - float(t)) + + stepper.advance() + t += dt + + min_val = min(uv.subfunctions[0].dat.data) + max_val = max(uv.subfunctions[0].dat.data) + + if min_val < -1.0: + bounds_violations.append(min_val) + if max_val > 1.0: + bounds_violations.append(max_val) + + return bounds_violations + + +@pytest.mark.parametrize('order', (1, 2, 3)) +@pytest.mark.parametrize('quad_rule', [None, GaussLobattoLegendreQuadratureLineRule]) +@pytest.mark.parametrize('basis_type', ('Bernstein', 'Lagrange')) +@pytest.mark.parametrize('bounds_type', ("galerkin", )) +def test_heat_CG_bounds(quad_rule, order, basis_type, bounds_type): + error_list = heat_CG(quad_rule, order, basis_type, bounds_type) + assert len(error_list) == 0 + + +@pytest.mark.parametrize('order', (1, 2, 3)) +@pytest.mark.parametrize('quad_rule', [None, RadauQuadratureLineRule]) +@pytest.mark.parametrize('basis_type', ('Bernstein', 'Lagrange')) +@pytest.mark.parametrize('bounds_type', ("galerkin", )) +def test_heat_DG_bounds(quad_rule, order, basis_type, bounds_type): + error_list = heat_DG(quad_rule, order, basis_type, bounds_type) + assert len(error_list) == 0 + + +@pytest.mark.parametrize('order', (1, 2, 3)) +@pytest.mark.parametrize('quad_rule', [None, GaussLobattoLegendreQuadratureLineRule]) +@pytest.mark.parametrize('basis_type', ('Bernstein', 'Lagrange')) +@pytest.mark.parametrize('bounds_type', ("galerkin", )) +def test_wave_CG_bounds(quad_rule, order, basis_type, bounds_type): + error_list = wave_H1_CG(quad_rule, order, basis_type, bounds_type) + assert len(error_list) == 0 + + +@pytest.mark.parametrize('order', (1, 2, 3)) +@pytest.mark.parametrize('quad_rule', [None, RadauQuadratureLineRule]) +@pytest.mark.parametrize('basis_type', ('Bernstein', 'Lagrange')) +@pytest.mark.parametrize('bounds_type', ("galerkin", )) +def test_wave_DG_bounds(quad_rule, order, basis_type, bounds_type): + error_list = wave_H1_DG(quad_rule, order, basis_type, bounds_type) + assert len(error_list) == 0 From ae9ff93f1b6227f110116fba28571dd0b51ffab0 Mon Sep 17 00:00:00 2001 From: Robert Kirby Date: Thu, 17 Jul 2025 11:43:27 -0400 Subject: [PATCH 3/6] Update irksome/discontinuous_galerkin_stepper.py Co-authored-by: Pablo Brubeck --- irksome/discontinuous_galerkin_stepper.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/irksome/discontinuous_galerkin_stepper.py b/irksome/discontinuous_galerkin_stepper.py index 68592542..1dc42d9c 100644 --- a/irksome/discontinuous_galerkin_stepper.py +++ b/irksome/discontinuous_galerkin_stepper.py @@ -200,7 +200,7 @@ def __init__(self, F, order, t, dt, u0, bcs=None, basis_type=None, self.update_b = vecconst(self.el.tabulate(0, (1.0,))[(0,)]) if bounds is not None: - assert (basis_type is None) or (basis_type == "Lagrange") or (basis_type == "Bernstein") + assert (basis_type is None) or (basis_type != "integral") super().__init__(F, t, dt, u0, num_stages, bcs=bcs, bounds=bounds, **kwargs) From cabe9ec2993927b38d04612b2a1c9b76093b1c1c Mon Sep 17 00:00:00 2001 From: John Stephens Date: Thu, 17 Jul 2025 13:27:50 -0500 Subject: [PATCH 4/6] updating tests --- tests/test_galerkinBounds.py | 110 ++++++++++++++++++++++++++++------- 1 file changed, 90 insertions(+), 20 deletions(-) diff --git a/tests/test_galerkinBounds.py b/tests/test_galerkinBounds.py index ef3d1088..5164a5b7 100644 --- a/tests/test_galerkinBounds.py +++ b/tests/test_galerkinBounds.py @@ -9,7 +9,7 @@ from ufl.algorithms import expand_derivatives from FIAT.quadrature import (RadauQuadratureLineRule, GaussLobattoLegendreQuadratureLineRule) -from FIAT import ufc_simplex +from FIAT import (ufc_simplex, Lagrange) lu_params = { "snes_type": "ksponly", @@ -20,7 +20,7 @@ vi_params = { "snes_type": "vinewtonrsls", - "snes_max_it": 300, + "snes_max_it": 150, "snes_atol": 1.e-8, "ksp_type": "preonly", "pc_type": "lu" @@ -70,9 +70,11 @@ def heat_CG(quad_rule, order, basis_type, bounds_type): bc = BoundsConstrainedDirichletBC(V, uexact, "on_boundary", (lb, ub), solver_parameters=vi_params) if quad_rule is not None: - quad_rule = quad_rule(ufc_simplex(1), order+1) + quad = quad_rule(ufc_simplex(1), order+1) + else: + quad = None - stepper = GalerkinTimeStepper(F_c, order, t, dt, u_c, quadrature=quad_rule, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) + stepper = GalerkinTimeStepper(F_c, order, t, dt, u_c, quadrature=quad, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) violations_for_constrained_method = [] @@ -130,9 +132,73 @@ def heat_DG(quad_rule, order, basis_type, bounds_type): bc = BoundsConstrainedDirichletBC(V, uexact, "on_boundary", (lb, ub), solver_parameters=vi_params) if quad_rule is not None: - quad_rule = quad_rule(ufc_simplex(1), order+1) + quad = quad_rule(ufc_simplex(1), order+1) + else: + quad = None - stepper = DiscontinuousGalerkinTimeStepper(F_c, order, t, dt, u_c, quadrature=quad_rule, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) + stepper = DiscontinuousGalerkinTimeStepper(F_c, order, t, dt, u_c, quadrature=quad, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) + + violations_for_constrained_method = [] + + for _ in range(5): + stepper.advance() + t += dt + min_value_c = min(u_c.dat.data) + print(min_value_c) + if min_value_c < 0: + violations_for_constrained_method.append(min_value_c) + + return violations_for_constrained_method + + +def heat_DG(quad_rule, order, basis_type, bounds_type): + + N = 16 + msh = UnitSquareMesh(N, N) + V = FunctionSpace(msh, "Lagrange", 1) + + MC = MeshConstant(msh) + dt = MC.Constant(2 / N) + t = MC.Constant(0.0) + + x, y = SpatialCoordinate(msh) + + uexact = 0.5 * exp(-t) * (1 + (tanh((0.1 - sqrt((x - 0.5) ** 2 + (y - 0.5) ** 2)) / 0.015))) + + rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact)) + + v = TestFunction(V) + u_init = Function(V) + + G = inner(u_init - uexact, v) * dx + + nlvp = NonlinearVariationalProblem(G, u_init) + nlvs = NonlinearVariationalSolver(nlvp, solver_parameters=vi_params) + + lb = Function(V) + ub = Function(V) + + ub.assign(np.inf) + lb.assign(0.0) + + nlvs.solve(bounds=(lb, ub)) + u_c = Function(V) + u_c.assign(u_init) + + v_c = TestFunction(V) + + F_c = (inner(Dt(u_c), v_c) * dx + inner(grad(u_c), grad(v_c)) * dx - inner(rhs, v_c) * dx) + + bounds = (bounds_type, lb, ub) + + bc = BoundsConstrainedDirichletBC(V, uexact, "on_boundary", (lb, ub), solver_parameters=vi_params) + + if quad_rule is not None: + quad = quad_rule(ufc_simplex(1), order+1) + else: + quad = None + + stepper = DiscontinuousGalerkinTimeStepper(F_c, order, t, dt, u_c, quadrature=quad, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) violations_for_constrained_method = [] @@ -196,9 +262,11 @@ def wave_H1_CG(quad_rule, order, basis_type, bounds_type): bounds = (bounds_type, lower, upper) if quad_rule is not None: - quad_rule = quad_rule(ufc_simplex(1), order+1) + quad = quad_rule(ufc_simplex(1), order+1) + else: + quad = None - stepper = GalerkinTimeStepper(F, order, t, dt, uv, quadrature=quad_rule, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) + stepper = GalerkinTimeStepper(F, order, t, dt, uv, quadrature=quad, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) bounds_violations = [] @@ -268,9 +336,11 @@ def wave_H1_DG(quad_rule, order, basis_type, bounds_type): bounds = (bounds_type, lower, upper) if quad_rule is not None: - quad_rule = quad_rule(ufc_simplex(1), order+1) + quad = quad_rule(ufc_simplex(1), order+1) + else: + quad = None - stepper = DiscontinuousGalerkinTimeStepper(F, order, t, dt, uv, quadrature=quad_rule, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) + stepper = DiscontinuousGalerkinTimeStepper(F, order, t, dt, uv, quadrature=quad, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) bounds_violations = [] @@ -292,37 +362,37 @@ def wave_H1_DG(quad_rule, order, basis_type, bounds_type): return bounds_violations -@pytest.mark.parametrize('order', (1, 2, 3)) +@pytest.mark.parametrize('order', (1, 2)) @pytest.mark.parametrize('quad_rule', [None, GaussLobattoLegendreQuadratureLineRule]) -@pytest.mark.parametrize('basis_type', ('Bernstein', 'Lagrange')) -@pytest.mark.parametrize('bounds_type', ("galerkin", )) +@pytest.mark.parametrize('basis_type', ('Bernstein', 'gll')) +@pytest.mark.parametrize('bounds_type', ('stage', )) def test_heat_CG_bounds(quad_rule, order, basis_type, bounds_type): error_list = heat_CG(quad_rule, order, basis_type, bounds_type) assert len(error_list) == 0 -@pytest.mark.parametrize('order', (1, 2, 3)) +@pytest.mark.parametrize('order', (1, 2)) @pytest.mark.parametrize('quad_rule', [None, RadauQuadratureLineRule]) @pytest.mark.parametrize('basis_type', ('Bernstein', 'Lagrange')) -@pytest.mark.parametrize('bounds_type', ("galerkin", )) +@pytest.mark.parametrize('bounds_type', ('stage', )) def test_heat_DG_bounds(quad_rule, order, basis_type, bounds_type): error_list = heat_DG(quad_rule, order, basis_type, bounds_type) assert len(error_list) == 0 -@pytest.mark.parametrize('order', (1, 2, 3)) +@pytest.mark.parametrize('order', (1, 2)) @pytest.mark.parametrize('quad_rule', [None, GaussLobattoLegendreQuadratureLineRule]) -@pytest.mark.parametrize('basis_type', ('Bernstein', 'Lagrange')) -@pytest.mark.parametrize('bounds_type', ("galerkin", )) +@pytest.mark.parametrize('basis_type', ('Bernstein', 'gll')) +@pytest.mark.parametrize('bounds_type', ('stage', )) def test_wave_CG_bounds(quad_rule, order, basis_type, bounds_type): error_list = wave_H1_CG(quad_rule, order, basis_type, bounds_type) assert len(error_list) == 0 -@pytest.mark.parametrize('order', (1, 2, 3)) +@pytest.mark.parametrize('order', (1, 2)) @pytest.mark.parametrize('quad_rule', [None, RadauQuadratureLineRule]) @pytest.mark.parametrize('basis_type', ('Bernstein', 'Lagrange')) -@pytest.mark.parametrize('bounds_type', ("galerkin", )) +@pytest.mark.parametrize('bounds_type', ('stage', )) def test_wave_DG_bounds(quad_rule, order, basis_type, bounds_type): error_list = wave_H1_DG(quad_rule, order, basis_type, bounds_type) assert len(error_list) == 0 From 6d481584f91099ca28968e7a653c672c628b9ab2 Mon Sep 17 00:00:00 2001 From: John Stephens Date: Thu, 17 Jul 2025 15:12:37 -0500 Subject: [PATCH 5/6] updated tests --- irksome/discontinuous_galerkin_stepper.py | 2 +- irksome/galerkin_stepper.py | 2 +- tests/test_galerkinBounds.py | 12 ++++++------ 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/irksome/discontinuous_galerkin_stepper.py b/irksome/discontinuous_galerkin_stepper.py index 1dc42d9c..06ed4340 100644 --- a/irksome/discontinuous_galerkin_stepper.py +++ b/irksome/discontinuous_galerkin_stepper.py @@ -200,7 +200,7 @@ def __init__(self, F, order, t, dt, u0, bcs=None, basis_type=None, self.update_b = vecconst(self.el.tabulate(0, (1.0,))[(0,)]) if bounds is not None: - assert (basis_type is None) or (basis_type != "integral") + assert (basis_type is None) or (basis_type != "integral"), "The Bernstein basis or a pointwise basis is required for bounds constraints" super().__init__(F, t, dt, u0, num_stages, bcs=bcs, bounds=bounds, **kwargs) diff --git a/irksome/galerkin_stepper.py b/irksome/galerkin_stepper.py index d6830416..cc3ed9fd 100644 --- a/irksome/galerkin_stepper.py +++ b/irksome/galerkin_stepper.py @@ -188,7 +188,7 @@ def __init__(self, F, order, t, dt, u0, bcs=None, basis_type=None, assert np.size(quadrature.get_points()) >= order if bounds is not None: - assert (basis_type is None) or (basis_type == "Lagrange") or (basis_type == "Bernstein") + assert (basis_type is None) or (basis_type != "integral"), "The Bernstein basis or a pointwise basis is required for bounds constraints" super().__init__(F, t, dt, u0, order, bcs=bcs, bounds=bounds, **kwargs) diff --git a/tests/test_galerkinBounds.py b/tests/test_galerkinBounds.py index 5164a5b7..fe4c7a84 100644 --- a/tests/test_galerkinBounds.py +++ b/tests/test_galerkinBounds.py @@ -9,7 +9,7 @@ from ufl.algorithms import expand_derivatives from FIAT.quadrature import (RadauQuadratureLineRule, GaussLobattoLegendreQuadratureLineRule) -from FIAT import (ufc_simplex, Lagrange) +from FIAT import ufc_simplex lu_params = { "snes_type": "ksponly", @@ -362,7 +362,7 @@ def wave_H1_DG(quad_rule, order, basis_type, bounds_type): return bounds_violations -@pytest.mark.parametrize('order', (1, 2)) +@pytest.mark.parametrize('order', (1, 2, 3)) @pytest.mark.parametrize('quad_rule', [None, GaussLobattoLegendreQuadratureLineRule]) @pytest.mark.parametrize('basis_type', ('Bernstein', 'gll')) @pytest.mark.parametrize('bounds_type', ('stage', )) @@ -371,7 +371,7 @@ def test_heat_CG_bounds(quad_rule, order, basis_type, bounds_type): assert len(error_list) == 0 -@pytest.mark.parametrize('order', (1, 2)) +@pytest.mark.parametrize('order', (1, 2, 3)) @pytest.mark.parametrize('quad_rule', [None, RadauQuadratureLineRule]) @pytest.mark.parametrize('basis_type', ('Bernstein', 'Lagrange')) @pytest.mark.parametrize('bounds_type', ('stage', )) @@ -380,7 +380,7 @@ def test_heat_DG_bounds(quad_rule, order, basis_type, bounds_type): assert len(error_list) == 0 -@pytest.mark.parametrize('order', (1, 2)) +@pytest.mark.parametrize('order', (1, 2, 3)) @pytest.mark.parametrize('quad_rule', [None, GaussLobattoLegendreQuadratureLineRule]) @pytest.mark.parametrize('basis_type', ('Bernstein', 'gll')) @pytest.mark.parametrize('bounds_type', ('stage', )) @@ -389,9 +389,9 @@ def test_wave_CG_bounds(quad_rule, order, basis_type, bounds_type): assert len(error_list) == 0 -@pytest.mark.parametrize('order', (1, 2)) +@pytest.mark.parametrize('order', (1, 2, 3)) @pytest.mark.parametrize('quad_rule', [None, RadauQuadratureLineRule]) -@pytest.mark.parametrize('basis_type', ('Bernstein', 'Lagrange')) +@pytest.mark.parametrize('basis_type', ('Bernstein', )) @pytest.mark.parametrize('bounds_type', ('stage', )) def test_wave_DG_bounds(quad_rule, order, basis_type, bounds_type): error_list = wave_H1_DG(quad_rule, order, basis_type, bounds_type) From d9b5e4d0e972ae699ce0309429dcc4cb965c17fb Mon Sep 17 00:00:00 2001 From: "Robert C. Kirby" Date: Thu, 17 Jul 2025 17:17:02 -0400 Subject: [PATCH 6/6] remove duplicate function --- tests/test_galerkinBounds.py | 62 ------------------------------------ 1 file changed, 62 deletions(-) diff --git a/tests/test_galerkinBounds.py b/tests/test_galerkinBounds.py index fe4c7a84..dfa7f283 100644 --- a/tests/test_galerkinBounds.py +++ b/tests/test_galerkinBounds.py @@ -151,68 +151,6 @@ def heat_DG(quad_rule, order, basis_type, bounds_type): return violations_for_constrained_method -def heat_DG(quad_rule, order, basis_type, bounds_type): - - N = 16 - msh = UnitSquareMesh(N, N) - V = FunctionSpace(msh, "Lagrange", 1) - - MC = MeshConstant(msh) - dt = MC.Constant(2 / N) - t = MC.Constant(0.0) - - x, y = SpatialCoordinate(msh) - - uexact = 0.5 * exp(-t) * (1 + (tanh((0.1 - sqrt((x - 0.5) ** 2 + (y - 0.5) ** 2)) / 0.015))) - - rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact)) - - v = TestFunction(V) - u_init = Function(V) - - G = inner(u_init - uexact, v) * dx - - nlvp = NonlinearVariationalProblem(G, u_init) - nlvs = NonlinearVariationalSolver(nlvp, solver_parameters=vi_params) - - lb = Function(V) - ub = Function(V) - - ub.assign(np.inf) - lb.assign(0.0) - - nlvs.solve(bounds=(lb, ub)) - u_c = Function(V) - u_c.assign(u_init) - - v_c = TestFunction(V) - - F_c = (inner(Dt(u_c), v_c) * dx + inner(grad(u_c), grad(v_c)) * dx - inner(rhs, v_c) * dx) - - bounds = (bounds_type, lb, ub) - - bc = BoundsConstrainedDirichletBC(V, uexact, "on_boundary", (lb, ub), solver_parameters=vi_params) - - if quad_rule is not None: - quad = quad_rule(ufc_simplex(1), order+1) - else: - quad = None - - stepper = DiscontinuousGalerkinTimeStepper(F_c, order, t, dt, u_c, quadrature=quad, basis_type=basis_type, bounds=bounds, bcs=bc, solver_parameters=vi_params) - - violations_for_constrained_method = [] - - for _ in range(5): - stepper.advance() - t += dt - min_value_c = min(u_c.dat.data) - print(min_value_c) - if min_value_c < 0: - violations_for_constrained_method.append(min_value_c) - - return violations_for_constrained_method - - def wave_H1_CG(quad_rule, order, basis_type, bounds_type): msh = UnitSquareMesh(16, 16, name='mesh')