Tracking/saving fields outside of those involved in a PDE

[iamlll]

Hello,

I'm working with a system of (2D spatial) equations that looks something like this:

$$ \frac{\partial u}{\partial t} = \nabla^2 v$$

$$ \frac{\partial v}{\partial t} = -\nabla^2 u $$

I'd also like to keep track of quantities such as the derivatives of these fields (e.g. d_dx(u), d_dx(v)) even though they are not explicitly included in the PDEs to be solved, but could not find any examples in the py-PDE documentation to show how to do this. Is such a task currently possible, or should calculation of gradients be done in post-processing?

[david-zwicker]

You could use the transformation argument of the tracker method of FileStorage. Alternatively, you could of course re-create the derivatives from the stored fields in post-processing.

[iamlll]

Yes - apologies, I was not clear in my previous response. I have been running
result = eq.solve(state,t_range=tmax,dt=solvestep,tracker=writer.tracker(interval=dt,transformation=fun)),
which produces errors like
ValueError: Expected 2 fields in state, but got 3 ones

[david-zwicker]

Could you provide a self-contained example showing the error? This would help me to debug things!

[iamlll]

Of course, thank you so much!!

def demo():
    L = 20
    resolution = 5
    bc = "auto_periodic_neumann"
    grid = pde.CartesianGrid([[-int(L/2), int(L/2)], [-int(L/2),int(L/2)]], [int(resolution*L), int(resolution*L)], periodic=True)

    def fun(state, t):
        div = pde.ScalarField(grid,label='div')
        div.data = state['u']._apply_operator('d_dx',bc=bc) + state['v']._apply_operator('d_dy',bc=bc)
        newstate = state.fields
        newstate.append(div)
        return pde.FieldCollection(newstate)

    expr = {
        'u': f'laplace(v)',
        'v': f'-laplace(u)',
    }
    print(expr)
    u = pde.ScalarField(grid,data=0,label='u')
    v = pde.ScalarField(grid,data=0,label='v')
    eq = pde.PDE(expr,bc=bc)
    state = pde.FieldCollection([u,v])

    path = 'demo.hdf5'
    writer = pde.FileStorage(path,write_mode='truncate')
    result = eq.solve(state,t_range=10,dt=0.005,tracker=writer.tracker(interval=0.1,transformation=fun))
    return

if __name__ == '__main__':
    demo()

[david-zwicker]

Thank you for providing the code! I can reproduce the error, which is caused by a subtle behavior of py-pde. In essence, new state.append(div) actually also modifies state.fields and leads to the error. You could instead use

newstate = state.fields[:]   # create a copy of the list
newstate.append(div)

or just

newstate = state.fields[:] + [div]

to prevent this error!

[iamlll]

That worked for me, thank you so much!