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)
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 jac
Am I missing anything here (e.g., I'm not sure the BCs and ICs are taken into account)?
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?
Thanks! Dominique
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?
r.
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
On 04/21/10 00:43, Dominique wrote:
On Apr 20, 2:43 pm, Robert Cimrmancimr...@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 blocksThe 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...
Yes, your way is ok. It could be made easier though - there could be a helper function directly in sfepy.
Maybe an example of input for assemble_by_blocks() would help too:
variables = { 'u' : ('unknown field', 'displacement_Y2', 0), 'v' : ('test field', 'displacement_Y2', 'u'), 'phi' : ('unknown field', 'potential_Y2', 0), 'psi' : ('test field', 'potential_Y2', 'phi'), }
equations = { 'K,v,u' : """dw_lin_elastic.i1.Y2( inclusion.C, v, u )""",
'B,psi,u' :
"""dw_piezo_coupling.i1.Y2( inclusion.B, u, psi )""",
'D,psi,phi' :
"""dw_diffusion.i1.Y2( inclusion.D, psi, phi )""",
'M,v,u' : """dw_mass_vector.i1.Y2( inclusion.density, v, u )""",
}
... mtx = assemble_by_blocks(conf.equations, problem, ebcs = conf.ebcs, epbcs = conf.epbcs) ... mtx is then a dict with 'K', ..., 'M' keys. Of course, the corresponding fields, materials, integrals, regions etc. have to be defined properly.
Thanks for all the help. Sorry for bombarding you with questions.
You are welcome! I am glad this mailing list is alive (finally). Do not hesitate to ask, the code is not by any means perfect, and I would like to know about the various (often unexpected) ways people would like to use it.
cheers, r.
On Apr 21, 8:14 am, Robert Cimrman cimr...@ntc.zcu.cz wrote:
On 04/21/10 00:43, Dominique wrote:
On Apr 20, 2:43 pm, Robert Cimrmancimr...@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 blocksThe 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...
Yes, your way is ok. It could be made easier though - there could be a helper function directly in sfepy.
Maybe an example of input for assemble_by_blocks() would help too:
variables = { 'u' : ('unknown field', 'displacement_Y2', 0), 'v' : ('test field', 'displacement_Y2', 'u'), 'phi' : ('unknown field', 'potential_Y2', 0), 'psi' : ('test field', 'potential_Y2', 'phi'),
}
equations = { 'K,v,u' : """dw_lin_elastic.i1.Y2( inclusion.C, v, u )""",
'B,psi,u' : """dw_piezo_coupling.i1.Y2( inclusion.B, u, psi )""",
'D,psi,phi' : """dw_diffusion.i1.Y2( inclusion.D, psi, phi )""",
'M,v,u' : """dw_mass_vector.i1.Y2( inclusion.density, v, u )""",
}
... mtx = assemble_by_blocks(conf.equations, problem, ebcs = conf.ebcs, epbcs = conf.epbcs)
Ah that's it! I did find out that the variables must appear in the equation name, but I was setting my equations as actual equations (i.e., """something = something""") and that wasn't working properly for some reason.
Is there a convention in SfePy as to how the equations are to be interpreted when only their components are listed (as in your example above), i.e., does SfePy assume that the first bit is a left-hand side, the second bit is a right-hand side, and so on? Or is that left to the programmer?
mtx is then a dict with 'K', ..., 'M' keys. Of course, the corresponding fields, materials, integrals, regions etc. have to be defined properly.
Thanks for all the help. Sorry for bombarding you with questions.
You are welcome! I am glad this mailing list is alive (finally). Do not hesitate to ask, the code is not by any means perfect, and I would like to know about the various (often unexpected) ways people would like to use it.
Cheers, Dominique
On 04/21/10 09:51, Dominique wrote:
On Apr 21, 8:14 am, Robert Cimrmancimr...@ntc.zcu.cz wrote:
On 04/21/10 00:43, Dominique wrote:
On Apr 20, 2:43 pm, Robert Cimrmancimr...@ntc.zcu.cz wrote: 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 blocksThe 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...
Yes, your way is ok. It could be made easier though - there could be a helper function directly in sfepy.
Maybe an example of input for assemble_by_blocks() would help too:
variables = { 'u' : ('unknown field', 'displacement_Y2', 0), 'v' : ('test field', 'displacement_Y2', 'u'), 'phi' : ('unknown field', 'potential_Y2', 0), 'psi' : ('test field', 'potential_Y2', 'phi'),
}
equations = { 'K,v,u' : """dw_lin_elastic.i1.Y2( inclusion.C, v, u )""",
'B,psi,u' : """dw_piezo_coupling.i1.Y2( inclusion.B, u, psi )""", 'D,psi,phi' : """dw_diffusion.i1.Y2( inclusion.D, psi, phi )""", 'M,v,u' : """dw_mass_vector.i1.Y2( inclusion.density, v, u )""",}
... mtx = assemble_by_blocks(conf.equations, problem, ebcs = conf.ebcs, epbcs = conf.epbcs)
Ah that's it! I did find out that the variables must appear in the equation name, but I was setting my equations as actual equations (i.e., """something = something""") and that wasn't working properly for some reason.
Is there a convention in SfePy as to how the equations are to be interpreted when only their components are listed (as in your example above), i.e., does SfePy assume that the first bit is a left-hand side, the second bit is a right-hand side, and so on? Or is that left to the programmer?
The key should be a string with three items separated by commas: <block name>, <test variable name>, <unknown variable name>
Each weak formulation term has exactly one test variable argument, but can take several unknowns, for example (nonexistent term...)
equation = { 'A,v,u' : 'dw_a_term.i.Omega( mat, v, u, p ) 'B,v,p' : 'dw_a_term.i.Omega( mat, v, u, p ) }
would build two blocks: A = d/du dw_a_term, B = d/dp dw_a_term.
mtx is then a dict with 'K', ..., 'M' keys. Of course, the corresponding fields, materials, integrals, regions etc. have to be defined properly.
cheers, r.
participants (2)
-
Dominique -
Robert Cimrman