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
>     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
>     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-12
ScipyDirect({'method':'umpack'})   7.1 s     3e-12
PyAMGSolver({})                       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 lower
sfepy: 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 anyway
sfepy: 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.