mathLab / pi-DoMUS

Parallel Deal.II MUltiphysics Solver
http://mathlab.github.io/pi-DoMUS/
GNU Lesser General Public License v2.1
14 stars 16 forks source link

'Maximum number of cells' parameter unused. #216

Open agzimmerman opened 7 years ago

agzimmerman commented 7 years ago

At least for the Poisson problem interface (I haven't tested others), the 'Maximum number of cells' parameter in the Refinement subsection appears to be unused. I mentioned this to @asartori86. who said he also experienced some issues with this parameter.

I already have made what I would call a "work around" rather than an ideal solution. I will request a pull, since I think that's clearer than explaining more here.

It's entirely possible that I'm just using pi-DoMUS incorrectly. Does this parameter work for anyone? Is there already a test that uses it successfully?

asartori86 commented 7 years ago

the 'Maximum number of cells' parameter in the Refinement subsection appears to be unused.

Have you changed the refinement strategy from fraction to number?

agzimmerman commented 7 years ago

I tried this now and the program enters what appears to be an infinite loop once max cells (here 10,000) is exceeded:

  ################ restart ######### 
max_kelly > threshold
0.00253967 >  0.0001
######################################
Number of active cells: 10.208 (on 7 levels)
Number of degrees of freedom: 10.832(10832)

  ################ restart ######### 
max_kelly > threshold
0.00253967 >  0.0001
######################################
Number of active cells: 10.208 (on 7 levels)
Number of degrees of freedom: 10.832(10832)

  ################ restart ######### 
max_kelly > threshold
0.00253967 >  0.0001
######################################
Number of active cells: 10.208 (on 7 levels)
Number of degrees of freedom: 10.832(10832)

(this keeps going until I stop the program manually)

At the end of this message I'll paste used_parameters.prm. Does anything look wrong? I'm new to the community and I don't see a way to attach a file (except for hosting it somewhere and linking to it); so please do let me know the appropriate way to share these things :)

I haven't yet looked in to what's causing the loop, and it might be a while before I can take time to dig into this.

As a side note: It seems odd to me that the "refinement strategy" would change whether or not the max_cells parameter is used. I would think the parameter should work for both the fraction and number strategies. Maybe I'm missing something. Is it just not implemented yet, or am I misunderstanding these strategies?

Thanks for your help.

used_parameters.prm:

# Parameter file generated with 
# D2K_GIT_BRANCH=       master
# D2K_GIT_SHORTREV=     6b09041
# DEAL_II_GIT_BRANCH=   
# DEAL_II_GIT_SHORTREV= e3ff743
subsection Dirichlet boundary conditions
  set IDs and component masks = 0=ALL % 1=ALL
  set IDs and expressions     = 0=cold % 1=hot
  set Known component names   = u
  set Used constants          = cold=-1.0, hot=1.0
end
subsection Domain
  set Colorize                      = false
  set Copy boundary to manifold ids = true
  set Copy material to manifold ids = false
  set Create default manifolds      = false
  set Grid to generate              = hyper_cube_with_cylindrical_hole
  set Input grid file name          = 
  set Manifold descriptors          = 1=HyperBallBoundary
  set Mesh smoothing alogrithm      = none
  set Optional Point<spacedim> 1    = 0,0
  set Optional Point<spacedim> 2    = 1,1
  set Optional double 1             = 0.25
  set Optional double 2             = 0.5
  set Optional double 3             = 1
  set Optional int 1                = 1
  set Optional int 2                = 2
  set Optional vector of dim int    = 1,1
  set Output grid file name         = 
end
subsection Error Tables
  set Compute error            = false
  set Error file format        = tex
  set Error precision          = 3
  set Output error tables      = true
  set Solution names           = u
  set Solution names for latex = u
  set Table names              = error
  set Write error files        = false
  subsection Table 0
    set Add convergence rates          = true
    set Extra terms                    = cells,dofs
    set Latex table caption            = error
    set List of error norms to compute = L2,H1
    set Rate key                       = 
  end
end
subsection Exact solution
  set Function constants  = 
  set Function expression = 0
  set Variable names      = x,y,t
end
subsection Forcing terms
  set IDs and component masks = 
  set IDs and expressions     = 
  set Known component names   = u
  set Used constants          = 
end
subsection IDA Solver Parameters
  set Absolute error tolerance                      = 1e-4
  set Final time                                    = 1
  set Ignore algebraic terms for error computations = false
  set Initial condition Newton max iterations       = 5
  set Initial condition Newton parameter            = 0.33
  set Initial condition type                        = use_y_diff
  set Initial condition type after restart          = use_y_dot
  set Initial step size                             = 1e-4
  set Initial time                                  = 0
  set Maximum number of nonlinear iterations        = 10
  set Maximum order of BDF                          = 5
  set Min step size                                 = 5e-5
  set Relative error tolerance                      = 1e-3
  set Seconds between each output                   = 1e-1
  set Show output of time steps                     = true
  set Use local tolerances                          = false
end
subsection IMEX Parameters
  set Absolute error tolerance                     = 0.000001
  set Final time                                   = 1.000000
  set Initial time                                 = 0.000000
  set Intervals between outputs                    = 1
  set Maximum number of inner nonlinear iterations = 3
  set Maximum number of outer nonlinear iterations = 5
  set Method used                                  = fixed_alpha
  set Newton relaxation parameter                  = 1.000000
  set Number of elements in backtracking sequence  = 5
  set Print useful informations                    = false
  set Relative error tolerance                     = 0.000000
  set Step size                                    = 1e-2
  set Update continuously Jacobian                 = true
  set Use the KINSOL solver                        = true
end
subsection Initial solution
  set Function constants  = cold=-1.0
  set Function expression = cold
  set Variable names      = x,y,t
end
subsection Initial solution_dot
  set Function constants  = 
  set Function expression = 0
  set Variable names      = x,y,t
end
subsection KINSOL for IMEX
  set Level of verbosity of the KINSOL solver            = 0
  set Maximum number of iteration before Jacobian update = 10
  set Maximum number of iterations                       = 200
  set Step tolerance                                     = 1e-11
  set Strategy                                           = newton
  set Tolerance for residuals                            = 1e-9
  set Use internal KINSOL direct solver                  = false
end
subsection Neumann boundary conditions
  set IDs and component masks = 
  set IDs and expressions     = 
  set Known component names   = u
  set Used constants          = 
end
subsection Output Parameters
  set Files to save in run directory = 
  set Incremental run prefix         = 
  set Output format                  = vtu
  set Output partitioning            = false
  set Problem base name              = solution
  set Solution names                 = u
  set Subdivisions                   = 1
end
subsection Poisson problem
  set Block of differential components = 1
  set Blocking of the finite element   = u
  set Finite element space             = FESystem[FE_Q(1)]
end
subsection Refinement
  set Bottom fraction                        = 0.100000
  set Maximum number of cells (if available) = 10000
  set Order (optimize)                       = 2
  set Refinement strategy                    = number
  set Top fraction                           = 0.300000
end
subsection Time derivative of Dirichlet boundary conditions
  set IDs and component masks = 
  set IDs and expressions     = 
  set Known component names   = u
  set Used constants          = 
end
subsection Zero average constraints
  set Known component names        = u
  set Zero average on boundary     = 
  set Zero average on whole domain = 
end
subsection pi-DoMUS
  set Adaptive refinement                            = true
  set Enable finer preconditioner                    = false
  set Initial global refinement                      = 1
  set Jacobian solver tolerance                      = 1e-8
  set Max iterations                                 = 50
  set Max iterations finer prec.                     = 0
  set Max tmp vectors                                = 30
  set Max tmp vectors for finer system               = 50
  set Number of cycles                               = 1
  set Overwrite Newton's iterations                  = true
  set Print some useful informations about processes = true
  set Refine mesh during transient                   = true
  set Threshold for solver's restart                 = 1e-4
  set Time stepper                                   = euler
  set Use direct solver if available                 = true
end
asartori86 commented 7 years ago

@alexanderzimmerman, I've updated your comment with proper style. When you need to include code snippets or text files use ``` before the first line of code and after the last one.

RE the infinite loop, the point is that you have set a threshold for refinement that is quite small

set Threshold for solver's restart                 = 1e-4

and set a maximum number of cells that is not-so-big

set Maximum number of cells (if available) = 10000

Our solvers try to refine the mesh until the infinity norm of the error estimation is below the given threshold. So you can play with three actions in order to avoid the infinite loop

  1. decrease the threshold
  2. increase the maximum number of cells
  3. change how the error is estimated
agzimmerman commented 7 years ago

@asartori86, I don't understand your answer. What is the "maximum number of cells" parameter supposed to do? I have been moving forward with my "work-around" for the past couple of weeks, but currently it's my only difference from the master branch, so I'd like to figure this out. Thanks for your help :)

Edit: To clarify, I essentially want to specify some maximum number of cells and some accuracy, and I want the code to adaptively refine the grid until either A. it exceeds the maximum number of cells or B. it meets the accuracy requirement. I think this is a typical capability that people would generally want. I think I could dig up a deal.II tutorial where they do this (unless I'm remembering it wrong).

agzimmerman commented 7 years ago

Update: I see the same behavior with the deal2lkit heat equation example. Is there a way for me to reference this pi-DoMUS issue directly if I raise the issue at the deal2lkit repo?

luca-heltai commented 7 years ago

@alexanderzimmerman , pi-DoMUS uses deal2lkit for this, so yes, this is how it goes. You can refer this issue by its html address in deal2lkit issues, and github will act automatically: https://github.com/mathLab/pi-DoMUS/issues/216 (this was obtained by copy-pasting the address).

The deal2lkit interface is a wrapper to the functionality in deal.II. What you are asking is not implemented in deal.II, at the moment.

The "maximum number of cells" is supposed to do exactly what it says, but not all strategies are supported in deal.II for this.

Ideally what you want to achieve coincides with the goals of this class. As you found out, however, it is not trivial to put into practice what you suggest. For example, if you exceed the number of cells, what should happen next step? Should one increase the accuracy requirements? Should this be done automatically? If you don't increase the accuracy requirements, the refinement strategy will always propose something that breaks the first rule (# of cells). This is something that we should add, but so far we had no stringent need for this.

If you come up with an idea to obtain this in a robust manner, please do provide a PR and we'll merge it in deal2lkit (pi-domus will pick it up automatically).

luca-heltai commented 7 years ago

Ps: @asartori86 , we should document this much better in deal2lkit. I'll open an issue.