date:20170918

Re: [deal.II] restrict a solution in some range?

2017-09-18 Thread Wolfgang Bangerth


On 09/18/2017 11:39 AM, Yiyang Zhang wrote:
I am following, e.g. Step-15, to solve a nonlinear system with Newton 
method.


The system is with the scalar potential V(phi)=(phi^2-1)^2.

The solution satisfies the boundary condition phi=1 at infinity, such 
that the total energy is finite.
Actually, we expect phi goes from 0 to 1 when the when the coordinate 
goes from center to infinity. phi should never exceeds 1.


That may be true for the *exact* solution of the PDE. But is it also 
true for the *discrete* solution?



However, when I do Newton iteration, phi always exceeds 1 near the 
boundary (my early description is wrong, it is not a local minima...), 
then things blow up...


Previously, I also solved such problems in 1D with finite difference. 
the field will also blow up at some point. In that case, after every 
iteration, I examine the values of the field at each discretized point. 
Wherever it exceed 1, I will manually replace it with 1. That will solve 
the problem.


Thus I wonder if I can do the same thing in dealii.


Yes, you can do that. You just have to loop over the elements of your 
solution vector and reset values that are too large or too small.


But you may lose the quadratic convergence of the Newton method this way 
(because the solution at the end of this step no longer satisfies the 
Newton update condition du=-J(u)^-1 F(u) ).


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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] restrict a solution in some range?

2017-09-18 Thread Yiyang Zhang

I am following, e.g. Step-15, to solve a nonlinear system with Newton 
method.

The system is with the scalar potential V(phi)=(phi^2-1)^2.

The solution satisfies the boundary condition phi=1 at infinity, such that 
the total energy is finite.
Actually, we expect phi goes from 0 to 1 when the when the coordinate goes 
from center to infinity. phi should never exceeds 1.

However, when I do Newton iteration, phi always exceeds 1 near the boundary 
(my early description is wrong, it is not a local minima...), then things 
blow up...

Previously, I also solved such problems in 1D with finite difference. the 
field will also blow up at some point. In that case, after every iteration, 
I examine the values of the field at each discretized point. Wherever it 
exceed 1, I will manually replace it with 1. That will solve the problem.

Thus I wonder if I can do the same thing in dealii. 

Thank you.

Best,
Yiyang

On Monday, September 18, 2017 at 11:15:02 AM UTC-5, Wolfgang Bangerth wrote:
>
> On 09/18/2017 09:33 AM, Yiyang Zhang wrote: 
> > 
> > For some problems, we might know in advance that the solution is in some 
> > range, e.g. [-1,1]. 
>
> You mean the exact solution u(x)? 
>
>
> > If we do not impose this restriction, the iteration may 
> > easily go to some other local minima. 
>
> Correct. But do you also know that the solution u_h(x) of the 
> *discretized* 
> equation satisfies the restriction that -1 <= u_h(x) <= 1 ? In general, 
> this 
> is not an easy thing to prove -- in fact, it will in general not be true. 
> And 
> if that statement is not true for the discrete solution u_h(x), then it is 
> not 
> correct to also *impose* this restriction. 
>
>
> > In dealii, how do we impose such a restriction on the solution vector 
> after 
> > each iteration? 
>
> You mean within the SolverCG, for example? Why would you want to impose it 
> in 
> each iteration, as opposed to just on the final solution of the linear 
> system? 
> Or do you mean in a nonlinear iteration? 
>
> Best 
>   W. 
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] restrict a solution in some range?

2017-09-18 Thread Wolfgang Bangerth


On 09/18/2017 09:33 AM, Yiyang Zhang wrote:


For some problems, we might know in advance that the solution is in some 
range, e.g. [-1,1].


You mean the exact solution u(x)?


If we do not impose this restriction, the iteration may 
easily go to some other local minima.


Correct. But do you also know that the solution u_h(x) of the *discretized* 
equation satisfies the restriction that -1 <= u_h(x) <= 1 ? In general, this 
is not an easy thing to prove -- in fact, it will in general not be true. And 
if that statement is not true for the discrete solution u_h(x), then it is not 
correct to also *impose* this restriction.



In dealii, how do we impose such a restriction on the solution vector after 
each iteration?


You mean within the SolverCG, for example? Why would you want to impose it in 
each iteration, as opposed to just on the final solution of the linear system? 
Or do you mean in a nonlinear iteration?


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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] restrict a solution in some range?

2017-09-18 Thread Yiyang Zhang

