11from sympy import Eq
22from devito .symbolics import retrieve_indexed , retrieve_dimensions
33from devito .petsc .types .equation import ConstrainBC
4- from devito .types .dimension import CustomBoundSubDimension
4+ from devito .types .dimension import CustomBoundSubDimension , SpaceDimension , CustomBoundSpaceDimension
55from devito import Min , Max
66
77
@@ -26,13 +26,17 @@ def constrain_essential_bcs(expressions, **kwargs):
2626 new_exprs .append (e )
2727 continue
2828
29+
30+ # this needs to be applied to all space dimensions and all subdims that are "MIDDLE"
31+
2932 indexeds = retrieve_indexed (e )
3033 dims = retrieve_dimensions ([i for j in indexeds for i in j .indices ], mode = 'unique' )
3134 # implicit_dims = set(e.implicit_dims)
3235 dims .update (e .implicit_dims )
3336 # from IPython import embed; embed()
3437 dims = [d for d in dims if d .is_Sub and not d .local ]
3538
39+
3640 if not dims :
3741 new_exprs .append (e )
3842 continue
@@ -80,4 +84,137 @@ def constrain_essential_bcs(expressions, **kwargs):
8084 new_e = new_e ._rebuild (implicit_dims = implicit_dims_new )
8185 new_exprs .append (new_e )
8286
83- return new_exprs
87+ return new_exprs
88+
89+
90+
91+ # def constrain_essential_bcs(expressions, **kwargs):
92+ # """TODO: improve docs ..Modify the subdims used in ConstrainEssentialBC equations ... to locally
93+ # constrain nodes (including non owned halo nodes) ....."""
94+
95+ # mapper = {}
96+ # new_exprs = []
97+
98+ # # build mapper
99+ # for e in expressions:
100+ # if not isinstance(e, ConstrainBC):
101+ # # new_exprs.append(e)
102+ # continue
103+ # # this needs to be applied to all space dimensions and all subdims that are "MIDDLE"
104+
105+ # indexeds = retrieve_indexed(e)
106+ # dims = retrieve_dimensions([i for j in indexeds for i in j.indices], mode='unique')
107+ # # from IPython import embed; embed()
108+ # # implicit_dims = set(e.implicit_dims)
109+ # dims.update(e.implicit_dims)
110+ # # from IPython import embed; embed()
111+
112+ # # Collect non-local subdims (which would be generated from a "middle" subdomain)
113+ # # "local" subdims are generated from 'left' and 'right" subdomains and do not need to be considered
114+ # # because by construction, left and right subdimensions cannot cross ranks
115+ # subdims = [d for d in dims if d.is_Sub and not d.local]
116+
117+ # space_dims = [d for d in dims if isinstance(d, SpaceDimension)]
118+
119+ # relevant_dims = subdims + space_dims
120+ # # relevant_dims = subdims
121+
122+
123+ # if not relevant_dims:
124+ # continue
125+
126+
127+ # # first map the subdims
128+ # for d in subdims:
129+ # # replace the dim with a new one that has a different symbolic_min and symbolic_max
130+
131+ # # obvs shouldn't be obtained from indexeds[0]
132+ # # USE e.lhs function -> the one that the BC is being applied to
133+ # # from IPython import embed; embed()
134+ # # f._size_nodomain.left
135+ # # from IPython import embed; embed()
136+ # # halo_size_left = indexeds[0].function._size_halo[d].left
137+ # # halo_size_right = indexeds[0].function._size_halo[d].right
138+
139+ # halo_size_left = 2
140+ # halo_size_right = 2
141+
142+ # from devito.petsc.types.dimension import SubDimMax, SubDimMin
143+
144+ # # TODO: change name..
145+
146+ # # in theory this class shoulod just take in d
147+ # # TODO: use unique name
148+ # sregistry = kwargs.get('sregistry')
149+ # subdim_max = SubDimMax(sregistry.make_name(prefix=d.name + '_max'), subdim=d)
150+ # subdim_min = SubDimMin(sregistry.make_name(prefix=d.name + '_min'), subdim=d)
151+ # # from IPython import embed; embed()
152+ # # unique_name
153+ # new_dim = CustomBoundSubDimension(
154+ # name=d.name,
155+ # parent=d.parent,
156+ # thickness=d.thickness,
157+ # local=d.local,
158+ # custom_left=Max(subdim_min, d.parent.symbolic_min - halo_size_left),
159+ # custom_right=Min(subdim_max, d.parent.symbolic_max + halo_size_right)
160+ # )
161+ # mapper[d] = new_dim
162+
163+ # # then tackle space dims
164+ # for d in space_dims:
165+ # mapper[d] = d
166+ # # for d in space_dims:
167+ # # halo_size_left = 2
168+ # # halo_size_right = 2
169+
170+ # # from devito.petsc.types.dimension import SpaceDimMax, SpaceDimMin
171+
172+ # # # TODO: change name..
173+
174+ # # # in theory this class shoulod just take in d
175+ # # # TODO: use unique name
176+ # # sregistry = kwargs.get('sregistry')
177+ # # space_dim_max = SpaceDimMax(sregistry.make_name(prefix=d.name + '_max'), space_dim=d)
178+ # # space_dim_min = SpaceDimMin(sregistry.make_name(prefix=d.name + '_min'), space_dim=d)
179+
180+ # # # unique_name
181+ # # new_dim = CustomBoundSpaceDimension(
182+ # # name=d.name,
183+ # # custom_left=Max(space_dim_min, d.symbolic_min - halo_size_left),
184+ # # custom_right=Min(space_dim_max, d.symbolic_max + halo_size_right)
185+ # # )
186+ # # mapper[d] = new_dim
187+ # # # from IPython import embed; embed()
188+
189+
190+ # # from IPython import embed; embed()
191+
192+
193+ # for e in expressions:
194+
195+ # if not isinstance(e, ConstrainBC):
196+ # new_exprs.append(e)
197+ # continue
198+
199+ # indexeds = retrieve_indexed(e)
200+ # dims = retrieve_dimensions([i for j in indexeds for i in j.indices], mode='unique')
201+ # # implicit_dims = set(e.implicit_dims)
202+ # dims.update(e.implicit_dims)
203+ # # from IPython import embed; embed()
204+ # subdims = [d for d in dims if d.is_Sub and not d.local]
205+
206+ # space_dims = [d for d in dims if isinstance(d, SpaceDimension)]
207+
208+ # relevant_dims = subdims + space_dims
209+
210+ # if not relevant_dims:
211+ # new_exprs.append(e)
212+ # continue
213+
214+ # new_e = e.subs(mapper)
215+ # if e.implicit_dims:
216+ # implicit_dims_new = tuple(mapper.get(d, d) for d in e.implicit_dims)
217+ # new_e = new_e._rebuild(implicit_dims=implicit_dims_new)
218+ # new_exprs.append(new_e)
219+
220+ # return new_exprs
0 commit comments