Hi all.
I have already asked my question at http://math.stackexchange.com/questions/422928/modeling-gravity-field-with-f... and http://stackoverflow.com/questions/17160835/modeling-gravity-field-with-fini... but seems there are no people there, who knows sfepy and finite elements, in particular.
I want to model the gravity field on 3D rectangular area. It can be described by the equation: div(G) = rho. Here G is vector unknown function, rho is scalar parameter, which is constant at the selected point and fully determined by material.
The weak form is: int( div(G) * g ) = int( rho * g ). Here g is vector test function. Using sfepy syntax: {'Gravity' : """dw_div_grad.1.Omega( g, G ) = dw_volume_lvf.1.Omega( m.rho, g )"""}
I used gmsh to build .mesh file. My .geo, .mesh and .py files attached to this message.
When I try to run sfepy, it says that matrix is singular, so the solving process does not converge.
What I have to do to correct the model?
Thanks in advance.
Hi Alexander
On 06/19/2013 04:42 AM, Alexander Mihailov wrote:
Hi all.
I have already asked my question at http://math.stackexchange.com/questions/422928/modeling-gravity-field-with-f... and http://stackoverflow.com/questions/17160835/modeling-gravity-field-with-fini... but seems there are no people there, who knows sfepy and finite elements, in particular.
I want to model the gravity field on 3D rectangular area. It can be described by the equation: div(G) = rho. Here G is vector unknown function, rho is scalar parameter, which is constant at the selected point and fully determined by material.
The weak form is: int( div(G) * g ) = int( rho * g ). Here g is vector test function. Using sfepy syntax: {'Gravity' : """dw_div_grad.1.Omega( g, G ) = dw_volume_lvf.1.Omega( m.rho, g )"""}
The terms you use do not correspond to the equation. As div(G) is a scalar, it cannot be used with a vector test function, but a scalar test function - the term that does that is dw_stokes, but then the matrix would not be square.
Try using a different approach - use the scalar gravity potential u as your function, such that G = grad(u). As div(grad(u)) = rho is the Poisson equation, use dw_laplace and dw_volume_lvf terms with scalar unknown and test fields. Check the equation in [1]. Then G could be computed in a postprocess hook function by evaluating ev_grad term - see [2] for an example of a similar hook function (post_process()).
r.
[1] http://sfepy.org/doc-devel/examples/diffusion/poisson_parametric_study.html [2] http://sfepy.org/doc-devel/examples/biot/biot_npbc.html
I used gmsh to build .mesh file. My .geo, .mesh and .py files attached to this message.
When I try to run sfepy, it says that matrix is singular, so the solving process does not converge.
What I have to do to correct the model?
Thanks in advance.
Thank you, Robert.
Adding a new variable is a great way! Very likely I can just add another equation to bind u=grad(G) and G fields instead of post_process, I think that's a little bit more explicit and evident than hook magic.
field_1 = { 'name' : 'gravity', 'dtype' : nm.float64, 'shape' : (3,), 'region' : 'Omega', 'approx_order' : 1, } field_2 = { 'name' : 'gravity_potential', 'dtype' : nm.float64, 'shape' : (1,), 'region' : 'Omega', 'approx_order' : 1, } variables = { 'G' : ('unknown field', 'gravity', 1), 'g' : ('test field', 'gravity', 'G'), 'u' : ('unknown field', 'gravity_potential', 0), 'v' : ('test field', 'gravity_potential', 'u'), } equations = { 'Gravity_potential' : """dw_laplace.1.Omega( v, u ) = dw_volume_lvf.1.Omega( m.rho, v )""", 'Gravity' : """ dw_volume_dot.2.Omega( g, G ) = dw_v_dot_grad_s.2.Omega( g, u ) """, }
But here is another issue... When I run an updated task, the answer is u(x) = 0, G(x) = [0,0,0] for all x.
I think it's the wrong specifying or material, or the boundary conditions issue.
<http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6hunxpA0mh8e7...>
среда, 19 июня 2013 г., 14:11:37 UTC+7 пользователь Robert Cimrman написал:
Hi Alexander
On 06/19/2013 04:42 AM, Alexander Mihailov wrote:
Hi all.
I have already asked my question at
http://math.stackexchange.com/questions/422928/modeling-gravity-field-with-f...
and
http://stackoverflow.com/questions/17160835/modeling-gravity-field-with-fini...
seems there are no people there, who knows sfepy and finite elements, in particular.
I want to model the gravity field on 3D rectangular area. It can be described by the equation: div(G) = rho. Here G is vector unknown function, rho is scalar parameter, which is constant at the selected point and fully determined by material.
The weak form is: int( div(G) * g ) = int( rho * g ). Here g is vector test function. Using sfepy syntax: {'Gravity' : """dw_div_grad.1.Omega( g, G ) = dw_volume_lvf.1.Omega( m.rho, g )"""}
The terms you use do not correspond to the equation. As div(G) is a scalar, it cannot be used with a vector test function, but a scalar test function - the term that does that is dw_stokes, but then the matrix would not be square.
Try using a different approach - use the scalar gravity potential u as your function, such that G = grad(u). As div(grad(u)) = rho is the Poisson equation, use dw_laplace and dw_volume_lvf terms with scalar unknown and test fields. Check the equation in [1]. Then G could be computed in a postprocess hook function by evaluating ev_grad term - see [2] for an example of a similar hook function (post_process()).
r.
[1] http://sfepy.org/doc-devel/examples/diffusion/poisson_parametric_study.html [2] http://sfepy.org/doc-devel/examples/biot/biot_npbc.html
I used gmsh to build .mesh file. My .geo, .mesh and .py files attached to this message.
When I try to run sfepy, it says that matrix is singular, so the solving process does not converge.
What I have to do to correct the model?
Thanks in advance.
On 06/19/2013 11:21 AM, Alexander Mihailov wrote:
Thank you, Robert.
Adding a new variable is a great way! Very likely I can just add another equation to bind u=grad(G) and G fields instead of post_process, I think that's a little bit more explicit and evident than hook magic.
But it is not the same thing:
- the problem is more difficult to solve (1 more equation)
- the resulting G = grad(u) only in the weak sense, which is probably not what you want.
field_1 = { 'name' : 'gravity', 'dtype' : nm.float64, 'shape' : (3,), 'region' : 'Omega', 'approx_order' : 1, } field_2 = { 'name' : 'gravity_potential', 'dtype' : nm.float64, 'shape' : (1,), 'region' : 'Omega', 'approx_order' : 1, } variables = { 'G' : ('unknown field', 'gravity', 1), 'g' : ('test field', 'gravity', 'G'), 'u' : ('unknown field', 'gravity_potential', 0), 'v' : ('test field', 'gravity_potential', 'u'), } equations = { 'Gravity_potential' : """dw_laplace.1.Omega( v, u ) = dw_volume_lvf.1.Omega( m.rho, v )""", 'Gravity' : """ dw_volume_dot.2.Omega( g, G ) = dw_v_dot_grad_s.2.Omega( g, u ) """, }
But here is another issue... When I run an updated task, the answer is u(x) = 0, G(x) = [0,0,0] for all x.
I think it's the wrong specifying or material, or the boundary conditions issue.
First, check your region definitions by using:
./simple.py my.py --save-regions-as-groups --solve-not
and inspect the resulting my_regions.vtk - it seems that the regions used in the material 'm' definition are empty. Maybe you want to use 'elements of group i' selector, where i = 1, 2, 3?
r.
<http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6hunxpA0mh8e7...>
среда, 19 июня 2013 г., 14:11:37 UTC+7 пользователь Robert Cimrman написал:
Hi Alexander
On 06/19/2013 04:42 AM, Alexander Mihailov wrote:
Hi all.
I have already asked my question at
http://math.stackexchange.com/questions/422928/modeling-gravity-field-with-f...
and
http://stackoverflow.com/questions/17160835/modeling-gravity-field-with-fini...
seems there are no people there, who knows sfepy and finite elements, in particular.
I want to model the gravity field on 3D rectangular area. It can be described by the equation: div(G) = rho. Here G is vector unknown function, rho is scalar parameter, which is constant at the selected point and fully determined by material.
The weak form is: int( div(G) * g ) = int( rho * g ). Here g is vector test function. Using sfepy syntax: {'Gravity' : """dw_div_grad.1.Omega( g, G ) = dw_volume_lvf.1.Omega( m.rho, g )"""}
The terms you use do not correspond to the equation. As div(G) is a scalar, it cannot be used with a vector test function, but a scalar test function - the term that does that is dw_stokes, but then the matrix would not be square.
Try using a different approach - use the scalar gravity potential u as your function, such that G = grad(u). As div(grad(u)) = rho is the Poisson equation, use dw_laplace and dw_volume_lvf terms with scalar unknown and test fields. Check the equation in [1]. Then G could be computed in a postprocess hook function by evaluating ev_grad term - see [2] for an example of a similar hook function (post_process()).
r.
[1] http://sfepy.org/doc-devel/examples/diffusion/poisson_parametric_study.html [2] http://sfepy.org/doc-devel/examples/biot/biot_npbc.html
I used gmsh to build .mesh file. My .geo, .mesh and .py files attached to this message.
When I try to run sfepy, it says that matrix is singular, so the solving process does not converge.
What I have to do to correct the model?
Thanks in advance.
Yes, Robert, you are right, regions were empty.
After correcting regions and setting post_process_hook, I have got some solution, slowly but surely it becomes believable, but the matrix is still singular and residual is high...
sfepy: nls: iter: 0, residual: 1.528953e-02 (rel: 1.000000e+00) warning: (almost) singular matrix! (estimated cond. number: 6.06e+15) sfepy: rezidual: 0.01 [s] sfepy: solve: 0.01 [s] sfepy: matrix: 0.00 [s] sfepy: linear system not solved! (err = 1.782581e-02 < 1.000000e-10) sfepy: linesearch: iter 1, (3.34778e-02 < 1.52894e-02) (new ls: 1.000000e-01) sfepy: nls: iter: 1, residual: 1.435926e-02 (rel: 9.391562e-01) sfepy: equation "tmp": sfepy: ev_grad.1.Omega( u ) sfepy: updating materials... sfepy: ...done in 0.00 s
There are still some subtleties that need to be considered?
<http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6husEB3Ixmr7K...>
среда, 19 июня 2013 г., 16:58:44 UTC+7 пользователь Robert Cimrman написал:
On 06/19/2013 11:21 AM, Alexander Mihailov wrote:
Thank you, Robert.
Adding a new variable is a great way! Very likely I can just add another equation to bind u=grad(G) and G fields instead of post_process, I think that's a little bit more explicit and evident than hook magic.
But it is not the same thing:
- the problem is more difficult to solve (1 more equation)
- the resulting G = grad(u) only in the weak sense, which is probably not what you want.
field_1 = { 'name' : 'gravity', 'dtype' : nm.float64, 'shape' : (3,), 'region' : 'Omega', 'approx_order' : 1, } field_2 = { 'name' : 'gravity_potential', 'dtype' : nm.float64, 'shape' : (1,), 'region' : 'Omega', 'approx_order' : 1, } variables = { 'G' : ('unknown field', 'gravity', 1), 'g' : ('test field', 'gravity', 'G'), 'u' : ('unknown field', 'gravity_potential', 0), 'v' : ('test field', 'gravity_potential', 'u'), } equations = { 'Gravity_potential' : """dw_laplace.1.Omega( v, u ) = dw_volume_lvf.1.Omega( m.rho, v )""", 'Gravity' : """ dw_volume_dot.2.Omega( g, G ) = dw_v_dot_grad_s.2.Omega( g, u ) """, }
But here is another issue... When I run an updated task, the answer is u(x) = 0, G(x) = [0,0,0] for all x.
I think it's the wrong specifying or material, or the boundary conditions issue.
First, check your region definitions by using:
./simple.py my.py --save-regions-as-groups --solve-not
and inspect the resulting my_regions.vtk - it seems that the regions used in the material 'm' definition are empty. Maybe you want to use 'elements of group i' selector, where i = 1, 2, 3?
r.
< http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6hunxpA0mh8e7...>
среда, 19 июня 2013 г., 14:11:37 UTC+7 пользователь Robert Cimrman
Hi Alexander
On 06/19/2013 04:42 AM, Alexander Mihailov wrote:
Hi all.
I have already asked my question at
http://math.stackexchange.com/questions/422928/modeling-gravity-field-with-f...
and
http://stackoverflow.com/questions/17160835/modeling-gravity-field-with-fini...
seems there are no people there, who knows sfepy and finite elements, in particular.
I want to model the gravity field on 3D rectangular area. It can be described by the equation: div(G) = rho. Here G is vector unknown function, rho is scalar parameter, which is constant at the selected point and fully determined by material.
The weak form is: int( div(G) * g ) = int( rho * g ). Here g is vector test function. Using sfepy syntax: {'Gravity' : """dw_div_grad.1.Omega( g, G ) = dw_volume_lvf.1.Omega( m.rho, g )"""}
The terms you use do not correspond to the equation. As div(G) is a scalar, it cannot be used with a vector test function, but a scalar test function
написал: -
the term that does that is dw_stokes, but then the matrix would not be square.
Try using a different approach - use the scalar gravity potential u as your function, such that G = grad(u). As div(grad(u)) = rho is the Poisson equation, use dw_laplace and dw_volume_lvf terms with scalar unknown and test fields. Check the equation in [1]. Then G could be computed in a postprocess hook function by evaluating ev_grad term - see [2] for an example of a similar hook function (post_process()).
r.
[1]
http://sfepy.org/doc-devel/examples/diffusion/poisson_parametric_study.html
[2] http://sfepy.org/doc-devel/examples/biot/biot_npbc.html
I used gmsh to build .mesh file. My .geo, .mesh and .py files attached to this message.
When I try to run sfepy, it says that matrix is singular, so the solving process does not converge.
What I have to do to correct the model?
Thanks in advance.
On 06/19/2013 12:56 PM, Alexander Mihailov wrote:
Yes, Robert, you are right, regions were empty.
After correcting regions and setting post_process_hook, I have got some solution, slowly but surely it becomes believable, but the matrix is still singular and residual is high...
sfepy: nls: iter: 0, residual: 1.528953e-02 (rel: 1.000000e+00) warning: (almost) singular matrix! (estimated cond. number: 6.06e+15) sfepy: rezidual: 0.01 [s] sfepy: solve: 0.01 [s] sfepy: matrix: 0.00 [s] sfepy: linear system not solved! (err = 1.782581e-02 < 1.000000e-10) sfepy: linesearch: iter 1, (3.34778e-02 < 1.52894e-02) (new ls: 1.000000e-01) sfepy: nls: iter: 1, residual: 1.435926e-02 (rel: 9.391562e-01) sfepy: equation "tmp": sfepy: ev_grad.1.Omega( u ) sfepy: updating materials... sfepy: ...done in 0.00 s
There are still some subtleties that need to be considered?
The scalar field needs to be fixed somewhere, at least in a single vertex, as:
regions = { 'Node' : ('node 0', {}), }
ebcs = { 'fix' : ('Node', {'u.0' : 0.0}), }
Does it help?
r.
<http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6husEB3Ixmr7K...>
среда, 19 июня 2013 г., 16:58:44 UTC+7 пользователь Robert Cimrman написал:
On 06/19/2013 11:21 AM, Alexander Mihailov wrote:
Thank you, Robert.
Adding a new variable is a great way! Very likely I can just add another equation to bind u=grad(G) and G fields instead of post_process, I think that's a little bit more explicit and evident than hook magic.
But it is not the same thing:
- the problem is more difficult to solve (1 more equation)
- the resulting G = grad(u) only in the weak sense, which is probably not what you want.
field_1 = { 'name' : 'gravity', 'dtype' : nm.float64, 'shape' : (3,), 'region' : 'Omega', 'approx_order' : 1, } field_2 = { 'name' : 'gravity_potential', 'dtype' : nm.float64, 'shape' : (1,), 'region' : 'Omega', 'approx_order' : 1, } variables = { 'G' : ('unknown field', 'gravity', 1), 'g' : ('test field', 'gravity', 'G'), 'u' : ('unknown field', 'gravity_potential', 0), 'v' : ('test field', 'gravity_potential', 'u'), } equations = { 'Gravity_potential' : """dw_laplace.1.Omega( v, u ) = dw_volume_lvf.1.Omega( m.rho, v )""", 'Gravity' : """ dw_volume_dot.2.Omega( g, G ) = dw_v_dot_grad_s.2.Omega( g, u ) """, }
But here is another issue... When I run an updated task, the answer is u(x) = 0, G(x) = [0,0,0] for all x.
I think it's the wrong specifying or material, or the boundary conditions issue.
First, check your region definitions by using:
./simple.py my.py --save-regions-as-groups --solve-not
and inspect the resulting my_regions.vtk - it seems that the regions used in the material 'm' definition are empty. Maybe you want to use 'elements of group i' selector, where i = 1, 2, 3?
r.
< http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6hunxpA0mh8e7...>
среда, 19 июня 2013 г., 14:11:37 UTC+7 пользователь Robert Cimrman
Hi Alexander
On 06/19/2013 04:42 AM, Alexander Mihailov wrote:
Hi all.
I have already asked my question at
http://math.stackexchange.com/questions/422928/modeling-gravity-field-with-f...
andhttp://stackoverflow.com/questions/17160835/modeling-gravity-field-with-fini...
seems there are no people there, who knows sfepy and finite elements, in particular.
I want to model the gravity field on 3D rectangular area. It can be described by the equation: div(G) = rho. Here G is vector unknown function, rho is scalar parameter, which is constant at the selected point and fully determined by material.
The weak form is: int( div(G) * g ) = int( rho * g ). Here g is vector test function. Using sfepy syntax: {'Gravity' : """dw_div_grad.1.Omega( g, G ) = dw_volume_lvf.1.Omega( m.rho, g )"""}
The terms you use do not correspond to the equation. As div(G) is a scalar, it cannot be used with a vector test function, but a scalar test function
написал: -
the term that does that is dw_stokes, but then the matrix would not be square.
Try using a different approach - use the scalar gravity potential u as your function, such that G = grad(u). As div(grad(u)) = rho is the Poisson equation, use dw_laplace and dw_volume_lvf terms with scalar unknown and test fields. Check the equation in [1]. Then G could be computed in a postprocess hook function by evaluating ev_grad term - see [2] for an example of a similar hook function (post_process()).
r.
[1]
http://sfepy.org/doc-devel/examples/diffusion/poisson_parametric_study.html
[2] http://sfepy.org/doc-devel/examples/biot/biot_npbc.html
I used gmsh to build .mesh file. My .geo, .mesh and .py files attached to this message.
When I try to run sfepy, it says that matrix is singular, so the solving process does not converge.
What I have to do to correct the model?
Thanks in advance.
On Tue, Jun 18, 2013 at 8:42 PM, Alexander Mihailov <alexander....@gmail.com> wrote:
Hi all.
I have already asked my question at http://math.stackexchange.com/questions/422928/modeling-gravity-field-with-f... and http://stackoverflow.com/questions/17160835/modeling-gravity-field-with-fini... but seems there are no people there, who knows sfepy and finite elements, in particular.
Try to ask such questions on this site:
http://scicomp.stackexchange.com/
That's a better fit.
Ondrej
Robert,
I have set boundary conditions, but the matrix is still singular.
regions = {
'Top' : ('nodes in (z > 0.9999)', {}),
'Bottom' : ('nodes in (z < -0.9999)', {}),
'Right' : ('nodes in (y > 0.9999)', {}),
'Left' : ('nodes in (y < -0.9999)', {}),
'Far' : ('nodes in (x > 0.9999)', {}),
'Near' : ('nodes in (x < -0.9999)', {}),
}
ebcs = {
'fix0' : ('Top', {'u.0' : 0.0},),
'fix1' : ('Bottom', {'u.0' : 0.0},),
'fix2' : ('Left', {'u.0' : 0.0},),
'fix3' : ('Right', {'u.0' : 0.0},),
'fix4' : ('Far', {'u.0' : 0.0},),
'fix5' : ('Near', {'u.0' : 0.0},),
}
sfepy: nls: iter: 0, residual: 1.528953e-01 (rel: 1.000000e+00) warning: (almost) singular matrix! (estimated cond. number: 3.76e+15) sfepy: rezidual: 0.00 [s] sfepy: solve: 0.00 [s] sfepy: matrix: 0.00 [s] sfepy: nls: iter: 1, residual: 4.692475e-17 (rel: 3.069077e-16) sfepy: equation "tmp": sfepy: ev_grad.1.Omega( u ) sfepy: updating materials... sfepy: ...done in 0.00 s
среда, 19 июня 2013 г., 19:04:08 UTC+7 пользователь Robert Cimrman написал:
On 06/19/2013 12:56 PM, Alexander Mihailov wrote:
Yes, Robert, you are right, regions were empty.
After correcting regions and setting post_process_hook, I have got some solution, slowly but surely it becomes believable, but the matrix is still singular and residual is high...
sfepy: nls: iter: 0, residual: 1.528953e-02 (rel: 1.000000e+00) warning: (almost) singular matrix! (estimated cond. number: 6.06e+15) sfepy: rezidual: 0.01 [s] sfepy: solve: 0.01 [s] sfepy: matrix: 0.00 [s] sfepy: linear system not solved! (err = 1.782581e-02 < 1.000000e-10) sfepy: linesearch: iter 1, (3.34778e-02 < 1.52894e-02) (new ls: 1.000000e-01) sfepy: nls: iter: 1, residual: 1.435926e-02 (rel: 9.391562e-01) sfepy: equation "tmp": sfepy: ev_grad.1.Omega( u ) sfepy: updating materials... sfepy: ...done in 0.00 s
There are still some subtleties that need to be considered?
The scalar field needs to be fixed somewhere, at least in a single vertex, as:
regions = { 'Node' : ('node 0', {}), }
ebcs = { 'fix' : ('Node', {'u.0' : 0.0}), }
Does it help?
r.
< http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6husEB3Ixmr7K...>
среда, 19 июня 2013 г., 16:58:44 UTC+7 пользователь Robert Cimrman
написал:
On 06/19/2013 11:21 AM, Alexander Mihailov wrote:
Thank you, Robert.
Adding a new variable is a great way! Very likely I can just add another equation to bind u=grad(G) and G fields instead of post_process, I think that's a little bit more explicit and evident than hook magic.
But it is not the same thing:
- the problem is more difficult to solve (1 more equation)
- the resulting G = grad(u) only in the weak sense, which is probably
not
what you want.
field_1 = { 'name' : 'gravity', 'dtype' : nm.float64, 'shape' : (3,), 'region' : 'Omega', 'approx_order' : 1, } field_2 = { 'name' : 'gravity_potential', 'dtype' : nm.float64, 'shape' : (1,), 'region' : 'Omega', 'approx_order' : 1, } variables = { 'G' : ('unknown field', 'gravity', 1), 'g' : ('test field', 'gravity', 'G'), 'u' : ('unknown field', 'gravity_potential', 0), 'v' : ('test field', 'gravity_potential', 'u'), } equations = { 'Gravity_potential' : """dw_laplace.1.Omega( v, u ) = dw_volume_lvf.1.Omega( m.rho, v )""", 'Gravity' : """ dw_volume_dot.2.Omega( g, G ) = dw_v_dot_grad_s.2.Omega( g, u ) """, }
But here is another issue... When I run an updated task, the answer is u(x) = 0, G(x) = [0,0,0] for all x.
I think it's the wrong specifying or material, or the boundary conditions issue.
First, check your region definitions by using:
./simple.py my.py --save-regions-as-groups --solve-not
and inspect the resulting my_regions.vtk - it seems that the regions used in the material 'm' definition are empty. Maybe you want to use 'elements of group i' selector, where i = 1, 2, 3?
r.
<
http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6hunxpA0mh8e7...>
среда, 19 июня 2013 г., 14:11:37 UTC+7 пользователь Robert Cimrman
написал:
Hi Alexander
On 06/19/2013 04:42 AM, Alexander Mihailov wrote:
Hi all.
I have already asked my question at
http://math.stackexchange.com/questions/422928/modeling-gravity-field-with-f...
andhttp://stackoverflow.com/questions/17160835/modeling-gravity-field-with-fini...
the term that does that is dw_stokes, but then the matrix would not be square.
Try using a different approach - use the scalar gravity potential u as your function, such that G = grad(u). As div(grad(u)) = rho is the Poisson equation, use dw_laplace and dw_volume_lvf terms with scalar unknown and test fields. Check the equation in [1]. Then G could be computed in a postprocess hook function by evaluating ev_grad term - see [2] for an example of a similar hook function (post_process()).
r.
[1]
seems there are no people there, who knows sfepy and finite elements, in particular.
I want to model the gravity field on 3D rectangular area. It can be described by the equation: div(G) = rho. Here G is vector unknown function, rho is scalar parameter, which is constant at the selected point and fully determined by material.
The weak form is: int( div(G) * g ) = int( rho * g ). Here g is vector test function. Using sfepy syntax: {'Gravity' : """dw_div_grad.1.Omega( g, G ) = dw_volume_lvf.1.Omega( m.rho, g )"""}
The terms you use do not correspond to the equation. As div(G) is a scalar, it cannot be used with a vector test function, but a scalar test function
http://sfepy.org/doc-devel/examples/diffusion/poisson_parametric_study.html
[2] http://sfepy.org/doc-devel/examples/biot/biot_npbc.html
I used gmsh to build .mesh file. My .geo, .mesh and .py files attached to this message.
When I try to run sfepy, it says that matrix is singular, so the solving process does not converge.
What I have to do to correct the model?
Thanks in advance.
Thanks, Ondrej.
Nice board, and very relevant to my question. It's pity that I haven't found it a little earlier...
среда, 19 июня 2013 г., 22:11:09 UTC+7 пользователь Ondrej Certik написал:
On Tue, Jun 18, 2013 at 8:42 PM, Alexander Mihailov <alexan...@gmail.com <javascript:>> wrote:
Hi all.
I have already asked my question at
http://math.stackexchange.com/questions/422928/modeling-gravity-field-with-f...
and
http://stackoverflow.com/questions/17160835/modeling-gravity-field-with-fini...
but seems there are no people there, who knows sfepy and finite elements, in particular.
Try to ask such questions on this site:
http://scicomp.stackexchange.com/
That's a better fit.
Ondrej
Hmm, it might be caused by the second equation. I would really just solve the Poisson equation and used the post-processing. But the solver converged even if it complained - did you look at the solution?
r.
On 06/19/2013 05:40 PM, Alexander Mihailov wrote:
Robert,
I have set boundary conditions, but the matrix is still singular.
regions = { 'Top' : ('nodes in (z > 0.9999)', {}), 'Bottom' : ('nodes in (z < -0.9999)', {}), 'Right' : ('nodes in (y > 0.9999)', {}), 'Left' : ('nodes in (y < -0.9999)', {}), 'Far' : ('nodes in (x > 0.9999)', {}), 'Near' : ('nodes in (x < -0.9999)', {}), } ebcs = { 'fix0' : ('Top', {'u.0' : 0.0},), 'fix1' : ('Bottom', {'u.0' : 0.0},), 'fix2' : ('Left', {'u.0' : 0.0},), 'fix3' : ('Right', {'u.0' : 0.0},), 'fix4' : ('Far', {'u.0' : 0.0},), 'fix5' : ('Near', {'u.0' : 0.0},), }sfepy: nls: iter: 0, residual: 1.528953e-01 (rel: 1.000000e+00) warning: (almost) singular matrix! (estimated cond. number: 3.76e+15) sfepy: rezidual: 0.00 [s] sfepy: solve: 0.00 [s] sfepy: matrix: 0.00 [s] sfepy: nls: iter: 1, residual: 4.692475e-17 (rel: 3.069077e-16) sfepy: equation "tmp": sfepy: ev_grad.1.Omega( u ) sfepy: updating materials... sfepy: ...done in 0.00 s
среда, 19 июня 2013 г., 19:04:08 UTC+7 пользователь Robert Cimrman написал:
On 06/19/2013 12:56 PM, Alexander Mihailov wrote:
Yes, Robert, you are right, regions were empty.
After correcting regions and setting post_process_hook, I have got some solution, slowly but surely it becomes believable, but the matrix is still singular and residual is high...
sfepy: nls: iter: 0, residual: 1.528953e-02 (rel: 1.000000e+00) warning: (almost) singular matrix! (estimated cond. number: 6.06e+15) sfepy: rezidual: 0.01 [s] sfepy: solve: 0.01 [s] sfepy: matrix: 0.00 [s] sfepy: linear system not solved! (err = 1.782581e-02 < 1.000000e-10) sfepy: linesearch: iter 1, (3.34778e-02 < 1.52894e-02) (new ls: 1.000000e-01) sfepy: nls: iter: 1, residual: 1.435926e-02 (rel: 9.391562e-01) sfepy: equation "tmp": sfepy: ev_grad.1.Omega( u ) sfepy: updating materials... sfepy: ...done in 0.00 s
There are still some subtleties that need to be considered?
The scalar field needs to be fixed somewhere, at least in a single vertex, as:
regions = { 'Node' : ('node 0', {}), }
ebcs = { 'fix' : ('Node', {'u.0' : 0.0}), }
Does it help?
r.
< http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6husEB3Ixmr7K...>
среда, 19 июня 2013 г., 16:58:44 UTC+7 пользователь Robert Cimrman
написал:
On 06/19/2013 11:21 AM, Alexander Mihailov wrote:
Thank you, Robert.
Adding a new variable is a great way! Very likely I can just add another equation to bind u=grad(G) and G fields instead of post_process, I think that's a little bit more explicit and evident than hook magic.
But it is not the same thing:
- the problem is more difficult to solve (1 more equation)
- the resulting G = grad(u) only in the weak sense, which is probably
not
what you want.
field_1 = { 'name' : 'gravity', 'dtype' : nm.float64, 'shape' : (3,), 'region' : 'Omega', 'approx_order' : 1, } field_2 = { 'name' : 'gravity_potential', 'dtype' : nm.float64, 'shape' : (1,), 'region' : 'Omega', 'approx_order' : 1, } variables = { 'G' : ('unknown field', 'gravity', 1), 'g' : ('test field', 'gravity', 'G'), 'u' : ('unknown field', 'gravity_potential', 0), 'v' : ('test field', 'gravity_potential', 'u'), } equations = { 'Gravity_potential' : """dw_laplace.1.Omega( v, u ) = dw_volume_lvf.1.Omega( m.rho, v )""", 'Gravity' : """ dw_volume_dot.2.Omega( g, G ) = dw_v_dot_grad_s.2.Omega( g, u ) """, }
But here is another issue... When I run an updated task, the answer is u(x) = 0, G(x) = [0,0,0] for all x.
I think it's the wrong specifying or material, or the boundary conditions issue.
First, check your region definitions by using:
./simple.py my.py --save-regions-as-groups --solve-not
and inspect the resulting my_regions.vtk - it seems that the regions used in the material 'm' definition are empty. Maybe you want to use 'elements of group i' selector, where i = 1, 2, 3?
r.
<
http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6hunxpA0mh8e7...>
среда, 19 июня 2013 г., 14:11:37 UTC+7 пользователь Robert Cimrman
написал:
Hi Alexander
On 06/19/2013 04:42 AM, Alexander Mihailov wrote: > Hi all. > > I have already asked my question at >
http://math.stackexchange.com/questions/422928/modeling-gravity-field-with-f...
> and >
http://stackoverflow.com/questions/17160835/modeling-gravity-field-with-fini...
the term that does that is dw_stokes, but then the matrix would not be square.
Try using a different approach - use the scalar gravity potential u as your function, such that G = grad(u). As div(grad(u)) = rho is the Poisson equation, use dw_laplace and dw_volume_lvf terms with scalar unknown and test fields. Check the equation in [1]. Then G could be computed in a postprocess hook function by evaluating ev_grad term - see [2] for an example of a similar hook function (post_process()).
r.
[1]
> seems there are no people there, who knows sfepy and finite elements, in > particular. > > I want to model the gravity field on 3D rectangular area. It can be > described by the equation: div(G) = rho. Here G is vector unknown function, > rho is scalar parameter, which is constant at the selected point and fully > determined by material. > > The weak form is: int( div(G) * g ) = int( rho * g ). Here g is vector test > function. > Using sfepy syntax: {'Gravity' : """dw_div_grad.1.Omega( g, G ) = > dw_volume_lvf.1.Omega( m.rho, g )"""}
The terms you use do not correspond to the equation. As div(G) is a scalar, it cannot be used with a vector test function, but a scalar test function
http://sfepy.org/doc-devel/examples/diffusion/poisson_parametric_study.html
[2] http://sfepy.org/doc-devel/examples/biot/biot_npbc.html
> I used gmsh to build .mesh file. My .geo, .mesh and .py files attached to > this message. > > When I try to run sfepy, it says that matrix is singular, so the solving > process does not converge. > > What I have to do to correct the model? > > Thanks in advance. >
Yes, I already using postprocessing to get the vector field (though it gives the sufficiently similar results as the second equation, I checked).
I reduced this problem to 2D case. I looked at the solution .vtk file and noticed that the computed gravity potential in the brick is extremely large (-1.508473e+13), so maybe it is something wrong here. The remaining field looks more or less adequate for me, but the field behaviour in general is strange, too:
<http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6hvOKiDh7XgW1...> Can it be the issue of discontinuous material coefficient? Or maybe improperly generated mesh?
среда, 19 июня 2013 г., 23:16:11 UTC+7 пользователь Robert Cimrman написал:
Hmm, it might be caused by the second equation. I would really just solve the Poisson equation and used the post-processing. But the solver converged even if it complained - did you look at the solution?
r.
On 06/19/2013 05:40 PM, Alexander Mihailov wrote:
Robert,
I have set boundary conditions, but the matrix is still singular.
regions = { 'Top' : ('nodes in (z > 0.9999)', {}), 'Bottom' : ('nodes in (z < -0.9999)', {}), 'Right' : ('nodes in (y > 0.9999)', {}), 'Left' : ('nodes in (y < -0.9999)', {}), 'Far' : ('nodes in (x > 0.9999)', {}), 'Near' : ('nodes in (x < -0.9999)', {}), } ebcs = { 'fix0' : ('Top', {'u.0' : 0.0},), 'fix1' : ('Bottom', {'u.0' : 0.0},), 'fix2' : ('Left', {'u.0' : 0.0},), 'fix3' : ('Right', {'u.0' : 0.0},), 'fix4' : ('Far', {'u.0' : 0.0},), 'fix5' : ('Near', {'u.0' : 0.0},), }sfepy: nls: iter: 0, residual: 1.528953e-01 (rel: 1.000000e+00) warning: (almost) singular matrix! (estimated cond. number: 3.76e+15) sfepy: rezidual: 0.00 [s] sfepy: solve: 0.00 [s] sfepy: matrix: 0.00 [s] sfepy: nls: iter: 1, residual: 4.692475e-17 (rel: 3.069077e-16) sfepy: equation "tmp": sfepy: ev_grad.1.Omega( u ) sfepy: updating materials... sfepy: ...done in 0.00 s
среда, 19 июня 2013 г., 19:04:08 UTC+7 пользователь Robert Cimrman написал:
On 06/19/2013 12:56 PM, Alexander Mihailov wrote:
Yes, Robert, you are right, regions were empty.
After correcting regions and setting post_process_hook, I have got
some
solution, slowly but surely it becomes believable, but the matrix is still singular and residual is high...
sfepy: nls: iter: 0, residual: 1.528953e-02 (rel: 1.000000e+00) warning: (almost) singular matrix! (estimated cond. number: 6.06e+15) sfepy: rezidual: 0.01 [s] sfepy: solve: 0.01 [s] sfepy: matrix: 0.00 [s] sfepy: linear system not solved! (err = 1.782581e-02 < 1.000000e-10) sfepy: linesearch: iter 1, (3.34778e-02 < 1.52894e-02) (new ls: 1.000000e-01) sfepy: nls: iter: 1, residual: 1.435926e-02 (rel: 9.391562e-01) sfepy: equation "tmp": sfepy: ev_grad.1.Omega( u ) sfepy: updating materials... sfepy: ...done in 0.00 s
There are still some subtleties that need to be considered?
The scalar field needs to be fixed somewhere, at least in a single vertex, as:
regions = { 'Node' : ('node 0', {}), }
ebcs = { 'fix' : ('Node', {'u.0' : 0.0}), }
Does it help?
r.
<
http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6husEB3Ixmr7K...>
среда, 19 июня 2013 г., 16:58:44 UTC+7 пользователь Robert Cimrman
написал:
On 06/19/2013 11:21 AM, Alexander Mihailov wrote:
Thank you, Robert.
Adding a new variable is a great way! Very likely I can just add another equation to bind u=grad(G) and G fields instead of post_process, I think that's a little bit more explicit
evident than hook magic.
But it is not the same thing:
- the problem is more difficult to solve (1 more equation)
- the resulting G = grad(u) only in the weak sense, which is
and probably
not
what you want.
field_1 = { 'name' : 'gravity', 'dtype' : nm.float64, 'shape' : (3,), 'region' : 'Omega', 'approx_order' : 1, } field_2 = { 'name' : 'gravity_potential', 'dtype' : nm.float64, 'shape' : (1,), 'region' : 'Omega', 'approx_order' : 1, } variables = { 'G' : ('unknown field', 'gravity', 1), 'g' : ('test field', 'gravity', 'G'), 'u' : ('unknown field', 'gravity_potential', 0), 'v' : ('test field', 'gravity_potential', 'u'), } equations = { 'Gravity_potential' : """dw_laplace.1.Omega( v, u ) = dw_volume_lvf.1.Omega( m.rho, v )""", 'Gravity' : """ dw_volume_dot.2.Omega( g, G ) = dw_v_dot_grad_s.2.Omega( g, u ) """, }
But here is another issue... When I run an updated task, the answer is u(x) = 0, G(x) = [0,0,0] for all x.
I think it's the wrong specifying or material, or the boundary conditions issue.
First, check your region definitions by using:
./simple.py my.py --save-regions-as-groups --solve-not
and inspect the resulting my_regions.vtk - it seems that the regions used in the material 'm' definition are empty. Maybe you want to use 'elements of group i' selector, where i = 1, 2, 3?
r.
<
http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6hunxpA0mh8e7...>
среда, 19 июня 2013 г., 14:11:37 UTC+7 пользователь Robert Cimrman
написал:
> > Hi Alexander > > On 06/19/2013 04:42 AM, Alexander Mihailov wrote: >> Hi all. >> >> I have already asked my question at >> >
http://math.stackexchange.com/questions/422928/modeling-gravity-field-with-f...
>> and >> >
>> seems there are no people there, who knows sfepy and finite elements, in >> particular. >> >> I want to model the gravity field on 3D rectangular area. It can be >> described by the equation: div(G) = rho. Here G is vector unknown > function, >> rho is scalar parameter, which is constant at the selected point and > fully >> determined by material. >> >> The weak form is: int( div(G) * g ) = int( rho * g ). Here g is vector > test >> function. >> Using sfepy syntax: {'Gravity' : """dw_div_grad.1.Omega( g, G ) = >> dw_volume_lvf.1.Omega( m.rho, g )"""} > > The terms you use do not correspond to the equation. As div(G) is a > scalar, it > cannot be used with a vector test function, but a scalar test function
> the > term that does that is dw_stokes, but then the matrix would not be square. > > Try using a different approach - use the scalar gravity potential u as > your > function, such that G = grad(u). As div(grad(u)) = rho is the Poisson > equation, > use dw_laplace and dw_volume_lvf terms with scalar unknown and test > fields. > Check the equation in [1]. Then G could be computed in a
http://stackoverflow.com/questions/17160835/modeling-gravity-field-with-fini... postprocess
hook
> function by evaluating ev_grad term - see [2] for an example of a similar > hook > function (post_process()). > > r. > > [1] >
http://sfepy.org/doc-devel/examples/diffusion/poisson_parametric_study.html
> [2] http://sfepy.org/doc-devel/examples/biot/biot_npbc.html > >> I used gmsh to build .mesh file. My .geo, .mesh and .py files attached > to >> this message. >> >> When I try to run sfepy, it says that matrix is singular, so the solving >> process does not converge. >> >> What I have to do to correct the model? >> >> Thanks in advance. >> > >
I got it!
<http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6hvU8yZbxn7ko...>
That really was a bad mesh file: elements of brick and the remaining ones were inconsistent.
Thanks all! :)
On 06/20/2013 11:28 AM, Alexander Mihailov wrote:
I got it!
<http://ft.trillian.im/b820c2e0aaf525034498f9740f19cbbd7b172056/6hvU8yZbxn7ko...>
That really was a bad mesh file: elements of brick and the remaining ones were inconsistent.
Thanks all! :)
Glad to hear that!
hth, r.
participants (3)
-
Alexander Mihailov -
Ondřej Čertík -
Robert Cimrman