On Apr 20, 2:43 pm, Robert Cimrman cimr...@ntc.zcu.cz wrote:
Hi Dominique,
On 04/20/10 15:29, Dominique wrote:
Hi Robert,
From a given configuration file, I think I am able to extract the Jacobian (called the "tangent matrix" in the code) and the right-hand side as follows:
def get_problem(fName=__file__): """ Read the configuration specified in the file specified in argument (the present file by default) and return an initialized problem constructed from it along with an initial state vector. """ from sfepy.base.conf import ProblemConf, get_standard_keywords from sfepy.fem import eval_term_op, ProblemDefinition
required, other = get_standard_keywords()
# Use this file as the input file to gather configuration. conf = ProblemConf.from_file(fName, required, other)
# Create a ProblemDefinition instance from current configuration. problem = ProblemDefinition.from_conf(conf, init_variables=True, init_equations=True)
There is a shortcut for this:
problem = ProblemDefinition.from_conf_file(fName,required=required, other=other)
state = problem.create_state_vector() problem.init_variables(state) problem.time_update(None)
return (problem, state)
def get_jacobian(problem, state=None): """ Return the Jacobian matrix of the system in SciPy CSR format. Input arguments are as returned by
get_problem(). """ # Instantiate a problem evaluator. evaluator = problem.get_evaluator(mtx=problem.mtx_a) if state is None: state = problem.create_state_vector() jac = evaluator.eval_tangent_matrix(state) return jacAm I missing anything here (e.g., I'm not sure the BCs and ICs are taken into account)?
This should be cared for by the problem.time_update() call.
Now I'd like to extract the contribution of a particular equation to this Jacobian. How can I do that? This is what I've tried but it seems to return the whole Jacobian:
def get_equation_jacobian(problem, equation, state=None): """ Return the Jacobian matrix of a specific equation in SciPy CSR format. """ from sfepy.fem.evaluate import assemble_matrix evaluator = problem.get_evaluator(mtx=problem.mtx_a) if state is None: state = problem.create_state_vector() problem.variables.data_from_state(state) jac = problem.mtx_a jac = assemble_matrix(jac, equation, problem.variables, problem.materials, **evaluator.data) return jac
In fact, this function returns None, but problem.mtx_a contains the entire Jacobian. I suppose I should specify a different matrix when I instantiate the evaluator, but where can I discover its size?
The sizes of the various components can be obtained by:
(Pdb) print problem.variables.di DofInfo:dof_info details: {'p': Struct:field_var_dof_details, 'u': Struct:field_var_dof_details} indx: {'p': slice(2484, 3312, None), 'u': slice(0, 2484, None)} n_dof: {'p': 828, 'u': 2484} name: dof_info names: ['u', 'p'] ptr: [ 0 2484 3312]
problem.variables.di is the full DOF info (including BC) problem.variables.adi is the stripped DOF info (without the fixed DOFs)
Using the 'indx', you can index a submatrix.
Look also at assemble_by_blocks() in sfepy/fem/evaluate.py.
Does that help?
Yes, although it seems I'm having more success with eval_term_op as follows:
def get_jacobian_blocks(problem, state=None):
"""
Return the blocks of the Jacobian matrix as a dictionary
of sparse matrices in SciPy CSR format.
Input arguments are as returned by get_problem().
"""
from sfepy.fem.evaluate import eval_term_op
if state is None:
state = problem.create_state_vector()
blocks = {}
for key in problem.conf.equations.keys():
equation = problem.conf.equations[key]
blocks[key] = eval_term_op(state, equation, problem,
dw_mode='matrix',
tangent_matrix=None)
return blocks
The difficulty is that all blocks have the size of the full Jacobian. I guess I could extract the actual blocks that I need using the dimensions given by problem.variables.di but I suspect there's a better way...
Thanks for all the help. Sorry for bombarding you with questions.
Dominique
-- You received this message because you are subscribed to the Google Groups "sfepy-devel" group. To post to this group, send email to sfepy...@googlegroups.com. To unsubscribe from this group, send email to sfepy-devel...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sfepy-devel?hl=en.