Hello!

For some problems, we might know in advance that the solution is in some 
range, e.g. [-1,1]. If we do not impose this restriction, the iteration may 
easily go to some other local minima. 

In dealii, how do we impose such a restriction on the solution vector after 
each iteration?

Thank you!

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: Refining a cylinder only in x

2017-09-18 Thread Bruno Turcksin

2017-09-18 10:32 GMT-04:00 Lucas Campos :

> Also, here is the new program, with the added if and
> prepare_coarsening_and_refinement(). You will find the interesting part
> between lines 85 and 104.
>
> You want to use prepare_coarsening_and_refinement after you have set the
flags not before.

On Monday, 18 September 2017 16:29:09 UTC+2, Lucas Campos wrote:
>>
>> Dear Bruno,
>>
>> You will find attached the resulting grids. While I originally expected
>> to find new divisions only along the x-direction, your previous answer
>> tells me it actually cuts in some local coordinate (whatever that might
>> be). In this case, what is the best way to refine the cylinder in a single
>> direction? I tried adding an if(cell->boundary_id() == 0) before setting
>> the flag, but it did not really help. The mesh is still being refined in
>> all directions.
>>
> Have you looked at step-30?

Best,

Bruno

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Refining a cylinder only in x

2017-09-18 Thread Lucas Campos

Also, here is the new program, with the added if and 
prepare_coarsening_and_refinement(). You will find the interesting part 
between lines 85 and 104.

On Monday, 18 September 2017 16:29:09 UTC+2, Lucas Campos wrote:
>
> Dear Bruno,
>
> You will find attached the resulting grids. While I originally expected to 
> find new divisions only along the x-direction, your previous answer tells 
> me it actually cuts in some local coordinate (whatever that might be). In 
> this case, what is the best way to refine the cylinder in a single 
> direction? I tried adding an if(cell->boundary_id() == 0) before setting 
> the flag, but it did not really help. The mesh is still being refined in 
> all directions.
>
> I ask this because I plan to do some transformations on this cylinder, to 
> turn it into a coil. In this case, the final result is much better if there 
> are enough cuts in the x-axis.
>
> Bests
>
> On Monday, 18 September 2017 14:44:40 UTC+2, Bruno Turcksin wrote:
>>
>> Lucas,
>>
>> On Monday, September 18, 2017 at 3:50:20 AM UTC-4, Lucas Campos wrote:
>>>
>>> I am trying to refine a cylinder in the x direction only, but it seems 
>>> that when I use 
>>>
>>> > cell->set_refine_flag(RefinementCase<3>::cut_x);
>>>
>>> the whole cylinder is refined. Indeed, if I run the above line for every 
>>> active cell *once* I would expect the number of cells to double. However, 
>>> they quintuple!
>>>
>> It looks good but how does the mesh look like? Here 
>> ,
>>  
>> you can see that cut_x is in the local coordinate system not the global 
>> one. So I would not be surprise if the mesh doesn't look like you think it 
>> should. Also if you do more than one refinement you should also use 
>> prepare_coarsening_and_refinement() before you refine your mesh.
>>
>> Best,
>>
>> Bruno
>>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 2013 - 2017 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Author: Timo Heister, Texas A University, 2013
 */

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 

#include 

using namespace dealii;

template 
void print_mesh_info(const Triangulation ,
 const std::string)
{
  std::cout << "Mesh info:" << std::endl
<< " dimension: " << dim << std::endl
<< " no. of cells: " << triangulation.n_active_cells() << std::endl;

  // Next loop over all faces of all cells and find how often each
  // boundary indicator is used (recall that if you access an element
  // of a std::map object that doesn't exist, it is implicitly created
  // and default initialized -- to zero, in the current case -- before
  // we then increment it):
  {
std::map boundary_count;
typename Triangulation::active_cell_iterator
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
  {
for (unsigned int face=0; faceface(face)->at_boundary())
  boundary_count[cell->face(face)->boundary_id()]++;
  }
  }

std::cout << " boundary indicators: ";
for (std::map::iterator it=boundary_count.begin();
 it!=boundary_count.end();
 ++it)
  {
std::cout << it->first << "(" << it->second << " times) ";
  }
std::cout << std::endl;
  }

  // Finally, produce a graphical representation of the mesh to an output
  // file:
  {
  std::ofstream out (filename + ".vtk");
  GridOut grid_out;
  grid_out.write_vtk (triangulation, out);
  }
  std::cout << " written to " << filename
