Saving update state output

mary · September 18, 2023, 2:31am

Hi I am using below function to compare a scalar value and then use update state to update the inputs of function and get the updated value.

function  new_EnergyState(ψPlusPrev_in,ψhPos_in)
    ψPlus_in = ψhPos_in
    if ψPlus_in  >=ψPlusPrev_in
        ψPlus_out =ψPlus_in
    else
        ψPlus_out=ψPlusPrev_in
    end
    true,ψPlus_out
end

update_state!( new_EnergyState ,ψPlusPrev ,ψhPos_in)

I want to store the output value of update but i cant do this as the return function of new energy state is true and value. i cant delete true because it seems that update state work with that. could you please tell me how can i save the energy state?
thank you

Joris_Pinkse · September 18, 2023, 4:24am

It would help a lot if you could provide a complete example, i.e. something minimal that runs.

Sevi · September 18, 2023, 7:07am

As @Joris_Pinkse said, a bit more information would be helpful.

Are you trying to do something like this?

function calculate_something(a, b)
    # Need to return a Bool and the result of the calculation
    return true, a + b
end

# This is defined somewhere else
function update_something!(f, a, b)
    flag, value = f(a, b)
    # Do some more stuff here
    return nothing # Cannot change the return of the update function
end

for i in 1:100
    update_something!(calculate_something, rand(), rand())
end

# Now, how to access the values of `calculate_something` ?

If you can’t change the update function, then you could use a variable that holds the results and gets updated like this:

function calculate_something_and_record_value(a, b, record)
    value = a + b
    push!(record, value)
    return true, value
end

function run_updates_and_save_results()
    # A vector to hold the results
    # If we know the length and/or element type, we should specify it as well.
    record = []

    # We need to match the argument signature expected by the update function,
    # so we wrap the recording function in a "closure"
    f(a, b) = calculate_something_and_record_value(a, b, record)

    for i in 1:100
        update_something!(f, rand(), rand())
    end

    return record
end

run_updates_and_save_results()
100-element Vector{Any}:
 0.9326293845337489
 1.0874839459729055
 1.0082781204954732
 0.5836277969172434
 1.4024936915908803
 0.6865361279464467
 0.9671876419825283
 0.6145749447883307
 0.49464888052707834
 ⋮
 0.12298307120303653
 0.2654180772308613
 1.3610828094748888
 0.5013004061965702
 1.104067998491948
 1.4623574406108641
 0.9781138176451926
 0.9795816463302204
 1.3227497030387108

mary · September 18, 2023, 6:32pm

thank you for your kind explanation. unfortunately the code is not minimal to share sth that run and honestly i dont know how to make it simple to run.@Joris_Pinkse.
it is mostly like the code that you provided at first however updatestate! is a defined function in julia gridap and i didnt define it by myself anywhere. i want to save the output of update_state.

mary · September 18, 2023, 6:55pm

I did like you said but it seems that it just save the record in some steps i have provided more detailed code below

function  new_EnergyState(ψPlusPrev_in,ψhPos_in, el)
    ψPlus_in = ψhPos_in
    if ψPlus_in  >=ψPlusPrev_in
        ψPlus_out =ψPlus_in
        push!(el, ψPlus_out)
    else
        ψPlus_out=ψPlusPrev_in
        push!(el, ψPlus_out)
    end
    true,ψPlus_out
end
function   stepPhaseField(fem_params, uh_in,ΨPlusPrev_in)
    A_PF(s,φ,uh_in,ΨPlusPrev_in) =Gc*ls*∇(φ)⋅∇(s) + (2*ΨPlusPrev_in*s*φ) + (Gc/ls)*s*φ
    a_PF(s,φ) = ∫(A_PF(s,φ,uh_in,ΨPlusPrev_in))fem_params.dΩ
    b_PF(φ) =∫((Gc/ls)*φ)fem_params.dΩ
    op_PF = AffineFEOperator(a_PF,b_PF ,U_PF ,V0_PF)
    sh_out = solve(op_PF)
