kronbichler/adaflo

adaflo, an adaptive finite element solver for incompressible fluid flow and two-phase flow

adaflo's Issues

Introduce namespace `adaflo`

In a student project, we use adaflo and the deal.II TimeStepping class (https://www.dealii.org/current/doxygen/deal.II/classTimeStepping_1_1TimeStepping.html). Since adaflo also has a TimeStepping class (see https://github.com/kronbichler/adaflo/blob/master/include/adaflo/time_stepping.h), we get a naming collision. @kronbichler would you agree, if we introduced an adaflo namespace around the whole code? It's a big incompatibility change, but it would help avoid naming collisions in the future.

@j-brotz FYI

1D: Nonlinear solver for the Navier-Stokes equations does not converge anymore

Unfortunately, with the current deal.II master (dealii/dealii@f6081ed) we observe convergence issues of the nonlinear solver in 1D. This issue is a follow-up to #84.

Using dealii/dealii@7e1098b, the 1D_flow example of adaflo converges as expected

Running a 1D flow using BDF-2, Q2/Q1 elements
Number of active cells: 2048.
Number of degrees of freedom (velocity/pressure): 6146 (4097 + 2049).
Approximate size last cell: 0.0012207

Time step #1, advancing from t_n-1 = 0 to t = 0.01 (dt = 0.01). 

  Nonlin Res     Prec Upd     Increment   Lin Iter     Lin Res
  ____________________________________________________________
  2.236e+00        ILUs       6.91e+01       29       6.94e-06
  6.940e-06        ---        7.16e-04       25       4.61e-10
  4.602e-10     converged.


Time step #2, advancing from t_n-1 = 0.01 to t = 0.02 (dt = 0.01). 

  Nonlin Res     Prec Upd     Increment   Lin Iter     Lin Res
  ____________________________________________________________
  2.470e-02        ILUs       2.56e-01        1       1.07e-09
  1.024e-09        ---        1.60e-08        1       3.78e-10
  3.779e-10     converged.


Time step #3, advancing from t_n-1 = 0.02 to t = 0.03 (dt = 0.01). 

  Nonlin Res     Prec Upd     Increment   Lin Iter     Lin Res
  ____________________________________________________________
  1.297e-09        ---        4.17e-08        1       4.15e-10
  4.149e-10     converged.


Time step #4, advancing from t_n-1 = 0.03 to t = 0.04 (dt = 0.01). 

  Nonlin Res     Prec Upd     Increment   Lin Iter     Lin Res
  ____________________________________________________________
  8.118e-10        ---        1.35e-07        3       3.16e-10
  3.165e-10     converged.


Time step #5, advancing from t_n-1 = 0.04 to t = 0.05 (dt = 0.01). 

  Nonlin Res     Prec Upd     Increment   Lin Iter     Lin Res
  ____________________________________________________________
  1.655e-09        ILUs       1.24e-07        2       4.69e-10
  4.687e-10     converged.

with the current master we obtain

Running a 1D flow using BDF-2, Q2/Q1 elements
 Number of active cells: 2048.
 Number of degrees of freedom (velocity/pressure): 6146 (4097 + 2049).
 Approximate size last cell: 0.0012207

Time step #1, advancing from t_n-1 = 0 to t = 0.01 (dt = 0.01). 

   Nonlin Res     Prec Upd     Increment   Lin Iter     Lin Res
   ____________________________________________________________
   2.236e+00        ILUs       4.26e+02       30       1.39e+00
   1.389e+00        ---        1.06e+02       30       1.38e+00
   1.382e+00        ---        3.81e+01       30       1.38e+00
   1.382e+00        ---        1.80e+01       30       1.38e+00
   1.381e+00        ---        9.05e+00       30       1.38e+00
   1.381e+00        ---        4.67e+00       30       1.38e+00
   1.381e+00        ---        2.48e+00       30       1.38e+00
   1.381e+00        ---        1.32e+00       30       1.38e+00
   1.381e+00        ---        7.15e-01       30       1.38e+00
   1.381e+00        ---        3.85e-01       30       1.38e+00
   1.381e+00    Warning: nonlinear iteration did not converge!

Time step #2, advancing from t_n-1 = 0.01 to t = 0.02 (dt = 0.01). 

   Nonlin Res     Prec Upd     Increment   Lin Iter     Lin Res
   ____________________________________________________________
   1.381e+00        ILUs       2.74e+02       30       1.35e+00
   1.349e+00        ---        1.31e+01       30       1.35e+00
   1.349e+00        ---        2.47e+00       30       1.35e+00
   1.349e+00        ---        1.59e-01       30       1.35e+00
   1.349e+00        ---        4.53e-02       30       1.35e+00
   1.349e+00        ---        5.61e-03       30       1.35e+00
   1.349e+00        ---        1.31e-03       30       1.35e+00
   1.349e+00        ---        2.46e-04       30       1.35e+00
   1.349e+00        ---        5.28e-05       30       1.35e+00
   1.349e+00        ---        1.24e-05       30       1.35e+00
   1.349e+00    Warning: nonlinear iteration did not converge!

Time step #3, advancing from t_n-1 = 0.02 to t = 0.03 (dt = 0.01). 

   Nonlin Res     Prec Upd     Increment   Lin Iter     Lin Res
   ____________________________________________________________
   1.381e+00        ---        2.60e+02       30       1.35e+00
   1.350e+00        ---        9.67e+00       30       1.35e+00
   1.350e+00        ---        2.73e+00       30       1.35e+00
   1.350e+00        ---        4.47e-01       30       1.35e+00
   1.350e+00        ---        1.83e-01       30       1.35e+00
   1.350e+00        ---        4.46e-02       30       1.35e+00
   1.350e+00        ---        1.70e-02       30       1.35e+00
   1.350e+00        ---        4.65e-03       30       1.35e+00
   1.350e+00        ---        1.70e-03       30       1.35e+00
   1.350e+00        ---        4.85e-04       30       1.35e+00
   1.350e+00    Warning: nonlinear iteration did not converge!

Time step #4, advancing from t_n-1 = 0.03 to t = 0.04 (dt = 0.01). 

   Nonlin Res     Prec Upd     Increment   Lin Iter     Lin Res
   ____________________________________________________________
   1.356e+00        ---        4.45e+00       30       1.36e+00
   1.356e+00        ---        1.08e+00       30       1.36e+00
   1.356e+00        ---        2.54e-01       30       1.36e+00
   1.356e+00        ---        1.81e-01       30       1.36e+00
   1.356e+00        ---        5.37e-02       30       1.36e+00
   1.356e+00        ---        3.61e-02       30       1.36e+00
   1.356e+00        ---        1.19e-02       30       1.36e+00
   1.356e+00        ---        7.55e-03       30       1.36e+00
   1.356e+00        ---        2.67e-03       30       1.36e+00
   1.356e+00        ---        1.62e-03       30       1.36e+00
   1.356e+00    Warning: nonlinear iteration did not converge!

Time step #5, advancing from t_n-1 = 0.04 to t = 0.05 (dt = 0.01). 

   Nonlin Res     Prec Upd     Increment   Lin Iter     Lin Res
   ____________________________________________________________
   1.366e+00        ---        4.24e+00       30       1.37e+00
   1.366e+00        ---        2.26e+00       30       1.37e+00
   1.366e+00        ---        1.11e+00       30       1.37e+00
   1.366e+00        ---        1.04e+00       30       1.37e+00
   1.366e+00        ---        6.06e-01       30       1.37e+00
   1.366e+00        ---        5.64e-01       30       1.37e+00
   1.366e+00        ---        3.54e-01       30       1.37e+00
   1.366e+00        ---        3.25e-01       30       1.37e+00
   1.366e+00        ---        2.14e-01       30       1.37e+00
   1.366e+00        ---        1.94e-01       30       1.37e+00
   1.366e+00    Warning: nonlinear iteration did not converge!

@kronbichler Before I start investigating the issues further, do you have any idea what change in dealii/dealii#16776 might be causing the problem? At first glance, I did not see anything obvious.

See over test and their run time

As noticed in #3 and detailed in #3 (comment), we look over the tests and make them run more quickly. Also, the output of iterations is brittle, so might need some work.

Update version of clang-format to v.11 to be consistent with deal.II master

@peterrum @kronbichler I suggest to update the CI's version of clang-format to v.11 for the indentation check to be consistent with current deal.II master. If you agree, I'll resolve this issue in a PR.

question about phase field solver

Hello Krobichler,
Is there any reference paper about phase field two-phase flow solver ?
In the todo.txt, you said that "Our simulations often crash
because the curvature gets very bad, which is in turn caused by the interface
profile that gets some small oscillations".
Is the phase field solver Robust enough to handle the big curvature problem?

Add parameter `output memory`

Compilation errors with `dim==1`

Currently, adaflo does not compile for dim==1 with current deal.II master due to the changes made in dealii/dealii#16770

• One issue is related to a call of FEEvaluation::get_divergence() for dim==1 here

  adaflo/source/navier_stokes_matrix.cc

  Line 954 in 6329eee

  pressure.submit_value(weight * velocity.get_divergence(q), q);

  which is disabled by dealii/dealii#16770.
• Another issue is related to the fact that we introduced utility functions for dim==1, n_components==dim which converts the
  gradient_type = Tensor<1,dim,VectorizedArray<double>> obtained from the FEEvaluation access to Tensor<2,dim,VectorizedArray<double>>, to obtain the same types as for dim>1 via

  adaflo/source/navier_stokes_matrix.cc

  Lines 703 to 704 in 6329eee

  	Tensor<2, dim, vector_t> grad_u =
  	convert_to_tensor<2, dim, vector_t>(velocity.get_gradient(q));

  and

  adaflo/include/adaflo/util.h

  Lines 146 to 179 in 6329eee

  	/**
  	* If dim == 1, convert a tensor of rank-1 to a tensor of rank. This function
  	* is useful to obtain equal, tensor-rank-valued return types of
  	* FEEvaluation-operations for dim == 1 and dim > 1.
  	*/
  	template <int rank_, int dim, typename VectorizedArrayType = VectorizedArray<double>>
  	static Tensor<rank_, dim, VectorizedArrayType>
  	convert_to_tensor(const Tensor<rank_ - 1, dim, VectorizedArrayType> &in)
  	{
  	Assert(dim == 1, ExcMessage("This operation is not permitted for dim>1."));
  	if (rank_ == 2)
  	{
  	Tensor<2, dim, VectorizedArrayType> tens;
  	tens[0][0] = in[0];
  	return tens;
  	}
  	else
  	{
  	Assert(false, ExcNotImplemented());
  	}
  	}
  	/**
  	* This function overloads the previous convert_to_tensor function, when the input
  	* tensor already has the desired tensor rank.
  	*/
  	template <int rank_, int dim, typename VectorizedArrayType = VectorizedArray<double>>
  	static Tensor<rank_, dim, VectorizedArrayType>
  	convert_to_tensor(const Tensor<rank_, dim, VectorizedArrayType> &in)
  	{
  	return in;
  	}

  . Our motivation for this was to straightforward loop over the components and make component-wise modifications. However, dealii/dealii#16770 requires to submit the original gradient_type, which makes e.g.

  adaflo/source/navier_stokes_matrix.cc

  Line 891 in 6329eee

  velocity.submit_gradient(grad_u, q);

  invalid.

My workaround for the described cases looks currently as follows

adaflo/source/navier_stokes_matrix.cc

Lines 982 to 986 in 2acdbf6

	if constexpr (dim == 1)
	pressure.submit_value(weight * velocity.get_gradient(q), q);
	else
	pressure.submit_value(weight * velocity.get_divergence(q), q);
	}

adaflo/source/navier_stokes_matrix.cc

Lines 916 to 919 in 2acdbf6

	if constexpr (dim == 1)
	velocity.submit_gradient(grad_u[0], q);
	else
	velocity.submit_gradient(grad_u, q);

@kronbichler do you have any suggestions how to implement such special cases with dim==1 in the future?
I want to make sure that we are on the same page here. Just as a remark: dealii/dealii#16776 requires even more changes.

Check the ghost state of vectors

While going through the code, I see many places where no ghost values updated although I would have expected at these places.

Let's introduce inside of MatrixFree::loop() asserts that values are not ghosted and find all relevant places.

Override warning [-Wsuggest-override]

I get a number of warnings concerning missing overrides with GCC 12.2.1. I'll fix them asap.