<< std::endl
<< std::endl;
}

void grid_cylinder()
{
  Triangulation<3> triangulation;
  static const CylindricalManifold<3> boundary;

  GridGenerator::cylinder (triangulation, 0.05, 2);
  triangulation.set_manifold (0, boundary);

[deal.II] Re: Refining a cylinder only in x

2017-09-18 Thread Lucas Campos

Dear Bruno,

You will find attached the resulting grids. While I originally expected to 
find new divisions only along the x-direction, your previous answer tells 
me it actually cuts in some local coordinate (whatever that might be). In 
this case, what is the best way to refine the cylinder in a single 
direction? I tried adding an if(cell->boundary_id() == 0) before setting 
the flag, but it did not really help. The mesh is still being refined in 
all directions.

I ask this because I plan to do some transformations on this cylinder, to 
turn it into a coil. In this case, the final result is much better if there 
are enough cuts in the x-axis.

Bests

On Monday, 18 September 2017 14:44:40 UTC+2, Bruno Turcksin wrote:
>
> Lucas,
>
> On Monday, September 18, 2017 at 3:50:20 AM UTC-4, Lucas Campos wrote:
>>
>> I am trying to refine a cylinder in the x direction only, but it seems 
>> that when I use 
>>
>> > cell->set_refine_flag(RefinementCase<3>::cut_x);
>>
>> the whole cylinder is refined. Indeed, if I run the above line for every 
>> active cell *once* I would expect the number of cells to double. However, 
>> they quintuple!
>>
> It looks good but how does the mesh look like? Here 
> ,
>  
> you can see that cut_x is in the local coordinate system not the global 
> one. So I would not be surprise if the mesh doesn't look like you think it 
> should. Also if you do more than one refinement you should also use 
> prepare_coarsening_and_refinement() before you refine your mesh.
>
> Best,
>
> Bruno
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] how to find cell in which point lies

2017-09-18 Thread Wolfgang Bangerth


On 09/18/2017 06:38 AM, Bryukhanov Ilya wrote:


I solved the finite element task for the elasticity problem and obtained 
displacement fields. Then I did the same as has been done in the step-18
concerning accumulated stress data as mean stress for quadrature points. Now I 
need the function that returns the stress tensor for a given point.

How to find the cell in which given point lies?


You will want to use something like the VectorTools::point_gradient() function 
for what you want to do. You should take a look at its implementation to see 
how it works and how it finds a cell.


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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] how to find cell in which point lies

2017-09-18 Thread Bryukhanov Ilya

Hi,

I solved the finite element task for the elasticity problem and obtained 
displacement fields. Then I did the same as has been done in the step-18 
concerning accumulated stress data as mean stress for quadrature points. 
Now I need the function that returns the stress tensor for a given point.
How to find the cell in which given point lies?

Regards,
Bryukhanov Ilya

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: about deal.II with CUDA C programming acceleration

2017-09-18 Thread Bruno Turcksin

2017-09-17 1:22 GMT-04:00 Denis Davydov :

>
> I though you could also get binaries (nvcc) and the rest by compiling
> CUDA, i.e. https://developer.nvidia.com/compute/cuda/8.0/Prod2/
> local_installers/cuda_8.0.61_375.26_linux-run
> but of course one can also set it as an externally provided package.
>
> Sure, you can but the big advantage of using a package is when you update
your drivers or your kernel, you need to everything to always be
synchronized. When you use their package, they take care of that.

Bruno

-- 
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.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: about deal.II with CUDA C programming acceleration

2017-09-18 Thread Chih-Che Chueh

Hi Bruno,




