在 2016年8月21日星期日 UTC+2下午10:37:37，Robert Cimrman写道：
Hi Ronghai,
On 08/19/2016 07:23 PM, Ronghai Wu wrote:
Hi Robert,
I am using sfepy for static linear elasticity problem in 2D. The mesh
size
is 200*200 (mesh file is attached) and equations are as following:
def D_function(self, ts, coors, mode, term=None, **kwargs): if mode != 'qp': return _, weights = term.integral.get_qp('2_4') self.n_qp = weights.shape[0] self.val_D = np.repeat(self.stiffness_tensor, self.n_qp, axis=0) return {'function' : self.val_D}
def prestress_function(self,ts, coors, mode, term=None, **kwargs): if mode != 'qp': return _, weights = term.integral.get_qp('2_4') self.n_qp = weights.shape[0] self.val_pre = np.repeat(self.prestress_value, self.n_qp, axis=0) return {'function' : self.val_pre}
self.m = Material('m', function=self.D_function) self.prestress_value = np.copy(self.stress_0) (self.prestress_value).shape = (-1, 3, 1)
self.prestress = Material('prestress', function=self.prestress_function) integral = Integral('i', order=2) self.t1 = Term.new('dw_lin_elastic(m.function, v, u )', integral,
self.omega, m=self.m, v=self.v, u=self.u)
self.t2 = Term.new('dw_lin_prestress(prestress.function, v )',integral,
self.omega, prestress=self.prestress, v=self.v)
self.t3 = Term.new('dw_surface_ltr(traction_lef.val, v)', integral,
self.lef, traction_lef=self.traction_lef, v=self.v)
self.t4 = Term.new('dw_surface_ltr(traction_rig.val, v)', integral,
self.rig, traction_rig=self.traction_rig, v=self.v)
self.t5 = Term.new('dw_surface_ltr(traction_top.val, v)', integral,
self.top, traction_top=self.traction_top, v=self.v)
self.t6 = Term.new('dw_surface_ltr(traction_bot.val, v)', integral,
self.bot, traction_bot=self.traction_bot, v=self.v)
self.fem_eq = Equation('balance', self.t1 - self.t2 - self.t3 - self.t4
- self.t5 - self.t6)
self.eqs = Equations([self.fem_eq]) self.ls = ScipyDirect({}) self.nls_status = IndexedStruct() self.nls = Newton({'i_max' : 1,'eps_a' : 1e-8,'problem' : 'nonlinear'},
lin_solver=self.ls, status=self.nls_status)
self.pb = Problem('elasticity', equations=self.eqs, nls=self.nls, ls=
self.ls)
self.pb.time_update(ebcs=Conditions([fix_u_lef_corn]))
self.pb.update_materials() self.vec = self.pb.solve()
When I run the simulation, the information below are printed. I have several questions about it. I tried to find answers from source codes
but
failed. Could you please give me some help, thanks. (1) I do understand why "matrix shape: (79999, 79999)" and "matrix structural nonzeros: 1439965". what exactly this "matrix" is? what is
its
relation with mesh size (200*200)?
You have 201 * 201 * 2 = 80802 degrees of freedom (DOFs) - two displacements in each vertex, and some DOFs are probably fixed by fix_u_lef_corn boundary condition (what are your regions?), so the shape is (79999, 79999). Concerning the number of non-zeros, it corresponds to the mesh connectivity (how many nodes are connected to a node), and indicates how much memory the matrix storage would take.
(2) I do not understand why "updating materials..." did twice.
Because you called:
self.pb.update_materials()
Problem.update_materials() is called automatically in Problem.solve().
Thanks, that clarifies a lot.
(3) "solve: 6.18 [s]". I guess it should not take such a long time for a
mesh size (200*200) problem. Is it because I generate mesh file in the wrong way or because of other reasons?
This depends on the linear solver used. You use ScipyDirect which is either superLU (rather slow), or umfpack (faster), provided scikit-umfpack and the umfpack libraries are installed.
You can try another solver (e.g. pyamg, or an iterative solver from petsc) to speed things up.
I tried the following solvers (2D mesh 512*512), and the solve time and residual are: *ScipyDirect({'method':'superlu'}) 42 s *
*3e-12ScipyDirect({'method':'umpack'}) 7.1 s
3e-12PyAMGSolver({}) 12 s 1.e0*
*PETScKrylovSolver({}) 2.5 s 5.e-1*The *PETScKrylovSolver({}) *is the fastest, but the solution precision is low. If I try to improve the prcision by solving it iteratively with "*nls = Newton({'i_max' : 3, 'problem' : 'nonlinear'}, lin_solver=ls, status=nls_status)*", I got the following warning massage and residual remains the same:
*sfepy: warning: linear system solution precision is lowersfepy: then the value set in solver options! (err = 5.188237e-01 < 1.000000e-10)sfepy: linesearch: iter 2, (5.18824e-01 < 5.18819e-01) (new ls: 1.000000e-01)sfepy: linesearch: iter 2, (5.18824e-01 < 5.18819e-01) (new ls: 1.000000e-02)sfepy: linesearch: iter 2, (5.18824e-01 < 5.18819e-01) (new ls: 1.000000e-03)sfepy: linesearch: iter 2, (5.18824e-01 < 5.18819e-01) (new ls: 1.000000e-04)sfepy: linesearch: iter 2, (5.18824e-01 < 5.18819e-01) (new ls: 1.000000e-05)sfepy: linesearch: iter 2, (5.18824e-01 < 5.18819e-01) (new ls: 1.000000e-06)sfepy: linesearch: iter 2, (5.18824e-01 < 5.18819e-01) (new ls: 1.000000e-07)sfepy: linesearch failed, continuing anywaysfepy: nls: iter: 2, residual: 5.188237e-01 (rel: 2.155866e-02)*
But I do not understand why the iterative does not work. An additional question: is it possible to solve my case parallelly, following the example "diffusion/poisson_parallel_interactive.py".
Regards Ronghai
(4) In which source code file the following information are printed?
You can use 'git grep' to find the patterns if you have the git version of sfepy.
r.