Add missing IterateOPs methods

src/floridyn_cl/iterate.jl

@@ -213,6 +213,87 @@ This is the high-performance version of the observation point iteration algorith
     return nothing
 end
 
+"""
+    iterateOPs!(::IterateOPs_average, wf::WindFarm, sim::Sim, floris::Floris, floridyn::FloriDyn, buffers::IterateOPsBuffers) -> Nothing
+
+Observation point iteration with wind field state averaging.
+
+This method mirrors the legacy MATLAB `iterateOPs` implementation where the wind
+field (`States_WF`) is advanced using an exponential moving average of the previous
+and newly shifted states controlled by weights `sim.dyn.op_iter_weights`.
+
+Weights convention:
+- `wNew = sim.dyn.op_iter_weights[1]` (weight for newly shifted state)
+- `wOld = sim.dyn.op_iter_weights[2]` (weight for existing state)
+
+The rest of the algorithm (advection, crosswind deflection, coordinate transform,
+reordering) matches the basic iteration behavior.
+"""
+@views function iterateOPs!(::IterateOPs_average, wf::WindFarm, sim::Sim, floris::Floris, floridyn::FloriDyn, 
+                             buffers::IterateOPsBuffers)
+    # Save turbine OP, turbine, and wind-field states at start indices
+    @inbounds for i in 1:size(wf.StartI, 2)
+        start_idx = wf.StartI[1, i]
+        @views buffers.tmpOPStates[i, :] .= wf.States_OP[start_idx, :]
+        @views buffers.tmpTStates[i, :]   .= wf.States_T[start_idx, :]
+        @views buffers.tmpWFStates[i, :]  .= wf.States_WF[start_idx, :]
+    end
+
+    # Downwind step (same as basic)
+    @inbounds @simd for i in 1:size(wf.States_OP, 1)
+        buffers.step_dw[i] = sim.time_step * sim.dyn.advection * wf.States_WF[i, 1]
+    end
+    @inbounds @simd for i in 1:size(wf.States_OP, 1)
+        wf.States_OP[i, 4] += buffers.step_dw[i]
+    end
+
+    # Crosswind deflection (reuse centerline!)
+    centerline!(buffers.deflection, wf.States_OP, wf.States_T, wf.States_WF, floris, wf.D[1])
+    @inbounds for i in 1:size(wf.States_OP, 1)
+        buffers.step_cw[i, 1] = buffers.deflection[i, 1] - wf.States_OP[i, 5]
+        buffers.step_cw[i, 2] = buffers.deflection[i, 2] - wf.States_OP[i, 6]
+        wf.States_OP[i, 5] = buffers.deflection[i, 1]
+        wf.States_OP[i, 6] = buffers.deflection[i, 2]
+    end
+
+    # Transform to world coordinates
+    @inbounds @simd for i in axes(wf.States_OP,1)
+        phiwi = angSOWFA2world(wf.States_WF[i, 2])
+        cphiwi = cos(phiwi)
+        sphiwi = sin(phiwi)
+        ai = buffers.step_cw[i, 1]
+        sdwi = buffers.step_dw[i]
+        wf.States_OP[i, 1] += cphiwi * sdwi - sphiwi * ai
+        wf.States_OP[i, 2] += sphiwi * sdwi + cphiwi * ai
+        wf.States_OP[i, 3] += buffers.step_cw[i, 2]
+    end
+
+    # Circular shift (manual) for OPs and restore turbine starts
+    _circshift_and_restore!(wf.States_OP, buffers.tmpOPStates, wf.StartI, buffers.temp_states_op)
+    _circshift_and_restore!(wf.States_T,  buffers.tmpTStates,  wf.StartI, buffers.temp_states_t)
+
+    # Wind field circular shift with averaging weights
+    wNew, wOld = sim.dyn.op_iter_weights  # Expect length 2 vector
+    # Perform shift into temp buffer
+    _circshift_and_restore!(wf.States_WF, buffers.tmpWFStates, wf.StartI, buffers.temp_states_wf)
+    # Average: existing (already shifted & restored) is treated as "new" part
+    @inbounds @simd for i in 1:size(wf.States_WF, 1)
+        @inbounds for j in 1:size(wf.States_WF, 2)
+            wf.States_WF[i, j] = wOld * buffers.temp_states_wf[i, j] + wNew * wf.States_WF[i, j]
+        end
+    end
+    # Restore the turbine starting rows again (MATLAB restored after averaging)
+    @inbounds for i in 1:size(wf.StartI, 2)
+        start_idx = wf.StartI[1, i]
+        @views wf.States_WF[start_idx, :] .= buffers.tmpWFStates[i, :]
+    end
+
+    # Reorder if downstream positions got unsorted
+    _reorder_ops!(wf, buffers)
+
+    return nothing
+end
+
 """
     _circshift_and_restore!(data, initial_states, start_indices, buffer)
 
@@ -281,3 +362,67 @@ function _reorder_ops!(wf::WindFarm, buffers::IterateOPsBuffers)
         end
     end
 end
+
+# function T = iterateOPs(T,Sim,paramFLORIS,paramFLORIDyn)
+# %ITERATEOPS propagates the OPs downstream and deletes the oldest ones.
+
+
+# %% Save turbine OPs
+# tmpOPStates = T.States_OP(T.StartI,:);
+# tmpTStates  = T.States_T(T.StartI,:);
+# tmpWFSTates = T.States_WF(T.StartI,:);
+# %% Shift states
+# % Calculate downwind step and apply to wake coordinate system
+# step_dw = Sim.TimeStep * T.States_WF(:,1) * Sim.Dyn.Advection;
+# T.States_OP(:,4) = T.States_OP(:,4) + step_dw;
+
+# % Calculate crosswind step and apply to wake coordinate system
+# %   TODO If turbines with different diameters are used, this has to be run
+# %   individually (in parallel) for all turbines.
+# deflection  = Centerline(...
+#     T.States_OP,T.States_T,T.States_WF,paramFLORIS,T.D(1));
+# step_cw     = (deflection - T.States_OP(:,5:6))*1;
+
+# T.States_OP(:,5:6) = deflection*1;
+
+# % Apply dw and cw step to the world coordinate system
+# phiW = angSOWFA2world(T.States_WF(:,2));
+
+# T.States_OP(:,1) = T.States_OP(:,1) + ...
+#     cos(phiW) .* step_dw - ...
+#     sin(phiW) .* step_cw(:,1);
+# T.States_OP(:,2) = T.States_OP(:,2) + ...
+#     sin(phiW) .* step_dw + ...
+#     cos(phiW) .* step_cw(:,1);
+# T.States_OP(:,3) = T.States_OP(:,3) + step_cw(:,2);
+
+# %% Circshift & Init first OPs
+# %   OPs
+# T.States_OP = circshift(T.States_OP,1);
+# T.States_OP(T.StartI,:) = tmpOPStates;
+# %   Turbines
+# T.States_T  = circshift(T.States_T,1);
+# T.States_T(T.StartI,:)  = tmpTStates;
+# %   Wind Field
+# % Getting weights for averaging
+# wNew = Sim.Dyn.OPiterWeights(1);
+# wOld = Sim.Dyn.OPiterWeights(2);
+# T.States_WF = wOld*T.States_WF + wNew*circshift(T.States_WF,1);
+# T.States_WF(T.StartI,:) = tmpWFSTates;
+
+# %% Check if OPs are in order
+# for iT=1:T.nT
+#     [~,indOP] = sort(T.States_OP(T.StartI(iT)+(0:(T.nOP-1)),4));
+#     if ~issorted(indOP)
+#         warning(['OPs overtaking, consider increasing the weight on ' ...
+#             'the old wind field state in setup > OP / wind field propagation.'])
+#         T.States_OP(T.StartI(iT)+(0:(T.nOP-1)),:) = ...
+#             T.States_OP(T.StartI(iT) + indOP - 1,:);
+#         T.States_T(T.StartI(iT)+(0:(T.nOP-1)),:) = ...
+#             T.States_T(T.StartI(iT) + indOP - 1,:);
+#         T.States_WF(T.StartI(iT)+(0:(T.nOP-1)),:) = ...
+#             T.States_WF(T.StartI(iT) + indOP - 1,:);
+#     end
+# end
+
+# end

test/runtests.jl

@@ -55,6 +55,7 @@ all_test_files = [
     "test_tit.jl",
     "test_vel.jl",
     "test_iterate.jl",
+    "test_iterate_average.jl",
     "test_floridyn_cl.jl",
     "test_runfloridyn.jl",
     "test_floris.jl",
@@ -93,6 +94,7 @@ suppress_warning_files = [
     "test_tit.jl", 
     "test_vel.jl",
     "test_iterate.jl",
+    # average iteration currently produces no warnings; keep separate if needed later
     "test_floridyn_cl.jl",
     "test_prepare_simulation.jl",
     "test_getDataVel_branches.jl",

test/test_iterate_average.jl

@@ -0,0 +1,96 @@
+# Copyright (c) 2025 Uwe Fechner
+# SPDX-License-Identifier: BSD-3-Clause
+
+using FLORIDyn, Test, LinearAlgebra
+
+"""
+    create_test_setup_average()
+
+Create a wind farm setup for testing iterateOPs! with IterateOPs_average.
+Returns (wf, sim, floris, floridyn, set, buffers).
+"""
+function create_test_setup_average()
+    settings_file = "data/2021_9T_Data.yaml"
+    wind, sim, con, floris, floridyn, ta = setup(settings_file)
+    set = Settings(wind, sim, con, Threads.nthreads() > 1, Threads.nthreads() > 1)
+    wf, wind, sim, con, floris = prepareSimulation(set, wind, con, floridyn, floris, ta, sim)
+    buffers = IterateOPsBuffers(wf)
+    return wf, sim, floris, floridyn, set, buffers
+end
+
+# Helper to set weights (mutate existing vector in-place)
+function set_iter_weights!(sim::Sim, wNew::Float64, wOld::Float64)
+    sim.dyn.op_iter_weights[1] = wNew
+    sim.dyn.op_iter_weights[2] = wOld
+    return sim
+end
+
+@testset "iterateOPs_average Unit Tests" begin
+
+    @testset "Weight Identity (old state)" begin
+        wf, sim, floris, floridyn, set, buffers = create_test_setup_average()
+        # Store original wind field
+        original_wf = copy(wf.States_WF)
+        # Set weights so result should equal old (wNew=0, wOld=1)
+        set_iter_weights!(sim, 0.0, 1.0)
+        @test_nowarn iterateOPs!(IterateOPs_average(), wf, sim, floris, floridyn, buffers)
+        @test wf.States_WF ≈ original_wf atol=1e-12
+    end
+
+    @testset "Weight Identity (shifted state matches basic)" begin
+        # Two identical setups
+        wf_a, sim_a, floris_a, floridyn_a, set_a, buffers_a = create_test_setup_average()
+        wf_b = deepcopy(wf_a); sim_b = deepcopy(sim_a); floris_b = deepcopy(floris_a); floridyn_b = floridyn_a; buffers_b = IterateOPsBuffers(wf_b)
+        # Configure weights wNew=1 (shifted), wOld=0
+        set_iter_weights!(sim_a, 1.0, 0.0)
+        iterateOPs!(IterateOPs_average(), wf_a, sim_a, floris_a, floridyn_a, buffers_a)
+        iterateOPs!(IterateOPs_basic(),   wf_b, sim_b, floris_b, floridyn_b, buffers_b)
+        @test wf_a.States_WF ≈ wf_b.States_WF atol=1e-12
+    end
+
+    @testset "General Weighted Combination" begin
+        wf, sim, floris, floridyn, set, buffers = create_test_setup_average()
+        # Capture original state
+        old = copy(wf.States_WF)
+        # Build initial_states matrix for start rows (as done in iterate)
+        nT = size(wf.StartI,2)
+        initial_states = similar(buffers.tmpWFStates)
+        for i in 1:nT
+            start_idx = wf.StartI[1,i]
+            @views initial_states[i,:] .= old[start_idx, :]
+        end
+        # Prepare a copy to produce the shifted+restored matrix S
+        shifted = copy(old)
+        temp_buffer = similar(shifted)
+        # Use internal helper (not exported)
+        FLORIDyn._circshift_and_restore!(shifted, initial_states, wf.StartI, temp_buffer)
+        # Choose non-trivial weights
+        wNew = 0.3; wOld = 0.7
+        set_iter_weights!(sim, wNew, wOld)
+        # Expected following Julia implementation order:
+        expected = wOld .* old .+ wNew .* shifted
+        # Restore start rows (Julia restores after averaging)
+        for i in 1:nT
+            start_idx = wf.StartI[1,i]
+            @views expected[start_idx, :] .= old[start_idx, :]
+        end
+        iterateOPs!(IterateOPs_average(), wf, sim, floris, floridyn, buffers)
+        @test wf.States_WF ≈ expected rtol=1e-12 atol=1e-12
+    end
+
+    @testset "Deterministic Behavior" begin
+        wf, sim, floris, floridyn, set, buffers = create_test_setup_average()
+        set_iter_weights!(sim, 0.4, 0.6)
+        results = Matrix{Float64}[]
+        for i in 1:3
+            wf_copy = deepcopy(wf)
+            sim_copy = deepcopy(sim)
+            buffers_copy = IterateOPsBuffers(wf_copy)
+            iterateOPs!(IterateOPs_average(), wf_copy, sim_copy, floris, floridyn, buffers_copy)
+            push!(results, copy(wf_copy.States_WF))
+        end
+        @test results[1] ≈ results[2]
+        @test results[2] ≈ results[3]
+    end
+end
+nothing