> You definitely don't want to update your OS. I am also using ubuntu and I
> have several versions of gcc and clang installed. Lately, I have been using
> spack to install everything for me (https://github.com/LLNL/spack). This
> is becoming the standard way to install new programs on clusters at the DOE
> labs. Last time I installed gcc myself was about a year ago and I just
> followed the instruction here https://gcc.gnu.org/install/configure.html
> If your admin wants to shoot me an email with the errors he is getting, I
> will be happy to help.
>


Thanks for your kind help.

I have conveyed what you say above to our system administrator, who will
surely follow the instruction you provide. If he has any trouble with this
and need your help, he will tell me (as his English is not as good as mine)
and then I will let you know.

Sincerely,

Chih-Che

-- 
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.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Refining a cylinder only in x

2017-09-18 Thread Lucas Campos

Dear all,

I am trying to refine a cylinder in the x direction only, but it seems that 
when I use 

> cell->set_refine_flag(RefinementCase<3>::cut_x);

the whole cylinder is refined. Indeed, if I run the above line for every 
active cell *once* I would expect the number of cells to double. However, 
they quintuple!

Any ideas on how to fix this?

I am using an adapted version of step-49. The offending function is 
outlined below, and the full program is attached, together with its output.

void grid_cylinder()
>
> {
>
>   Triangulation<3> triangulation;
>
>   static const CylindricalManifold<3> boundary;
>
>
>>   GridGenerator::cylinder (triangulation, 0.05, 2);
>
>   triangulation.set_manifold (0, boundary);
>
>
>>   print_mesh_info (triangulation, "pre_grid");
>
>   for (int i=0; i<1; i++) {
>
>   for(auto cell: triangulation.active_cell_iterators()) {
>
>   cell->set_refine_flag(RefinementCase<3>::cut_x);
>
>   }
>
>   triangulation.execute_coarsening_and_refinement();
>
>   }
>
>
>>   print_mesh_info (triangulation, "grid");
>
> }
>
>

-- 
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.
For more options, visit https://groups.google.com/d/optout.
/* -
 *
 * Copyright (C) 2013 - 2017 by the deal.II authors
 *
 * This file is part of the deal.II library.
 *
 * The deal.II library is free software; you can use it, redistribute
 * it, and/or modify it under the terms of the GNU Lesser General
 * Public License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * The full text of the license can be found in the file LICENSE at
 * the top level of the deal.II distribution.
 *
 * -

 *
 * Author: Timo Heister, Texas A University, 2013
 */

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

#include 
#include 

#include 

using namespace dealii;

template 
void print_mesh_info(const Triangulation ,
 const std::string)
{
  std::cout << "Mesh info:" << std::endl
<< " dimension: " << dim << std::endl
<< " no. of cells: " << triangulation.n_active_cells() << std::endl;

  // Next loop over all faces of all cells and find how often each
  // boundary indicator is used (recall that if you access an element
  // of a std::map object that doesn't exist, it is implicitly created
  // and default initialized -- to zero, in the current case -- before
  // we then increment it):
  {
std::map boundary_count;
typename Triangulation::active_cell_iterator
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
  {
for (unsigned int face=0; faceface(face)->at_boundary())
  boundary_count[cell->face(face)->boundary_id()]++;
  }
  }

std::cout << " boundary indicators: ";
for (std::map::iterator it=boundary_count.begin();
 it!=boundary_count.end();
 ++it)
  {
std::cout << it->first << "(" << it->second << " times) ";
  }
std::cout << std::endl;
  }

  // Finally, produce a graphical representation of the mesh to an output
  // file:
  {
  std::ofstream out (filename + ".vtk");
  GridOut grid_out;
  grid_out.write_vtk (triangulation, out);
  }
  std::cout << " written to " << filename
<< std::endl
<< std::endl;
}

void grid_cylinder()
{
  Triangulation<3> triangulation;
  static const CylindricalManifold<3> boundary;

  GridGenerator::cylinder (triangulation, 0.05, 2);
  triangulation.set_manifold (0, boundary);

  print_mesh_info (triangulation, "pre_grid");
  for (int i=0; i<1; i++) {
  for(auto cell: triangulation.active_cell_iterators()) {
  cell->set_refine_flag(RefinementCase<3>::cut_x);
  }
  triangulation.execute_coarsening_and_refinement();
  }

  print_mesh_info (triangulation, "grid");
}

int main ()
{
  try
{
  grid_cylinder ();
}
  catch (std::exception )
{
  std::cerr << std::endl << std::endl
<< ""
<< std::endl;
  std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< ""
<< std::endl;

  return 1;
}
  catch (...)
{
  std::cerr << std::endl << std::endl

Re: [deal.II] restrict a solution in some range?

Re: [deal.II] restrict a solution in some range?

Re: [deal.II] restrict a solution in some range?

[deal.II] restrict a solution in some range?

Re: [deal.II] Re: Refining a cylinder only in x

[deal.II] Re: Refining a cylinder only in x

[deal.II] Re: Refining a cylinder only in x

Re: [deal.II] how to find cell in which point lies

[deal.II] how to find cell in which point lies

Re: [deal.II] Re: about deal.II with CUDA C programming acceleration

Re: [deal.II] Re: about deal.II with CUDA C programming acceleration

[deal.II] Refining a cylinder only in x

12 matches

Site Navigation

Mail list logo

Footer information