#     return  sh_out
    return  get_free_dof_values(sh_out)

end
using Gridap.CellData
function  stepDisp(fem_params,pth, uh_in,sh_in, vApp)
    uApp1(x) = VectorValue(0.0,0.0)
    uApp2(x) = VectorValue(0.0,0.0)
    uApp3(x) = VectorValue(0.0,-vApp)
    U_Disp = TrialFESpace(V0_Disp ,[uApp1 ,uApp2 ,uApp3])
    A_Disp(u,v,pth,uh_in,sh_in) =  ((p->Em(p))∘pth) * (ε(v) ⊙ (σfun∘(ε(u), ε(uh_in), sh_in)))
    a_Disp(u,v) =∫(A_Disp(u,v,pth,uh_in,sh_in))fem_params.dΩ
    b_Disp(v) = 0.0
    op_Disp = AffineFEOperator(a_Disp ,b_Disp ,U_Disp ,fem_params.V0_Disp)
    uh_out = solve(op_Disp)
#     return  uh_out
    return  get_free_dof_values(uh_out)

end
innerMax = 10
count = 0
vApp = 0
ψPlusPrev = CellState(0.0,dΩ)
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)


sl = zeros(sr,co)
ul = zeros(ur,co)
vl = zeros(co)
el=[]
push!(el, ψPlusPrev)
while vApp .< vAppMax
    count = count .+ 1
    if vApp >= vAppmid
        delv = delvmid
    end
    vApp = vApp .+ delv
    print("\n 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)fem_params.dΩ-∫( (Gc/ls)*sh)fem_params.dΩ))/abs(sum(∫( (Gc/ls)*sh)fem_params.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)
        uAppp1(x) = VectorValue(0.0,0.0)
        uAppp2(x) = VectorValue(0.0,0.0)
        uAppp3(x) = VectorValue(0.0,-vApp)
        U_Disp = TrialFESpace(V0_Disp ,[uAppp1 ,uAppp2 ,uAppp3])
        uh = FEFunction(U_Disp, u_vec)
        ψhPos_in = ψPos∘(ε(uh))
        update_state!( new_EnergyState ,ψPlusPrev ,ψhPos_in)
        if   RelErr  < 1e-8
            sl[:,count] = sh_vec
            ul[:,count] = u_vec
            vl[count] = vApp
            println("Energy", el)
            break
        end
   end
end
writevtk(Ω,"resultstph",cellfields=["uh"=>uh, "sh"=>sh])

and the output is like below


 step0.001EnergyAny[CellState()]

 step0.002EnergyAny[CellState()]

 step0.003EnergyAny[CellState()]

 step0.004EnergyAny[CellState()]

 step0.005EnergyAny[CellState()]

 step0.006EnergyAny[CellState()]

 step0.007EnergyAny[CellState()]

 step0.008EnergyAny[CellState()]

 step0.009000000000000001EnergyAny[CellState()]

 step0.010000000000000002EnergyAny[CellState()]

 step0.011000000000000003EnergyAny[CellState()]

 step0.012000000000000004EnergyAny[CellState()]

 step0.013000000000000005EnergyAny[CellState()]

 step0.014000000000000005EnergyAny[CellState()]

 step0.015000000000000006EnergyAny[CellState()]

 step0.016000000000000007EnergyAny[CellState()]

 step0.017000000000000008EnergyAny[CellState()]

 step0.01800000000000001EnergyAny[CellState()]

 step0.01900000000000001EnergyAny[CellState()]

 step0.02000000000000001EnergyAny[CellState()]

 step0.02100000000000001EnergyAny[CellState()]

 step0.022000000000000013EnergyAny[CellState()]

 step0.023000000000000013EnergyAny[CellState()]

 step0.024000000000000014EnergyAny[CellState()]

 step0.025000000000000015EnergyAny[CellState()]

 step0.026000000000000016EnergyAny[CellState()]

 step0.027000000000000017EnergyAny[CellState()]

 step0.028000000000000018EnergyAny[CellState()]

 step0.02900000000000002EnergyAny[CellState()]

 step0.03000000000000002EnergyAny[CellState()]

 step0.03010000000000002EnergyAny[CellState()]

 step0.03020000000000002EnergyAny[CellState()]

 step0.030300000000000018EnergyAny[CellState()]

 step0.030400000000000017EnergyAny[CellState()]

 step0.030500000000000017EnergyAny[CellState()]

 step0.030600000000000016EnergyAny[CellState()]

 step0.030700000000000015EnergyAny[CellState()]

 step0.030800000000000015EnergyAny[CellState()]

 step0.030900000000000014EnergyAny[CellState()]

 step0.031000000000000014EnergyAny[CellState()]

 step0.031100000000000013EnergyAny[CellState()]

 step0.031200000000000012EnergyAny[CellState()]

 step0.031300000000000015EnergyAny[CellState()]

 step0.03140000000000002EnergyAny[CellState()]

 step0.03150000000000002EnergyAny[CellState()]

 step0.031600000000000024EnergyAny[CellState()]

 step0.03170000000000003EnergyAny[CellState()]

 step0.03180000000000003EnergyAny[CellState()]

 step0.03190000000000003EnergyAny[CellState()]

 step0.032000000000000035EnergyAny[CellState()]

 step0.03210000000000004EnergyAny[CellState()]

 step0.03220000000000004EnergyAny[CellState()]

 step0.032300000000000044EnergyAny[CellState()]

 step0.03240000000000005EnergyAny[CellState()]

 step0.03250000000000005EnergyAny[CellState()]

 step0.03260000000000005EnergyAny[CellState()]

 step0.032700000000000055EnergyAny[CellState()]

 step0.03280000000000006EnergyAny[CellState()]

 step0.03290000000000006EnergyAny[CellState()]

 step0.033000000000000064EnergyAny[CellState()]

 step0.03310000000000007EnergyAny[CellState()]

 step0.03320000000000007EnergyAny[CellState()]

 step0.03330000000000007EnergyAny[CellState()]

 step0.033400000000000075EnergyAny[CellState()]

 step0.03350000000000008EnergyAny[CellState()]

 step0.03360000000000008EnergyAny[CellState()]

 step0.033700000000000084EnergyAny[CellState()]

 step0.03380000000000009EnergyAny[CellState()]

 step0.03390000000000009EnergyAny[CellState()]

 step0.03400000000000009EnergyAny[CellState()]

 step0.034100000000000096EnergyAny[CellState()]

 step0.0342000000000001EnergyAny[CellState()]

 step0.0343000000000001EnergyAny[CellState()]

 step0.034400000000000104EnergyAny[CellState()]

 step0.03450000000000011EnergyAny[CellState()]

 step0.03460000000000011EnergyAny[CellState()]

 step0.03470000000000011EnergyAny[CellState()]

 step0.034800000000000116
 step0.03490000000000012EnergyAny[CellState()]

 step0.03500000000000012
 step0.035100000000000124EnergyAny[CellState()]

 step0.03520000000000013
 step0.03530000000000013EnergyAny[CellState()]

 step0.03540000000000013
 step0.035500000000000136
 step0.03560000000000014EnergyAny[CellState()]

 step0.03570000000000014
 step0.035800000000000144
 step0.03590000000000015
 step0.03600000000000015
 step0.03610000000000015
 step0.036200000000000156
 step0.03630000000000016
 step0.03640000000000016
 step0.036500000000000164
 step0.03660000000000017
 step0.03670000000000017
 step0.03680000000000017
 step0.036900000000000176
 step0.03700000000000018
 step0.03710000000000018
 step0.037200000000000184
 step0.03730000000000019
 step0.03740000000000019
 step0.03750000000000019
 step0.037600000000000196
 step0.0377000000000002
 step0.0378000000000002
 step0.037900000000000204
 step0.03800000000000021
 step0.03810000000000021
 step0.03820000000000021
 step0.038300000000000216
 step0.03840000000000022
 step0.03850000000000022
 step0.038600000000000224
 step0.03870000000000023
 step0.03880000000000023
 step0.03890000000000023
 step0.039000000000000236
 step0.03910000000000024
 step0.03920000000000024

Also it seems that in this way by using record inside the calculate sth it will save record in every iteration of for loop in inner counter. I need the final value of calculate sth in each vApp
thank you

Sevi · September 18, 2023, 9:24pm

Ah, that’s very useful information. After a bit of digging, I think the function update_state! is defined here

github.com

gridap/Gridap.jl/blob/84b29ada6df077434d1b47773b30658d4c3ca85d/src/CellData/CellStates.jl#L79-L115


      
          function update_state!(updater::Function,f::CellField...)
            ids = findall(map(i->isa(i,CellState),f))
            @assert length(ids) > 0 """\n
            At least one CellState object has to be given to the update_state! function
            """
            a = f[ids]
            x = first(a).points
            @assert all(map(i->i.points===x,a)) """\n
            All the CellState objects given to the update_state! function need to be
            defined on the same CellPoint.
            """
            fx = map(i->evaluate(i,x),f)
            if num_cells(x) > 0
              fxi = map(first,fx)
              fxiq = map(first,fxi)
              need_to_update, states = first_and_tail(updater(fxiq...))
              @assert isa(need_to_update,Bool) && isa(states,Tuple{Vararg{Number}}) """\n
              Wrong return value of the user-defined updater Function. The signature is
          
                  need_to_update, state = first_and_tail(updater(args...))

This file has been truncated. show original

At the very end, it returns nothing , so the basic situation is the one I tried to describe above (plus the inner loop).

I’m not sure if the code you posted would really run, since new_EnergyState now requires three arguments, but you only supply two in update_state!. Unless you still have the original version of new_EnergyState defined (you can check by running methods(new_EnergyState) and inspecting the output), there should be a MethodError – also when you restart the Julia session. That would also explain why el is always empty, because at the moment, it never gets updated.

Note the “trick” of wrapping the three-argument function in a two-argument function using a closure:

el = [] # Some variable holding the recorded values.
new_EnergyState(a, b, record) = ... # Needs three arguments (two from the simulation, plus the record).
closure_of_new_EnergyState(a, b) = new_EnergyState(a, b, el) # Needs two arguments, so we can put it in `update_state!`. The record vector is now fixed to be the variable `el`.

If you want to save the output of new_EnergyState at the last iteration of the inner loop we need some extra bit of logic that only pushes the value to the vector if it’s the last iteration, or that overwrites the last entry if we stay in the inner loop. Here is a suggestion:

function new_EnergyState(ψPlusPrev_in, ψhPos_in)
    # Original function
end

# One function pushes a new entry ...
function new_EnergyState_push(ψPlusPrev_in, ψhPos_in, record)
    flag, value = new_EnergyState(ψPlusPrev_in, ψhPos_in)
    push!(record, value)
    return flag, value
end

# ... and another one updates the last entry
function new_EnergyState_update_last(ψPlusPrev_in, ψhPos_in, record)
    flag, value = new_EnergyState(ψPlusPrev_in, ψhPos_in)
    record[end] = value
    return flag, value
end

# ...

# Need to prepare the record variable ...
el = []
# ... as well as the closures
newE_push(a, b) = new_EnergyState_push(a, b, el)
newE_update(a, b) = new_EnergyState_update_last(a, b, el)

while ...

    for inner = 1:innerMax

        # We need to check whether to write a new entry or update the last one
        # Note that we cannot use `new_EnergyState` directly, but need to use the closures which "carry" the variable `el` with them:
        if inner == 1
            # First try: create a new entry in the record
            update_state!(newE_push, ψPlusPrev, ψhPos_in)
        else
            # We already created an entry in the record => update it
            update_state!(newE_update, ψPlusPrev, ψhPos_in)
        end

        # ...

    end

end

mary · September 20, 2023, 9:13pm

thank you for your kind explanation. It really helps

Sevi · September 21, 2023, 5:35am

You’re welcome!