Calculate diffusive flux from solution in interactive mode

[Scott Calabrese Barton]

I've made a lot of progress learning sfepy this weekend but I'm stuck on post-processing.  I'm in interactive mode (using Jupyter) using a toy model that solves the laplace equation on a rectangle with dirichlet conditions on the ends (image below).  Now, I'd like to evaluate the flux at designated positions, say on a line or a set of points.

I'm able to evaluate the flux field using:

flux=pb.evaluate('ev_diffusion_velocity.3.Omega(s.D, c)', copy_materials=False, mode='qp')

where 'c' is the solution field variable and s is the material. This produces a set of flux vectors at quadrature points. What I can't seem to do is put this data in a proper field variable so that I can evaluate it with probes, say for example at the boundaries. I've looked at the examples but have found them too complex to reduce to a few lines of code. Can anyone suggest the easiest way to get this data into a field variable?

[Robert Cimrman]

Hi Scott,

check the example [1] - it does pretty much all that you described, i.e.,
putting a quadrature point-based quantity into a FE field, and then probing
that field. Let us know if you need more details, or help with reducing the
code to a your problem.

r.

[1] http://sfepy.org/doc-devel/examples/linear_elasticity/its2D_interactive.html

> <https://lh3.googleusercontent.com/-IGio8BtAJJE/WOF2SY3wmCI/AAAAAAAAAKI/q9S1IPyVfp0A9FXF37NQSv_JM1qYD9LFACLcB/s1600/rectangle.jpg>
>
>

[Scott Calabrese Barton]

Robert, thank you for your response. I've now written a function that calculates the gradient as a field variable from a scalar field variable:

def gradient(pb,region,scalar,optorder=1):

    from sfepy_scb.projections import project_by_component
    import numpy as nm

    

    gfield = Field.from_args('gfield', nm.float64, 2, region,

                            approx_order=optorder - 1)

    gradient = FieldVariable('gradient', 'parameter', dfield,

                           primary_var_name='(set-to-None)')

    cfield = Field.from_args('component', nm.float64, 1, region,

                             approx_order=optorder - 1)

    component = FieldVariable('component', 'parameter', cfield,

                              primary_var_name='(set-to-None)')

        

    ev = pb.evaluate

    order = 2 * (optorder - 1)

    theterm= 'ev_grad.%d.%s(%s)' % (order,region.name,scalar.name)

    gradient_qp=ev(theterm, mode='qp')

    gradient_qp=nm.swapaxes(gradient_qp,2,3) #make the row vectors into columns

    

    project_by_component(gradient, gradient_qp, component, order)

    return gradient

I can then use the output to evaluate line probes:

gradc= gradient(pb,omega,c)

p=LineProbe([0,0],[5,0],0)

gline=p(gradc)

Note that in the second line of gradient() I've imported my own version of project_by_component. That's because the original version expects a 3-component gradient and mine has only 2 components. I had to add a check for the vector dimension:

    dim = nm.shape(tensor_qp)[2]
    for ic in range(dim):

My next task is to submit that as an issue.

I would also like to make a similar function to calculate flux, using something like:

    theterm= 'ev_diffusion_velocity.%d.%s(s.D, %s)' % (order,region.name,scalar.name)
    flux_qp=ev(theterm, copy_materials=False, mode='qp')

The problem is that I don't know how to pass the material.diffusivity in as an argument. The best is could do is pass a string 's.D' because s.D does not evaluate to anything. Any suggestions?

[Robert Cimrman]

Hi Scott,

OK, thanks for submitting the issue - it seems like an easy-to-do improvement.

> I would also like to make a similar function to calculate flux, using
> something like:
>
> theterm= 'ev_diffusion_velocity.%d.%s(s.D, %s)' % (order,region.name,
> scalar.name)
> flux_qp=ev(theterm, copy_materials=False, mode='qp')
>
> The problem is that I don't know how to pass the material.diffusivity in as
> an argument. The best is could do is pass a string 's.D' because s.D does
> not evaluate to anything. Any suggestions?
>

Either you can pass the material as a keyword argument to Problem.evaluate()
(like "s=s", with s = Material(D=...)), or you can create the term instance
directly, as in tests/test_high_level.py: test_term_evaluation() - there is
unfortunately an example for this. Let us know if you need more help - in that
case send also your script so that it can be tried.

r.