# Phase field TO in gridap

**URL:** <https://discourse.julialang.org/t/phase-field-to-in-gridap/102607>\
**Category:** New to Julia\
**Tags:** question, gridap\
**Created:** [August 8, 2023, 2:27pm UTC](https://discourse.julialang.org/t/phase-field-to-in-gridap/102607 "2023-08-08T14:27:35Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![mary](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mary/32/53139_2.png) [@mary](https://discourse.julialang.org/u/mary)\
**Post date:** [August 8, 2023, 2:27pm UTC](https://discourse.julialang.org/t/phase-field-to-in-gridap/102607/1 "2023-08-08T14:27:35Z")

</div>

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

```julia
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.