Turbine/array modelling improvements

[cpjordan]

A few things:

[stephankramer]

Probably most importantly is to have an option to update the Physical ID on the mesh at the turbine locations so we can determine individual turbine powers more quickly.

Whatever approach you take, make sure you do some profiling on a real large scale case. I agree that it seems waste full to integrate over the entire domain, in particular if we have to do this for dozens of individual integrals, one for each turbine - but it's not necessarily the case that restricting it to a physical region id is going to make it massively faster. For example, if we have millions of nodes we'd probably be running in parallel with a fair number of cores, and the region would only be located in one of the local domain, thus the integral would only be done by one core, with all other cores idling. It might also be that we're not actually spending most of our time in the actual loop over elements, but in the "symbolic compilation and retrieving from cache" phase - and the reason it's expensive is just because of the number of integrals not the size of the mesh. If the latter is true, we could look at other ways to try to speed this up

[cpjordan]

Thanks @stephankramer as always for the useful comment. I've got (draft) PRs for the first two as they're straightforward. The last part has been a longstanding, more challenging issue which I have always dealt with by just post-processing instead...

Probably most importantly is to have an option to update the Physical ID on the mesh at the turbine locations so we can determine individual turbine powers more quickly.

Whatever approach you take, make sure you do some profiling on a real large scale case. I agree that it seems waste full to integrate over the entire domain, in particular if we have to do this for dozens of individual integrals, one for each turbine - but it's not necessarily the case that restricting it to a physical region id is going to make it massively faster.

What I tried

I generated a mesh with 235k vertices (predominantly where the array is), but with 199 subdomains for 199 turbines. I then generated a DiscreteTidalTurbineFarmOptions() for each and added these with:

for i, (phys_id, farm_options) in enumerate(zip(phys_ids, single_farm_options_list)):
    options.discrete_tidal_turbine_farms[phys_id] = [farm_options]

Checked I had my 199 turbines still (i.e. I've mapped my PhysIDs to the right turbines) and then added the callback as per usual with:

cb_turbines = turbines.TurbineFunctionalCallback(solver_obj)
solver_obj.add_callback(cb_turbines, 'timestep')

I think this is the equivalent of what I was suggesting we do (except I'd have temporarily updated the mesh to do it).

What happened

Sadly, it is taking 2300s (first iteration) / 50-130s per iteration (onward) to work as opposed to the 20s (first iteration) / 13s per iteration (onward) when we treat it as a farm using a near identical mesh but with the whole 'lease area' (100k+ vertices) as a subdomain. I was hoping the massively smaller subdomains was going to speed up the assemble() time proportionally to the inverse of the number of additional integrals I needed to do (so that even if only a limited number of cores were used, it would still be relatively fast).

... It might also be that we're not actually spending most of our time in the actual loop over elements, but in the "symbolic compilation and retrieving from cache" phase - and the reason it's expensive is just because of the number of integrals not the size of the mesh.

Judging by my results, it looks like you're probably spot on. My memory usage was quite a bit higher here too.

What can we do instead?

Before, I have just ended up extracting the entire field and processing in NumPy at each timestep i.e. just discretely integrating manually. I got close enough to the same results and in more or less the same time as the farm callback. However, this 1) is not ideal, 2) caused massive memory leakage to the extent is was not feasible and I just ended up exporting and saving the whole velocity & elevation fields for post-processing instead.

My problem/desired outcome is not necessarily speeding things up, but just wanting to get the results of individual turbines back, rather than the whole farm. Am I missing something? It isn't the calculation of the power that is time consuming, it's basically either having to a) iterate through for each turbine (slow) or b) the application of a mask such that I can see the powers of individual turbines instead of the whole farm at once (memory leakage).

That said, I am tempted to also check what speed up we can get with the farm by having the subdomain defined just over the turbine locations (but uniformly, not 1 for each), just to check if it has any significant impact.

[cpjordan]

Discussed over a call - using VertexOnlyMeshes might be the solution here