Re: [deal.II] How to sum up several LinearAlgebraPETSc::MPI::Vector completely_distributed_solution ?
Dear Prof. Wolfgang Bangerth, I have tried with the sum of the right hand side, but as I think, it did not work for my problem. I will try with your guidance. Thank you very much. Best regards, Kien Vào Th 7, 3 thg 8, 2019 vào lúc 09:42 Wolfgang Bangerth < bange...@colostate.edu> đã viết: > On 8/1/19 9:24 PM, Phạm Ngọc Kiên wrote: > > > > However, I do not know if I can sum up all the > completely_distributed_solution > > from the different right-hand-sides in order to get a vector of solution. > > Could you please tell me how to do this? > > You could just sum them up: > >AppropriateType sum_of_solutions (...); >for (unsigned int i=0; i<...; ++i) > sum_of_solutions += solution[i]; > > Of course, the sum of the solutions equals the solution of a single linear > system with the sum of the right hand sides and the same matrix -- which > would > by substantially cheaper to compute. > > Best > W. > > -- > > Wolfgang Bangerth email: bange...@colostate.edu > www: http://www.math.colostate.edu/~bangerth/ > > -- > The deal.II project is located at http://www.dealii.org/ > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > --- > You received this message because you are subscribed to the Google Groups > "deal.II User Group" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to dealii+unsubscr...@googlegroups.com. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/2871cd63-7236-a03d-cb02-8dc49a98a8db%40colostate.edu > . > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/CAAo%2BSZeZDo0dLUfG-xtvi1XdRi-BkJb5N0wRCowK52ZrUAa9EA%40mail.gmail.com.
Re: [deal.II] How to sum up several LinearAlgebraPETSc::MPI::Vector completely_distributed_solution ?
On 8/1/19 9:24 PM, Phạm Ngọc Kiên wrote: > > However, I do not know if I can sum up all the > completely_distributed_solution > from the different right-hand-sides in order to get a vector of solution. > Could you please tell me how to do this? You could just sum them up: AppropriateType sum_of_solutions (...); for (unsigned int i=0; i<...; ++i) sum_of_solutions += solution[i]; Of course, the sum of the solutions equals the solution of a single linear system with the sum of the right hand sides and the same matrix -- which would by substantially cheaper to compute. Best W. -- Wolfgang Bangerth email: bange...@colostate.edu www: http://www.math.colostate.edu/~bangerth/ -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/2871cd63-7236-a03d-cb02-8dc49a98a8db%40colostate.edu.
Re: [deal.II] Collapse vertex and apply constraint
On Tuesday, July 30, 2019 at 10:06:59 PM UTC+2, Wolfgang Bangerth wrote:
>
>
> You can't compute the determinant for anything but quite small matrices.
> It's not a numerically stable operation. The only way to determine
> whether a matrix is singular is to compute all eigenvalues and see
> whether one or more are "small" compared to the size of the other
> eigenvalues.
>
> Can you do that for the matrices you have, and plot the eigenvalues in
> the same plot? Is there a qualitative difference?
>
>
>
Hi Wolfgang,
I calculated the eigenvalues and the condition number. It seems that it is
ok to collapse one vertex and transform the quad into a triangle.
Interestingly, applying the constraints with AffineConstraints does not
change the condition number. Even, if I fix all the dofs with add_line()
the condition number does not change.
I collapsed one vertex in a 3D simulation with NedelecSZ and the solution
and the condition are ok. I have to do more tests with 3D simulations with
NelelecSZ.
In the simple simulation from step-6. I collapsed one vertex and then, if I
move another vertex and squeeze the quad into almost a line the surface
becomes very small and then the condition number becomes very large. Note
that I didn't completely squeeze the quad into a line. If the surface is
zero, then as expected the simulation breaks.
This means that it is ok to collapse a vertex as long as the surface/volume
does not become too small?
According to this publication it is ok to collapse vertices in quads/hexes
https://www.sciencedirect.com/science/article/pii/S1877705814016713
Enclosed you can find my program and the python script to analyze the data.
To run it: ./collapse_vertex>matrix.txt
You can also find enclosed the meshes. I added the condition number in the
image for each mesh.
Below you can find the detailed results :
- Original grid:
* surface_ratio = 1
* eigenvalues = [ 0.67 0.67 0.67 0.67
1.17200912 1.3
1.3 1.3 1.3 1.3 1.3 1.3
1.3 1.3 1.3 1.3 1.3
9.13881996
9.13881996 9.42425 12.79450083 21.95213004 21.95213004 31.83336
58.59408006]
* condition_number = 80
- Collapse vertex a
* surface_ratio = 0.5
* eigenvalues = [ 0.67 0.67 0.67 0.67
1.17122833
1.3 1.3 1.3 1.3 1.3
1.3 1.3 1.3 1.3 1.3
1.3 1.3 8.92744864 10.04655013 12.5454928
13.98137723 20.47958898 31.17725707 45.95448961 185.24586722]
* condition_number = 3e+02
- Collapse vertex a, move vertex b
* surface_ratio = 0.1
* eigenvalues = [ 0.67 0.67 0.67 0.67
1.16902972
1.3 1.3 1.3 1.3 1.3
1.3 1.3 1.3 1.3 1.33939
1.52754 1.53365.11586123 10.46272767 12.09280555
14.54544349 22.42497592 31.57412136 50.19072911 264.78640594]
* condition_number = 4e+02
-Collapse vertex a, almost collapse vertex b
* surface_ratio = 1e-8
* eigenvalues = [ 6.7000e-01 6.7000e-01 6.7000e-01
6.7000e-01
1.21307648e+00 1.3000e+00 1.3000e+00 1.3000e+00
1.3000e+00 1.3000e+00 1.3000e+00 1.3000e+00
1.3000e+00 1.3000e+00 1.35897000e+00 1.6000e+00
1.69231000e+00 8.83186491e+00 1.10689839e+01 1.37682911e+01
1.98164659e+01 2.84090528e+01 4.79480606e+01 2.92923117e+05
1.70718688e+06]
* condition_number = 3e+06
Best,
Dani
[image: original_grid.png][image: collapse_a.png][image:
collapse_a_move_b.png][image: collapse_a_close_to_collapse_b.png]
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/c680299c-8124-4b8a-94fd-614a6e8f0017%40googlegroups.com.
#!/usr/bin/env python3
import numpy as np
import numpy.linalg
matrix_size = 25
matrix = np.zeros((matrix_size,matrix_size))
hola_txt = open('build/matrix.txt','r')
for line in hola_txt:
if line[0] == '(':
value = float(line.split(' ')[1])
coordinates = line.split(' ')[0]
x = int(coordinates.split(',')[0][1:])
y = int(coordinates.split(',')[1][:-1])
matrix[x,y] = value
determinant = np.linalg.det(matrix)
eigenvalues, eigenvectors = np.linalg.eig(matrix)
condition_number = np.linalg.cond(matrix)
print('{:.2e}'.format(determinant))
print(np.sort(eigenvalues))
Re: [deal.II] Collapse vertex and apply constraint
On Tuesday, July 30, 2019 at 10:06:59 PM UTC+2, Wolfgang Bangerth wrote:
>
>
>
> You can't compute the determinant for anything but quite small matrices.
> It's not a numerically stable operation. The only way to determine
> whether a matrix is singular is to compute all eigenvalues and see
> whether one or more are "small" compared to the size of the other
> eigenvalues.
>
> Can you do that for the matrices you have, and plot the eigenvalues in
> the same plot? Is there a qualitative difference?
>
>
>
Hi Wolfgang,
I calculated the eigenvalues and the condition number. It seems that it is
ok to collapse one vertex and transform the quad into a triangle.
Interestingly, applying the constraints with AffineConstraint does not
change the condition number. Even if I fix all the dofs with add_line() the
condition number doesn't change.
I collapsed one vertex in a 3D simulation with NedelecSZ and the solution
and the condition number are ok. I have to do more tests with NelelecSZ.
In the simple simulation from step-6. I collapse one vertex and then if I
move another vertex and transform the quad into a line the surface becomes
very small and the condition number becomes very large. Note that I didn't
completely squeeze the quad into a line. If the surface is zero, then as
expected the simulation breaks.
This means that it is ok to collapse a vertex as long as the surface/volume
does not become too small?
According to this publication it is ok to collapse vertices
https://www.sciencedirect.com/science/article/pii/S1877705814016713
Enclosed you can find my program and the python script to analyze the data.
To run it: ./collapse_vertex>matrix.txt
You can also find enclosed the meshes. I added the condition number in the
image for each mesh.
Below you can find the detailed results :
- Original grid:
* surface_ratio = 1
* eigenvalues = [ 0.67 0.67 0.67 0.67
1.17200912 1.3
1.3 1.3 1.3 1.3 1.3 1.3
1.3 1.3 1.3 1.3 1.3
9.13881996
9.13881996 9.42425 12.79450083 21.95213004 21.95213004 31.83336
58.59408006]
* condition_number = 80
- Collapse vertex a
* surface_ratio = 0.5
* eigenvalues = [ 0.67 0.67 0.67 0.67
1.17122833
1.3 1.3 1.3 1.3 1.3
1.3 1.3 1.3 1.3 1.3
1.3 1.3 8.92744864 10.04655013 12.5454928
13.98137723 20.47958898 31.17725707 45.95448961 185.24586722]
* condition_number = 3e+02
- Collapse vertex a, move vertex b
* surface_ratio = 0.1
* eigenvalues = [ 0.67 0.67 0.67 0.67
1.16902972
1.3 1.3 1.3 1.3 1.3
1.3 1.3 1.3 1.3 1.33939
1.52754 1.53365.11586123 10.46272767 12.09280555
14.54544349 22.42497592 31.57412136 50.19072911 264.78640594]
* condition_number = 4e+02
-Collapse vertex a, almost collapse vertex b
* surface_ratio = 1e-8
* eigenvalues = [ 6.7000e-01 6.7000e-01 6.7000e-01
6.7000e-01
1.21307648e+00 1.3000e+00 1.3000e+00 1.3000e+00
1.3000e+00 1.3000e+00 1.3000e+00 1.3000e+00
1.3000e+00 1.3000e+00 1.35897000e+00 1.6000e+00
1.69231000e+00 8.83186491e+00 1.10689839e+01 1.37682911e+01
1.98164659e+01 2.84090528e+01 4.79480606e+01 2.92923117e+05
1.70718688e+06]
* condition_number = 3e+06
Best,
Dani
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/6b07998a-f493-4b4d-9373-2a01333585f3%40googlegroups.com.
#!/usr/bin/env python3
import numpy as np
import numpy.linalg
matrix_size = 25
matrix = np.zeros((matrix_size,matrix_size))
hola_txt = open('build/matrix.txt','r')
for line in hola_txt:
if line[0] == '(':
value = float(line.split(' ')[1])
coordinates = line.split(' ')[0]
x = int(coordinates.split(',')[0][1:])
y = int(coordinates.split(',')[1][:-1])
matrix[x,y] = value
determinant = np.linalg.det(matrix)
eigenvalues, eigenvectors = np.linalg.eig(matrix)
condition_number = np.linalg.cond(matrix)
print('{:.2e}'.format(determinant))
print(np.sort(eigenvalues))
print('{:.2e}'.format(condition_number))
/* -
*
* Copyright (C) 2000 - 2019 by the deal.II authors
*
* This
[deal.II] Creation of diagonal matrix for the heat equation in a matrix-free context
Following step-37 I tried to extend it for non-linear equations, using a
jacobian approximation for the "matrix" (which can be written as Jv = (F(u
+ ev) - F(u)) * 1/e, with e a small value (~1e-6), F() the residual of the
function, u the current solution, v the krylov vector obtained from the
GMRES-solver and J the jacobian matrix). Using this method I tried to
calculate the diagonal, which then is used for the preconditioner.
Based on the equation given above and the equation for the linear heat
equation, I get
F(u) = (u - u_0)/dt - (nabla^2u+nabla^2u_0)/2
F(u+ev) = F(u) + ev/dt - e(nabla^2v/2)
Thus Jv can be written as
Jv = (v/dt - nabla^2v/2)
or in weak form
Jv = v\phi/dt + 0.5*\nabla v\nabla\phi
Using that I tried to write the function for creating the diagonal
accordingly:
template
void JacobianOperator::local_compute_diagonal(
const MatrixFree & data,
vector_t ,
const unsigned int &,
const std::pair _range) const
{
FEEvaluation phi(data),
phi_old(data);
AlignedVector> diagonal(phi.dofs_per_cell);
for (unsigned int cell = cell_range.first; cell < cell_range.second; ++
cell)
{
phi_old.reinit(cell);
for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < phi.dofs_per_cell; ++j){
phi_old.submit_dof_value(VectorizedArray(), j);
}
phi_old.submit_dof_value(make_vectorized_array(1.), i);
phi_old.evaluate(true, true);
for (unsigned int q = 0; q < phi_old.n_q_points; ++q){
phi_old.submit_gradient(0.5
* (phi_old.get_gradient(q)),
q);
phi_old.submit_value((phi_old.get_value(q) / time_step),
q);
}
phi_old.integrate(true, true);
diagonal[i] = phi_old.get_dof_value(i);
}
for (unsigned int i = 0; i < phi_old.dofs_per_cell; ++i)
phi_old.submit_dof_value(diagonal[i], i);
phi_old.distribute_local_to_global(dst);
}
}
Still, the convergence behaviour is significantly worse than using an
Identity preconditioner. Thus I was wondering if this is related to the
wrong preconditioner, or due to wrong code? I am still trying to understand
how to create the preconditioner for that case most efficiently, thus there
might be some errors here.
Thanks!
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/cd41c84f-9441-4ca9-b48a-d8b9f120456f%40googlegroups.com.
[deal.II] Re: Thermoelastic Problem
Hi Daniel,
Thank you for the suggestion and I also think it can be one
of the way to successfully run both solutions by sharing information on
correct corresponding q_points.
I tried the same approach (as I am also using partitioned approach for
thermomechanical coupling) but so far my (mesh files being read for both
thermal and solid mechanics parts are same but ) triangulation are
different.
So I also want to use the same triangulation with different dof_handler for
both parts.
*The scenario of workflow in my case is following:*
The thermal code uses the *Triangulation triangulation_thermal* while
the solid mechanics code is using *parallel::distributed::Triangulation
triangulation_solid*. Yes the thermal part is programmed in serial and the
solid mechanics part in parallel. The solid mechanics code runs and makes
the grid with it's triangulation, As I want to use the same triangulation
so I try to copy this triangulation with function
*triangulation_solid.copy_triangulation(triangulation_thermal)
*and here the *triangulation_thermal *is the one which I will pass to the
name space and class of thermal analysis part before running it. Moreover
during my analysis of solid mechanics part I also refine the mesh (and
would also want to use the same refined mesh for the thermal analysis too
by copying it).
After describing the scenario following are my concerns where I need your
expert opinion :
*1)* when I use the *copy_triangulation() *function, I come across the
error (even here the mesh is not refined so far but still) mentioning :
* void dealii::Triangulation::copy_triangulation(const
dealii::Triangulation&) [with int dim = 2; int spacedim =
2]The violated condition was: (vertices.size() == 0) && (levels.size()
== 0) && (faces == nullptr)Additional information: You are trying to
perform an operation on a triangulation that is only allowed if the
triangulation is currently empty. However, it currently stores 100 vertices
and has cells on 1 levels.*
*2)* As far as the documentation of *tria.h* is concerned it also informs
that the *copy_triangulation()* cannot copy the refined mesh (where in my
case the mesh might also get refined during analysis based on certain
criteria and I would like to use same refined triangulation by copying it
for thermal analysis so that triangulation is same in both problems).
considering the above scenario as well as the concerns, I would be grateful
to receive any suggestion from your side. Hope I am clear in my description.
Waiting for your kind response. Thank you in advance!
*Best regards,*
Muhammad
On Thursday, March 3, 2016 at 10:49:37 AM UTC+1, Daniel Arndt wrote:
>
> Hello Sudharsan,
>
> what you want to do is to step through all the cells and reinitialize both
> DoFHandlers simultaneously, i.e. something like
>
> cell_elastic=dh_elastic.begin_active();
> cell_temperature=dh_temperature.begin_active();
> for (;cell_elastic != dh_elastic.end(); ++cell_elastic, ++cell_temperature)
> {
> fe_values_elastic.reinit(cell_elastic);
> fe_values_temperature.reinit(cell_temperature);
> fe_values_temperature.get_function_values(temperature,
> local_temperature);
> ...
> }
>
> This of course can only work if both DoFHandlers are based on the same
> Triangulation.
>
> Best regards,
> Daniel
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/1ca1a6d7-81a2-4353-af83-2d8e4d002c30%40googlegroups.com.
