Help to generate code

General Usage
1

Dear community,
I think I have a good use case for the learning of metaprogramming. I want to generate this code:

        for problem in get_problems(analysis)
            for element in get_elements(problem)
                for ip in get_integration_points(element)
                    material = ip("material", time)
                    material_preprocess_increment!(material, element, ip, time)
                end
            end
        end

for example using lines:

@integrate_material material_preprocess_analysis! 
@integrate_material material_preprocess_increment!
@integrate_material material_preprocess_iteration!

The original block of code is repeated multiple times, here it is in use.

I read here: Metaprogramming · The Julia Language and here: http://mikeinnes.github.io/MacroTools.jl/stable/, but I couldn’t figure out how to start also if someone can propose the right way of working document or idea.

2

If you’re repeating code multiple times with just a different function f being called in the innermost loop, why not extract it out into a function that takes f as an argument?

4 Likes
3

Thanks. I was fixed to the macro idea and didn’t see an obvious solutions. Here is my closure:

    function operate_integration_points!(func!, time)
        for problem in get_problems(analysis)
            for element in get_elements(problem)
                for ip in get_integration_points(element)
                    if string(func!) == "FEMMaterials.material_preprocess_analysis!"
                        material_type = getfield(Materials, problem.properties.material_model)
                        material = Material(material_type, tuple())
                        ip.fields["material"] = field(material)
                        func!(material, element, ip, time)
                    else
                        material = ip("material", time)
                        func!(material, element, ip, time)
                    end
                end
            end
        end
    end

All comments welcome, I’m just learning. I have opened a pull request here: https://github.com/JuliaFEM/FEMMaterials.jl/pull/7

1 Like
4

Actually, this comparison fails in julia 1.2. What would be the right way to do it?

5

I found a solution:

func_name = last(split(string(func!),"."))
if func_name == "material_preprocess_analysis!"

It seems that julia 1.0 and 1.1 returns full name with module and julia 1.2 just the method name.

6

I think everything would be cleaner if your higher-order function was responsible only for factorizing the code that can be factorized everywhere (and not more). This would avoid the need for a test within the nested loops.

For example:

function operate_integration_points!(f, time)
    for problem in get_problems(analysis)
        for element in get_elements(problem)
            for ip in get_integration_points(element)
                f(problem, element, ip, time)
            end
        end
    end
end

If you do this, the function f passed as argument should now be responsible for determining the material from other arguments. In the specific case of material_preprocess_analysis!, this would look like

function material_preprocess_analysis_aux!(problem, element, ip, time)
    material_type = getfield(Materials, problem.properties.material_model)
    material = Material(material_type, tuple())
    ip.fields["material"] = field(material)
    material_preprocess_analysis!(material, element, ip, time)
end

operate_integration_points!(material_preprocess_analysis_aux!, time)

while in the general case, you would have:

function material_preprocess_increment_aux!(problem, element, ip, time)
    material = ip("material", time)
    material_preprocess_increment!(material, element, ip, time)
end

operate_integration_points!(material_preprocess_increment_aux!, time)

Now, instead of defining auxilliary functions, you might want to use Julia’s special do syntax. The last block of code above could be written instead like this:

operate_integration_points!(time) do problem, element, ip, time
    material = ip("material", time)
    material_preprocess_increment!(material, element, ip, time)
end

Finally, if you have many preprocessing functions which fall in the general case, and you want to factor out the logic of getting the material first, before calling the preprocessing step on it, you can use something like this:

aux(preprocess!) = (problem, element, ip, time) -> begin
    material = ip("material", time)
    preprocess!(material, element, ip, time)
end

operate_integration_points!(aux(material_preprocess_increment!), time)
operate_integration_points!(aux(material_preprocess_iteration!), time)

Does this make sense?

2 Likes
7

Yes, thanks for the good ideas. This is performance critical code and we yet don’t have performance tests. Is there any differences in performance (which you can say without testing)?

8

There should not be any performance difference.

However, thanks to BenchmarkTools.jl, such things are really easy to check by yourself. So I would advise you to

  1. hand-write write one nested loop like this (you already have it IIUC)
  2. create a realistically large test case for it (this is where you perhaps have work to do, but you’ll have to do it at some point anyway: better do it early)
  3. benchmark it (almost no work involved)
  4. transform it using higher-order functions like explained above (this is what you’re working on anyway)
  5. if benchmark performance does not degrade generalize the approach to all preprocess_*functions
1 Like
9

Thanks a lot. I am just learning the structure, but my current understanding is that all of the preprocessing functions are empty. They are generalized developer API, which gives hooks to all phases of material integration. This way material model development is disconnected from FEM solver development.

Thus I am expecting compiler to throw away most of the code most of the cases. Only challenge is that remaining code is always in the hot loop.

I hope my explanation make sense, but I got the answer I was waiting for: no differences in performance just in readability/code elegance.