Trouble saving with evaluate_handlers inside time loop

[Sidney Holden]

Hey guys,

I'm having trouble with manually saving state after manipulating a field variable within a time loop. I'm trying to perform a time step, do some operation on the field, and then save the resulting state. 

Here's the idea:

...

problem = d3.IVP([u], namespace=locals())

...

solver = problem.build_solver(d3.RK222)

...

analysis = solver.evaluator.add_file_handler("data", max_writes=max_writes, mode="overwrite")

analysis.add_tasks(solver.state, layout='g')

while solver.proceed:

solver.step(time_step)

u['g'] = np.abs(u['g'])

if solver.iteration % save_cadence == 0:

solver.evaluator.evaluate_handlers([analysis], iteration=solver.iteration, sim_time=solver.sim_time) ...

But the evaluator still saves negative values of u.

Is there an obvious fix and/or better way of doing this kind of thing?

If not I'll make a minimal working example of exactly what I'm trying to do.

Cheers,
Sid

[Daniel Lecoanet]

Hi Sid,

It is not obvious to me what is going wrong. It would be useful to send a MWE.

Daniel

--
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 visit https://groups.google.com/d/msgid/dedalus-users/182a341d-1ac8-4b60-aa2f-357250e9742en%40googlegroups.com.

[Sidney Holden]

Thanks Daniel. I’ve attached a notebook with a MWE for debugging. It should run in ~10s. 

I think what’s going on is similar to what you mentioned here: https://groups.google.com/g/dedalus-users/c/LbZWI4N_t4Y/m/OSg-OmECCAAJ. Except I have a PSD matrix K I'm trying to solve for rather than a positive density.

I’ve played around with require_coeff_space and require_grid_space without any luck. The projection is doing the correct thing based off what gets printed directly, but the evaluator just doesn’t seem to save what I see on the grid. 

I’ve also tried solving for the components of a log-Cholesky factorization K = LL^T, but this blows up when the eigenvalues would otherwise drop below 0 (around T = 41.2). So the reason I’m trying to do this PSD projection each timestep is because the behavior also persists for high space and time resolutions, different timesteppers (RK443, SBDF2), etc

[Keaton Burns]

Hi Sid,

My guess is there’s a roundtrip transform before the output is saved, and this is affecting the values of K by removing the Nyquist component. The cone projection on the grid is likely introducing high frequencies since it has sharp thresholding.  Maybe replacing the maximum there with some version of softplus could help.

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 visit https://groups.google.com/d/msgid/dedalus-users/2184F56B-881C-442D-A706-9FA69C132F61%40flatironinstitute.org.

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/9364E311-8454-41A8-BAAF-978726294EE9%40northwestern.edu.

--
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 visit https://groups.google.com/d/msgid/dedalus-users/2184F56B-881C-442D-A706-9FA69C132F61%40flatironinstitute.org.

[Sidney Holden]

Thanks Keaton. Unfortunately I’ve found softplus to be pretty finicky here as well. 

I just realized I can solve, save, then project. That way I know exactly what’s going into the next time step, and in the post processing I just need to project K again. Good enough for now.

Cheers!

To view this discussion visit https://groups.google.com/d/msgid/dedalus-users/CADZXxBiZSAdfJLVKV9uE9gXA2Dn0bXe5V4Q-s47OSB-z1ZKJGw%40mail.gmail.com.