Calculate diffusive flux from solution in interactive mode
I've made a lot of progress learning sfepy this weekend but I'm stuck on postprocessing. 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?
Hi Scott,
check the example [1]  it does pretty much all that you described, i.e., putting a quadrature pointbased 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/docdevel/examples/linear_elasticity/its2D_interactive.html
On 04/03/2017 03:12 AM, Scott Calabrese Barton wrote:
I've made a lot of progress learning sfepy this weekend but I'm stuck on postprocessing. 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, 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='(settoNone)')
cfield = Field.from_args('component', nm.float64, 1, region,
approx_order=optorder  1)
component = FieldVariable('component', 'parameter', cfield,
primary_var_name='(settoNone)')
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 3component 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?
Hi Scott,
On 9.4.2017 20:56, Scott Calabrese Barton wrote:
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='(settoNone)') cfield = Field.from_args('component', nm.float64, 1, region, approx_order=optorder  1) component = FieldVariable('component', 'parameter', cfield, primary_var_name='(settoNone)') 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 3component 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.
OK, thanks for submitting the issue  it seems like an easytodo 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.
participants (2)

Robert Cimrman

Scott Calabrese Barton