Skip to content

Conversation

@ZoeLeibowitz
Copy link
Owner

@ZoeLeibowitz ZoeLeibowitz commented Feb 26, 2025

  • New machinery to solve for more than 1 field within a single matrix system.

TODOs:

  • Add more extensive tests
  • Utilise some of the new machinery to support Tensor Functions
  • Start wiki page/project board to track current implementation status

@JDBetteridge JDBetteridge self-requested a review March 5, 2025 17:38
Copy link
Collaborator

@JDBetteridge JDBetteridge left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There's a lot of code here most of which looks good, but I have two major concerns:

  1. I'm not sure that the (existing) tests include sufficient coverage for the new code being added.
  2. There isn't much in the way of documentation.
    Perhaps there needs to be a more general project page to cover integration?

Comment on lines 200 to 228
def _callback_builder(self):
sobjs = self.objbuilder.solver_objs
args = (self.injectsolve, self.objs, sobjs, self.timedep)
return (
CCBBuilder(*args, **self.kwargs)
if self.coupled else
CBBuilder(*args, **self.kwargs)
)

# Solver setup
self.solversetup = BaseSetup(
self.solver_objs, objs, injectsolve, self.cbbuilder
def _setup(self):
sobjs = self.objbuilder.solver_objs
args = (self.injectsolve, self.objs, sobjs, self.cbbuilder)
return (
CoupledSetup(*args, **self.kwargs)
if self.coupled else
BaseSetup(*args, **self.kwargs)
)

# Execute the solver
self.solve = Solver(
self.solver_objs, objs, injectsolve, iters,
self.cbbuilder, timedep=self.timedep
def _solver_execution(self):
sobjs = self.objbuilder.solver_objs
args = (
self.injectsolve, self.objs, sobjs,
self.iters, self.cbbuilder, self.timedep
)
return (
CoupledSolver(*args)
if self.coupled else
Solver(*args)
)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There's a lot of code repetition in these methods, can it be factored out?

@property
def _C_free_priority(self):
return 3
return 4
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, we need an enumeration or something, because what started off as a simple thing is now growing out of hand :)

@ZoeLeibowitz ZoeLeibowitz marked this pull request as ready for review March 6, 2025 18:41
if not issubclass(ctype, ctypes.Structure):
return None
except TypeError:
# E.g., `ctype` is of type `dtypes_lowering.CustomDtype`
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this comment to be moved inside the else?

# E.g., `ctype` is of type `dtypes_lowering.CustomDtype`
return None
if isinstance(obj, LocalCompositeObject):
ctype = obj
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

how can the ctype being a LocalCompositeObject? I struggle to follow here

ret.extend(self.pointer.free_symbols)
else:
assert issubclass(self.pointer._C_ctype, ctypes._Pointer)
assert isinstance(self.pointer, LocalCompositeObject) or \
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nitpicking, for code homogeneity : for multi-line statements, instead of the backslash, we use parenthesis for perfectly aligned indentation

assert (isinstance(....) or
             issubclass(...))

"""
# Create a symbol registry
kwargs.setdefault('sregistry', SymbolRegistry())
# Update based on changes in #2509
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We may need a more useful comment here

@property
def _C_free_priority(self):
return 3
return 4
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, we need an enumeration or something, because what started off as a simple thing is now growing out of hand :)

__str__ = __repr__


class CustomIntType(CustomDtype):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why a subclass and not just a specific instance?

raise ValueError("Symbolic RHS `%s` must be of type `int`, found "
"`%s` instead" % (rhs, rhs.dtype))

is_int_type = isinstance(rhs.dtype, type) and \
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe this should be moved inside a utility function?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To add: convention is to use parenthesis


is_commutative = True

# Set the operator priority higher than SymPy (10.0) to force the overridden
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can you elaborate what the issue is?

# TODO: Are users going to be able to use PETSc command line arguments?
# In firedrake, they have an options_prefix for each solver, enabling the use
# of command line options
Null = objs['Null']
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there any chance of this creating a KeyError? Ditto for two lines below

**self.kwargs
}

self.objbuilder = self._object_builder()
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could this just be a cached property rather than requiring a method? Ditto with some of the similar constructions further down in this init, which would reduce the number of lines of code.

self.args['timedep'] = self.timedep

self.cbbuilder = self._callback_builder()
self.args['cbbuilder'] = self.cbbuilder
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are any of these functions sensitive to the contents of self.args? I.e. do they need to be called in a specific order?

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes

Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why?

self.objbuilder = BaseObjectBuilder(injectsolve, **kwargs)
self.solver_objs = self.objbuilder.solver_objs
def _callback_builder(self):
return (CCBBuilder(**self.args) if self.coupled else CBBuilder(**self.args))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nitpick: unnecessary brackets? (Ditto for the next few methods)

def populate_matrix_context(efuncs, objs):
name = 'PopulateMatContext'

if name not in efuncs:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This hardcoded string check doesn't look great, is there no way to check for a class?

' + x_matvec_f[x + 2, y + 1]/lctx->h_y**2' + \
' - 2.0*x_matvec_f[x + 2, y + 2]/lctx->h_y**2' + \
' + x_matvec_f[x + 2, y + 3]/lctx->h_y**2'
'x_f[x + 1, y + 2]/ctx0->h_x**2' + \
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nitpick: You can use parenthesis to join strings for tidiness

# Non-coupled
petsc1 = PETScSolve(eq1, target=e)
petsc2 = PETScSolve(eq2, target=g)
with switchconfig(openmp=False, mpi=True):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Switchconfig for MPI shouldn't be needed due to the parallel tag on the test

enorm1 = norm(e)
gnorm1 = norm(g)

# reset
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nitpick: start comment with uppercase


@skipif('petsc')
def test_coupled_structs(self):
grid = Grid(shape=(11, 11))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This gets used in a lot of tests

grid = Grid(shape=(11, 11))
e = Function(name='e', grid=grid, space_order=2)
f = Function(name='f', grid=grid, space_order=2)
g = Function(name='g', grid=grid, space_order=2)
h = Function(name='h', grid=grid, space_order=2)

would a utility function to set these up be worthwhile?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or a list comprehension?

ccode, hcode = op.cinterface(force=True)

dirname = op._compiler.get_jit_dir()
assert os.path.isfile(os.path.join(dirname, "%s.c" % name))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fstring (see a couple of other places in this test)

@ZoeLeibowitz ZoeLeibowitz moved this to In Progress in PETSc Mar 12, 2025
@ZoeLeibowitz ZoeLeibowitz added the enhancement New feature or request label Mar 12, 2025
@ZoeLeibowitz ZoeLeibowitz merged commit c68dd4f into master Mar 19, 2025
29 checks passed
@github-project-automation github-project-automation bot moved this from In Progress to Done in PETSc Mar 19, 2025
ZoeLeibowitz added a commit that referenced this pull request Mar 21, 2025
* compiler/dsl: Add machinery to support coupled solvers
ZoeLeibowitz added a commit that referenced this pull request Mar 24, 2025
* compiler/dsl: Add machinery to support coupled solvers
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request

Projects

Status: Done

Development

Successfully merging this pull request may close these issues.

5 participants