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)) print('{:.2e}'.format(
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 file
[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 &dst, const unsigned int &, const std::pair &cell_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.