sfepy / sfepy

Main SfePy repository
http://sfepy.org
BSD 3-Clause "New" or "Revised" License
745 stars 183 forks source link

Seeking guidance on elastodynamics.py problem defnition #598

Closed ashwmk closed 4 years ago

ashwmk commented 4 years ago

I tried running simple.py on elastodynamics.py in the examples/linear_elasticity folder, and everything went fine. I checked the solution with postproc.py and the result seemed sensible. My issue is that I don't understand the equation to be solved:

equations = {
    'balance_of_forces' :  
    """dw_volume_dot.i.Omega(solid.rho, ddv, ddu)  
     + dw_zero.i.Omega(dv, du)  
     + dw_lin_elastic.i.Omega(solid.D, v, u) = 0""",  
}

First of all, the first term doesn't seem to match the one shown in the documentation: \int_{\Omega} \rho \ul{v} \pddiff{\ul{u}}{\ul{v}} Second of all, neither one of these is what I would have expected to pop out of the weak form of the wave equation: \int_{\Omega} \rho \ul{v} \ul{ddu} Finally, the equation definition doesn't even pass a unit analysis test. solid.rho, ddv, ddu, solid.D, v, and u have SI units of kg/m3, m/s2, m/s2, N/m2 (or kg/m/s2), m, and m, respectively, so the integrand for the first term ends up with units of kg/m/s4 and the integrand for the third term ends up with units of kg*m/s2.

Out of curiosity, I tried changing the first term in the equations definition to the one I thought made sense, and ran simple.py on it. The result was the following error message:

sfepy: left over: ['verbose', 'cl', 'mesh_hook', 'absolute_import', 'shape', 'cs', 'nm', 'UserMeshIO', '_filename', '__package__', 'post_process', 'lam', 'E', '__builtins__', 'H', '__file__', 'L', 'v0', 'gen_block_mesh', 'plane', 'get_ic', 'rho', '__name__', 'dt', 'dim', 'd', 'mc', 'dims', 'nu', '__doc__', 't1', 'mu']
sfepy: reading mesh [user] (function:mesh_hook)...
sfepy: ...done in 0.00 s
sfepy: creating regions...
sfepy:     Impact
sfepy:     Omega
sfepy:     Symmetry-z
sfepy:     Symmetry-y
sfepy: ...done in 0.01 s
sfepy: equation "balance_of_forces":
sfepy: dw_volume_dot.i.Omega(solid.rho, v, ddu)
     + dw_zero.i.Omega(dv, du)
     + dw_lin_elastic.i.Omega(solid.D, v, u) = 0
sfepy: using solvers:
                ts: tsn
               nls: newton
                ls: ls
sfepy: updating variables...
sfepy: ...done
sfepy: setting up dof connectivities...
sfepy: ...done in 0.00 s
sfepy: matrix shape: (6804, 6804)
sfepy: assembling matrix graph...
sfepy: ...done in 0.06 s
sfepy: matrix structural nonzeros: 421632 (9.11e-03% fill)
sfepy: updating materials...
sfepy:     solid
sfepy: ...done in 0.00 s
sfepy: ====== time 0.000000e+00 (step  1 of 19) =====
sfepy: updating variables...
sfepy: ...done
sfepy: updating materials...
sfepy:     solid
sfepy: ...done in 0.01 s
Traceback (most recent call last):
  File "simple.py", line 182, in <module>
    main()
  File "simple.py", line 179, in main
    app()
  File "C:\Python27\lib\site-packages\sfepy-2019.3-py2.7-win32.egg\sfepy\applications\application.py", line 29, in call_basic
    return self.call(**kwargs)
  File "C:\Python27\lib\site-packages\sfepy-2019.3-py2.7-win32.egg\sfepy\applications\pde_solver_app.py", line 226, in call
    post_process_hook_final=self.post_process_hook_final)
  File "C:\Python27\lib\site-packages\sfepy-2019.3-py2.7-win32.egg\sfepy\discrete\problem.py", line 1414, in solve
    status=status)
  File "C:\Python27\lib\site-packages\sfepy-2019.3-py2.7-win32.egg\sfepy\solvers\ts_solvers.py", line 34, in _standard_ts_call
    status=status, **kwargs)
  File "C:\Python27\lib\site-packages\sfepy-2019.3-py2.7-win32.egg\sfepy\solvers\ts_solvers.py", line 648, in __call__
    nls, vec0, init_fun, prestep_fun, poststep_fun)
  File "C:\Python27\lib\site-packages\sfepy-2019.3-py2.7-win32.egg\sfepy\solvers\ts_solvers.py", line 424, in get_initial_vec
    at = self.get_a0(nls, u0, v0)
  File "C:\Python27\lib\site-packages\sfepy-2019.3-py2.7-win32.egg\sfepy\solvers\ts_solvers.py", line 406, in get_a0
    a0 = nls.lin_solver(-r, mtx=M)
  File "C:\Python27\lib\site-packages\sfepy-2019.3-py2.7-win32.egg\sfepy\solvers\ls.py", line 84, in _standard_call
    context=context, **kwargs)
  File "C:\Python27\lib\site-packages\sfepy-2019.3-py2.7-win32.egg\sfepy\solvers\ls.py", line 192, in __call__
    self.presolve(mtx)
  File "C:\Python27\lib\site-packages\sfepy-2019.3-py2.7-win32.egg\sfepy\solvers\ls.py", line 203, in presolve
    self.solve = self.sls.factorized(mtx)
  File "C:\Python27\lib\site-packages\scipy\sparse\linalg\dsolve\linsolve.py", line 366, in factorized
    return splu(A).solve
  File "C:\Python27\lib\site-packages\scipy\sparse\linalg\dsolve\linsolve.py", line 242, in splu
    ilu=False, options=_options)
