"Proper" way to run a paremeter sweep for steady state simulations?

[Kevin Modica]

Hi All,

I've been using dedalus to solve for steady state solutions to an advection-diffusion equation with different parameters. The naive way of doing it would be to just run the entire dedalus workflow in a for loop, but in the past I had issues with the distributor being created too many times and spyder crashing. I also suspect a lot of time is being wasted in remaking the fields and fourier system.

I was wondering if there were any recommended ways of setting up the problem to allow for this kind of high throughput analysis.

I've included a simple example script with the naive implementation, and would appreciate any advice on how an expert in dedalus would go about this.

The example runs quickly because this is a simple 1D problem, but for more involved problems, even a slight speedup in the loop can save me a bunch of time down the road.

Thanks,

Kevin

[Keaton Burns]

Hi Kevin,

As long as the parameters being adjusted are entering the equations via parameter fields, as is the case here, then actually all of the objects, up to and including the solver, can be re-used.  The only thing that needs to change across the loops is that the LBVP LHS matrices need to rebuilt, since the parameter is acting as an NCC. This can be triggered by passing the rebuild_matrices=True flag to solver.solve , which will recreate the LHS matrices using the new NCC values.

For complicated problems, building the matrices can take a lot of time so this may or may not speed things up overall, but it will reuse as much as possible (field objects, transform plans, etc.).

Best,

-Keaton

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/8a24c773-bafd-41db-9caa-df83dfe4faf7n%40googlegroups.com.

[Ben Brown]

Kevin,

      I was drafting up a script for you when I saw Keaton's response, and it's 100% right.

Here are two versions of the script, your original version ("20240103_SimpleDedalusSetup_original.py"), and a version that does the rebuild_matrices=True approach ("20240103_SimpleDedalusSetup.py").  I threw in some logging outputs so that I could double check that they both seem to be doing the same thing.

The latter is about 2x faster than the former.  Some of the 

This kind of approach gets particularly powerful, in my experience, once you start doing NLBVPs.  Then you can use the prior solution to your variables (e.g., u) as the initial guess for the new solution with nearby but changed parameters.  It's a simple form of continuation, and one that Jeff Oishi, I and some others here have enjoyed playing around with.  It can be pretty powerful for the right problems.

--Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/CADZXxBjQeyUfQxC0s4ncFZrqvm9%3Dj-NCPUsdirVyfpyZhZ6w0Q%40mail.gmail.com.

[Kevin Modica]

Thank you both for your advice!

Ben, I had a couple of minor question about the sample you included, I noticed you put on lines 25 - 30:

x = dist.local_grids(xbasis)
xg = dist.Field(name='x', bases=xbasis)
xg['g'] = x
V = dist.Field(bases=xbasis)
A = dist.Field()
V = -A*np.cos(xg)

Can I ask why you chose to create an xg field instead of using the local grid from the distributor? I'm guessing that the line V = dist.Field(bases=xbasis) is overwritten by the following two lines, and instead you need to use xg to define V as a composite field that has a basis in x.

Is it fair to say that

V = dist.Field(bases=xbasis)

A = dist.Field()

V['g'] = -A*np.cos(x)

is equivalent to

xg = dist.Field(name='x', bases=xbasis)
xg['g'] = x
A = dist.Field()
V = -A*np.cos(xg)

[Ben Brown]

Kevin,

     I defined xg as a field so that later in the analysis I could use xg directly in defining the diagnostic quantity.  

You're exactly right that the line

V = dist.Field(bases=xbasis)

is redundant and can (and should) be pruned out; using

V = -A*np.cos(xg)

with A and xg fields already does the job of creating the right object for V.

A needs to be a field, and defined before V, because we're going to be updating it in the loop.

Yes, the later two blocks that you quote both have the same broad outcome for V.  

I prefer the xg version because of the output task.  I think this version, with V=A*np.cos(xg) will also render nicely as A*cos(x) in the debugging outputs or when printing V.name, or at least would do so if we properly gave A a name when defining it (A = dist.Field(name='A')).

Ah, yep, if you do:

A = dist.Field(name='A')

xg = dist.Field(name='x', bases=xb) 
V = -A*np.cos(xg)

V.evaluate().name

you get:

'-1*A*cos(x)'

 --Ben

To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/46e419a4-c398-41e9-99e2-79996b4da764n%40googlegroups.com.