Phase field TO in gridap

New to Julia
,
1

hi, I am using gridap to do topology optimization along with phase field model. it is mentioned that the adjoint vector should be calculated in each load increament in backward manner. It is confusing for me. by this i should solve the forward solution for displacement and save each output and then again solve for adjoint in backward manner? or i can calculate in each load step while solving for displacement as below

delv = 1e-3
const vAppMax = 0.1
innerMax = 10
count = 0
p0 = ones(fem_params.np)
pf_vec = Filter(p0;r, fem_params)
pfh = FEFunction(fem_params.Pf, pf_vec)
pth = (pf -> Threshold(pf; β, η)) ∘ pfh
sPrev = CellState(1.0,dΩ)
sh = project(sPrev ,model,dΩ,order)

w = zeros(2106,730)
wg = zeros(2106,730)

while vApp .< vAppMax
    count = count .+ 1
    if vApp >= 3e-2
        delv = 1e-4
    end
    vApp = vApp .+ delv
    print("\n Entering displacemtent step :", float(vApp))
  for  inner = 1: innerMax
        ψhPlusPrev = project(ψPlusPrev ,model ,dΩ,order)
        RelErr = abs(sum(∫(Gc*ls*∇(sh)⋅∇(sh) + 2*ψhPlusPrev*sh*sh + (Gc/ls)*sh*sh)*dΩ-∫( (Gc/ls)*sh)*dΩ))/abs(sum(∫( (Gc/ls)*sh)*dΩ))
        sh_vec = stepPhaseField(fem_params,uh,ψhPlusPrev)
        sh = FEFunction(fem_params.U_PF, sh_vec)
        u_vec = stepDisp(fem_params,pth, uh, sh, vApp)
        uh = FEFunction(fem_params.U0_Disp, u_vec)
        ψhPos_in = ψPos∘(ε(uh))
        update_state!( new_EnergyState ,ψPlusPrev ,ψhPos_in)
        if   RelErr  < 1e-20
            break
        end
    
        A_mat = MatrixA(pth;fem_params)
        G_mat = MatrixG(fem_params)
        O_mat = MatrixOf(fem_params)
        w_vec =  A_mat' \ (O_mat * u_vec)
        w_vecg =  G_mat' \ (Gc * ((sh_vec .- 1)/ls) + ls*sh_vec)
        wconjh = FEFunction(fem_params.U0_Disp, conj(w_vec)) 
        wconjhg = FEFunction(fem_params.U_PF, conj(w_vecg)) 
        l_temp(dp) = ∫( (-2*(DAdpf(wconjh,uh , pfh, uh, sh; β, η)))*dp )fem_params.dΩ
        l_tempg(dpg) = ∫( -DGdpf(wconjhg,sh , pfh; β, η) *dpg)fem_params.dΩ
        dgfdpf = assemble_vector(l_temp, fem_params.Pf)
        dggfdpf = assemble_vector(l_tempg, fem_params.Pf)
        w[:,count] = dgfdpf
        wg[:,count] = dggfdpf
  end
end

thank you. i appreciate your suggestions.