RuntimeError: Factor is exactly singular

I guess it's possible that I needed to change something else in the code, but in any case I'm still mystified as to where the original equation definition comes from and how it even works. Can anyone provide an explanation?

rc commented 4 years ago

The equation syntax, as it is, is the result of a design decision (not perfect), that allows solver access to the individual matrix blocks, i.e. the matrices M (mass), B (damping) and K (stiffness). Note that dv and ddv are not to be interpreted as derivatives of v, but simply as auxiliary variables - dv allows to define the damping block (with du) and ddv the inertia block (with ddu). In sfepy, each unknown variable has to have a corresponding test (virtual) variable. The elastodynamic solvers can then be abstract in a sense that they do not "see" the individual terms, and expect only the three blocks in the matrix/residual, defined by the three pairs of variables. Note that the equation in the example is not the equation that is actually solved - instead, each solver solves a set of equations corresponding to its time discretization scheme.

In retrospect, another definition might have been preferable:

equations = {
    'M' :
    """dw_volume_dot.i.Omega(solid.rho, ddv, ddu)""",
    'B' :
    """dw_zero.i.Omega(dv, du)""",
    'K' :
    """dw_lin_elastic.i.Omega(solid.D, v, u)""",
}

Does it help?

ashwmk commented 4 years ago

Thanks for the quick reply. Frankly, I'm still a little lost, but then again I'm not very well-versed with the nuts and bolts of FEM. I think, for me, a simpler follow-up question would be: where does the weak form of the equation given on the documentation page come from? Since this is just a differential equation on a continuous domain (i.e. it hasn't yet been discretized), I wouldn't expect it to be very specific to how your solvers work or what notation you use.

However, it still doesn't pass a simple unit analysis test: the first integrand has units of kg/m3 and the second integrand has units of N/m2 (if you understand u and v to be displacement vectors). Whatever method or notation you use in software, the underlying differential equation should reflect physical principals, so, at minimum, the units should work out. Why does it not work out in this case?

By the way, sorry if this discussion is going too far into the realm of FEM theory. If you could point me in the direction of any good references which are specific to the FEM algorithms implemented in SfePy, I'd really appreciate it.

rc commented 4 years ago

It's a typo in the docs, that I kept overlooking - the derivative in the first (inertia) term should be w.r.t. time (to have u''(t) = acceleration), not the test function. Thanks for spotting that! I will fix that. Let us know if you notice other problems, please.

So, in the equations string in the example, u = displacements, du = velocity, ddu = acceleration. You are right, the equation in the docs is on a continuous domain, the actual, discretized, equation depends on the time discretization scheme of a time-stepping solver (e.g. Newmark method).

As for algorithms, check our recent article [1] (or preprint [2]), and the references therein.

[1] https://link.springer.com/article/10.1007/s10444-019-09666-0 [2] https://arxiv.org/abs/1810.00674