Hi David,
could you send me the modified example?
r.
On 12/05/2017 02:04 PM, David Fernández wrote:
Hi, thank you for your reply, it has been very useful. However, there are still some things which are a bit strange to me.
1- I have modified "http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_interac..." to solve the system many times in a loop, updating only the body forces (see code bellow) and setting ls = ScipyDirect({'method':'umfpack','presolve':True}) outside the loop (and not defining nls solver). The problem is that the factorization still takes place in every iteration. I have looked for the reason and it lies in line 163 of ls.py (if self.mtx_id != id(mtx)). For some reason the system matrix is being updated even if I'm not mollifying it. By manually enforcing to no repeat the factorization (modifying ls.py) it works and it's much faster. Is it a bug or am I doing something wrong in the sample code?
2- Using nls = Newton({'is_linear' : True}, lin_solver=ls, status=nls_status) (commented lines in the sample code bellow) the results are not good (the residuals are huge and the solution are bad)
Best,
David.
#!/usr/bin/env python from __future__ import print_function from __future__ import absolute_import from argparse import ArgumentParser import numpy as nm
import sys sys.path.append('.')
from sfepy.base.base import IndexedStruct from sfepy.discrete import (FieldVariable, Material, Integral, Function, Equation, Equations, Problem) from sfepy.discrete.fem import Mesh, FEDomain, Field from sfepy.terms import Term from sfepy.discrete.conditions import Conditions, EssentialBC from sfepy.solvers.ls import ScipyDirect from sfepy.solvers.nls import Newton #from sfepy.postprocess.viewer import Viewer from sfepy.mechanics.matcoefs import stiffness_from_lame import time
def shift_u_fun(ts, coors, bc=None, problem=None, shift=0.0): """ Define a displacement depending on the y coordinate. """ val = shift * coors[:,1]**2
return val
helps = { 'show' : 'show the results figure', }
def main(): from sfepy import data_dir
parser = ArgumentParser() parser.add_argument('--version', action='version', version='%(prog)s') parser.add_argument('-s', '--show', action="store_true", dest='show', default=False, help=helps['show']) options = parser.parse_args()
mesh = Mesh.from_file(data_dir + '/meshes/2d/rectangle_tri.mesh') domain = FEDomain('domain', mesh)
min_x, max_x = domain.get_mesh_bounding_box()[:,0] eps = 1e-8 * (max_x - min_x) omega = domain.create_region('Omega', 'all') gamma1 = domain.create_region('Gamma1', 'vertices in x < %.10f' % (min_x + eps), 'facet') gamma2 = domain.create_region('Gamma2', 'vertices in x > %.10f' % (max_x - eps), 'facet')
field = Field.from_args('fu', nm.float64, 'vector', omega, approx_order=4)
u = FieldVariable('u', 'unknown', field) v = FieldVariable('v', 'test', field, primary_var_name='u')
m = Material('m', D=stiffness_from_lame(dim=2, lam=1.0, mu=1.0))
fix_u = EssentialBC('fix_u', gamma1, {'u.all' : 0.0})
bc_fun = Function('shift_u_fun', shift_u_fun, extra_args={'shift' : 0.01}) shift_u = EssentialBC('shift_u', gamma2, {'u.0' : bc_fun})
ls = ScipyDirect({'method':'umfpack','presolve':True})
t = time.time() for i in range(100):
f = Material('f', val=[[0.02+i*0.001], [0.01]])
integral = Integral('i', order=4)
t1 = Term.new('dw_lin_elastic(m.D, v, u)', integral, omega, m=m, v=v, u=u) t2 = Term.new('dw_volume_lvf(f.val, v)', integral, omega, f=f, v=v) eq = Equation('balance', t1 + t2) eqs = Equations([eq])
#nls_status = IndexedStruct() #nls = Newton({'is_linear' : True}, lin_solver=ls, status=nls_status) #pb = Problem('elasticity', equations=eqs, nls=nls, ls=ls)
pb = Problem('elasticity', equations=eqs, ls=ls)
pb.time_update(ebcs=Conditions([fix_u, shift_u])) vec = pb.solve() pb.save_state('/home/me/linear_elasticity_{}.vtk'.format(i), vec)
print(time.time()-t)
if __name__ == '__main__': main()
On 04/12/17 22:22, Robert Cimrman wrote:
Hello David,
On 12/04/2017 03:01 PM, David Fernández wrote:
Hello, I'm very new to sfepy and therefore it's being a bit hard to find in the documentation a solution to my problem, which first of all I would like to know if it is solvable within sfepy.
I'm solving the elastic equilibrium equations, which I need to solve many times during a single program run. However, in between calls to the solver, only the RHS of the FEM formulation is changing, i.e. the body forces or Neumann conditions. As in this case the stiffness matrix is constant, I would like to use a direct solver (e.g. umfpack) to calculate the matrix factorization (i.e. setup the solver) only once, an then simply call the solver many times with different body forces. I would really appreciate if some of you could point me in the right direction, or provide me with a couple of lines of code showing how this can be (more or less) achieved with this really nice module.
Yes, this could be done basically in two ways: either modify one of the non-stationary examples with a time-stepper (e.g. [1]) where you set the option
'quasistatic' : True,
for the 'ts.simple' time-stepper, and the option
'is_linear' : True,
for the 'nls.newton' nonlinear solver and the option
'presolve' : True,
for the 'ls.scipy_direct' linear solver. this will result in prefactorization of the matrix before the first time step, and reuse of the factors in all the time steps.
The other way is to take a script such as [2] and modify it as needed.
Feel free to ask more. Cheers, r.
[1] http://sfepy.org/doc-devel/examples/diffusion/time_advection_diffusion.html [2] http://sfepy.org/doc-devel/examples/linear_elasticity/linear_elastic_interac...
SfePy mailing list sfepy@python.org https://mail.python.org/mm3/mailman3/lists/sfepy.python.org/