mesh-adaptation / goalie

Goal-oriented error estimation and mesh adaptation for finite element problems solved using Firedrake
Other
1 stars 1 forks source link

Misleading time in demos #141

Open ddundo opened 4 months ago

ddundo commented 4 months ago

All demos have their get_solver defined as something like this:

def get_solver(mesh_seq):
    def solver(index, ic):
        ...
        # Time integrate from t_start to t_end
        P = mesh_seq.time_partition
        t_start, t_end = P.subintervals[index]
        dt = P.timesteps[index]
        t = t_start
        while t < t_end - 1.0e-05:
            solve(F == 0, u, ad_block_tag="u")
            u_.assign(u)
            t += dt
        return {"u": u}

    return solver

where the implication is that we are solving the problem for $t \in [t{start}, t{end})$ when it's actually $t \in (t{start}, t{end}]$. This doesn't matter in these demos, but it does in my bubble shear demos, where I had to change this to

t = t_start + dt
while t < t_end + 0.5 * dt:
    # update the background velocity field at the current timestep
    u.interpolate(velocity_expression(x, y, t))

in order to compute the velocity at the correct time.

How about refactoring the demos to instead look like this:

def get_solver(mesh_seq):
    def solver(index, ic):
        ...
        # Time integrate from t_start to t_end
        P = mesh_seq.time_partition
        num_timesteps = P.num_timesteps_per_subinterval[index]
        for _ in range(num_timesteps):
            solve(F == 0, u, ad_block_tag="u")
            u_.assign(u)
        return {"u": u}

    return solver
jwallwork23 commented 4 months ago

Sure, that works! Cheers

ddundo commented 3 months ago

@jwallwork23 actually I guess a better solution would be to more rigorously define tp.subintervals so they follow this $t \in (t{start}, t{end}]$. At the moment we have tp.subintervals be something like [(0, 1), (1, 2)] whereas we would need [(0+dt, 1), (1+dt, 2)]. What do you think?

Then we can keep t = t_start in demos and similarly MeshSeq.time (introduced in #137) can be initialised to the subinterval start time. So it would be nice to keep everything consistent.

jwallwork23 commented 3 months ago

I think rather than edit the subintervals attribute (which makes intuitive sense given how the TimePartition splits up the temporal domain), it would be better to introduce two new ones (something like):

Then you could even

for time in tp.timestep_ends:
    # solve
ddundo commented 2 months ago

Thanks @jwallwork23. I'm waiting to see how the use of time develops in #28 before addressing this. (no rush at all)

Initial thought is that I don't really see a need to introduce timestep_starts and timestep_ends - I think I would rather just make a clear (verbal) distinction between tp.subintervals[0][0] and MeshSeq.time, for both of which we currently use "start time"