Skip to content

Commit

Permalink
Removed if-else-statements for mesh type
Browse files Browse the repository at this point in the history
  • Loading branch information
lisawim committed Feb 23, 2024
1 parent e24ff2b commit 486cfbb
Show file tree
Hide file tree
Showing 2 changed files with 12 additions and 26 deletions.
31 changes: 10 additions & 21 deletions pySDC/projects/DAE/misc/ProblemDAE.py
Original file line number Diff line number Diff line change
Expand Up @@ -56,28 +56,17 @@ def solve_system(self, impl_sys, u0, t):
"""
me = self.dtype_u(self.init)

def implSysAsNumpy(unknowns, **kwargs):
me.diff[:] = unknowns[: np.size(me.diff)].reshape(me.diff.shape)
me.alg[:] = unknowns[np.size(me.diff) :].reshape(me.alg.shape)
def implSysFlatten(unknowns, **kwargs):
me[:] = unknowns.reshape(me.shape)

This comment has been minimized.

Copy link
@brownbaerchen

brownbaerchen Feb 23, 2024

Contributor

I am wondering if we can make this more efficient.
impl_sys wants a MultiComponentMesh, so we cannot just pass unknowns.reshape(self.init[0]). This would be nice because, if I understand correctly, it would not cause any copying of data, but create a view onto the same data with a different shape. Could we plug unknowns.reshape(self.init[0]).view(type(u0)) into impl_sys? Would this even avoid copying? What do you think @tlunet, @pancetta? If you don't think it's worth bothering, we can just go ahead and finally merge this!

This comment has been minimized.

Copy link
@lisawim

lisawim Feb 23, 2024

Author Collaborator

in implSysFlatten we can directly use unknowns, right? So, we'd have

def implSysFlatten(unknowns, **kwargs):
    sys = impl_sys(unknowns, **kwargs)
    return sys.flatten()

? This does also work. Ok, but unknowns still would be a numpy.1darray from the scipy routine..

This comment has been minimized.

Copy link
@lisawim

lisawim Feb 23, 2024

Author Collaborator

Forgot that.. I see your point!

This comment has been minimized.

Copy link
@lisawim

lisawim Feb 23, 2024

Author Collaborator

Hm, can we do something like unknowns.reshape(u0.copy().shape).view(type(u0))? Does that make sense? self.init[0] only stores nvars, but in case DAEMesh the shape that we need is (2,nvars).

I do not understand the point why we cannot use me.shape for that. Can you explain that to me?

This comment has been minimized.

Copy link
@brownbaerchen

brownbaerchen Feb 23, 2024

Contributor

You're right! me.shape is good! I wasn't reading the code throughly, my mistake!

This comment has been minimized.

Copy link
@lisawim

lisawim Feb 23, 2024

Author Collaborator

No worries! And anyway, you enforced me to do further improvements. Thanks a lot for that!

sys = impl_sys(me, **kwargs)
return np.append(sys.diff.flatten(), sys.alg.flatten()) # TODO: more efficient way?
return sys.flatten()

if type(me) == DAEMesh:
opt = root(
implSysAsNumpy,
np.append(u0.diff.flatten(), u0.alg.flatten()),
method='hybr',
tol=self.newton_tol,
)
me.diff[:] = opt.x[: np.size(me.diff)].reshape(me.diff.shape)
me.alg[:] = opt.x[np.size(me.diff) :].reshape(me.alg.shape)
else:
opt = root(
impl_sys,
u0,
method='hybr',
tol=self.newton_tol,
)
me[:] = opt.x
opt = root(
implSysFlatten,
u0.flatten(),
method='hybr',
tol=self.newton_tol,
)
me[:] = opt.x.reshape(me.shape)
self.work_counters['newton'].niter += opt.nfev
return me
7 changes: 2 additions & 5 deletions pySDC/projects/DAE/sweepers/fully_implicit_DAE.py
Original file line number Diff line number Diff line change
Expand Up @@ -107,11 +107,8 @@ def implSystem(params):
System to be solved as implicit function.
"""

if type(params) == DAEMesh:
params_mesh = P.dtype_f(params)
else:
params_mesh = P.dtype_f(P.init)
params_mesh[:] = params
params_mesh = P.dtype_f(P.init)
params_mesh[:] = params.reshape(params_mesh.shape)

# build parameters to pass to implicit function
local_u_approx = P.dtype_f(u_approx)
Expand Down

0 comments on commit 486cfbb

Please sign in to comment.