Rotating Circle
The second example is a typical benchmark problem for the level set equation: a rotating bubble. The problem involves initializing (see Theory) with a "bubble" of radius 0.15 at for . This bubble is then advected with the given velocity field , so that, at , the bubble should return to its original position.
Level Set Equation
Figure 1 show the results of solving the rotating bubble problem with the level set equation alone, which initially behaves in a consistent manner. However, near the end of the simulation node-to-node oscillations appear in the solution, which is evident in the contour lines shown in Figure 1. As expected, these oscillations influence the area of the region encapsulated by the 0.5 level set contour as shown in Figure 4 and discussed in the Area Comparison section.
The complete input file for running this portion of the example is included below and it may be executed as follows.
cd ~/projects/moose/module/level_set/examples/rotating_circle
../../level_set-opt -i rotating_circle.i
[Mesh]
type = GeneratedMesh
dim = 2
xmin = -1
xmax = 1
ymin = -1
ymax = 1
nx = 32
ny = 32
uniform_refine = 2
[]
[AuxVariables]
[./vel_x]
[../]
[./vel_y]
[../]
[]
[AuxKernels]
[./vel_x]
type = FunctionAux
function = 4*y
variable = vel_x
execute_on = initial
[../]
[./vel_y]
type = FunctionAux
function = -4*x
variable = vel_y
execute_on = initial
[../]
[]
[Variables]
[./phi]
[../]
[]
[Functions]
[./phi_exact]
type = LevelSetOlssonBubble
epsilon = 0.03
center = '0 0.5 0'
radius = 0.15
[../]
[]
[ICs]
[./phi_ic]
type = FunctionIC
function = phi_exact
variable = phi
[../]
[]
[Kernels]
[./time]
type = TimeDerivative
variable = phi
[../]
[./advection]
type = LevelSetAdvection
velocity_x = vel_x
velocity_y = vel_y
variable = phi
[../]
[]
[Postprocessors]
[./area]
type = LevelSetVolume
threshold = 0.5
variable = phi
location = outside
execute_on = 'initial timestep_end'
[../]
[./cfl]
type = LevelSetCFLCondition
velocity_x = vel_x
velocity_y = vel_y
execute_on = 'initial' #timestep_end'
[../]
[]
[Executioner]
type = Transient
solve_type = PJFNK
start_time = 0
end_time = 1.570796
scheme = crank-nicolson
petsc_options_iname = '-pc_type -pc_sub_type'
petsc_options_value = 'asm ilu'
[./TimeStepper]
type = PostprocessorDT
postprocessor = cfl
scale = 0.8
[../]
[]
[Outputs]
csv = true
exodus = true
[]
Level Set Equation with SUPG
Adding SUPG stabilization—set the theory for details—mitigates the oscillations present in the first step, as shown in Figure 2. Adding the SUPG stabilization is trivial simply add the time and advection SUPG kernels to the input file (circle_rotate_supg.i) shown previously, the kernels block will then appear as:
[AuxKernels]
[./vel_x]
type = FunctionAux
function = 4*y
variable = vel_x
execute_on = initial
[../]
[./vel_y]
type = FunctionAux
function = -4*x
variable = vel_y
execute_on = initial
[../]
[]
Adding the stabilization improve the numerical solution but it suffers from a loss of conservation of the phase field variable, as discussed in the Area Comparison section and shown in Figure 4.
Level Set Equation with Reinitialization
Adding reinitializtion, in this case the scheme proposed by Olsson et al. (2007), requires the use of the MOOSE MultiApp. The enable reinitialization two input files are required: a master and sub-application.
The master input file must add the necessary MultiApps and Transfers blocks. For the problem at hand (circle_rotate_master.i) this easily accomplished by adding the following to the input file from the first step (i.e., do not include the SUPG kernels).
[MultiApps]
[./reinit]
type = LevelSetReinitializationMultiApp
input_files = 'circle_rotate_sub.i'
execute_on = 'timestep_end'
[../]
[]
[Transfers]
[./to_sub]
type = MultiAppCopyTransfer
source_variable = phi
variable = phi
direction = to_multiapp
multi_app = reinit
execute_on = 'timestep_end'
[../]
[./to_sub_init]
type = MultiAppCopyTransfer
source_variable = phi
variable = phi_0
direction = to_multiapp
multi_app = reinit
execute_on = 'timestep_end'
[../]
[./from_sub]
type = MultiAppCopyTransfer
source_variable = phi
variable = phi
direction = from_multiapp
multi_app = reinit
execute_on = 'timestep_end'
[../]
[]
Next, the sub-application input file must be created, which is shown below. This input file mimics the master input file closely, with three notable exceptions. First, the Kernels block utilize the time derivative and a new object, LevelSetOlssonReinitialization, that implements the reinitialization scheme of Olsson et al. (2007). Second, the Problem is set to use the LevelSetReinitializationProblem. Finally, the UserObjects block includes a terminator, LevelSetOlssonTerminator, which is responsible for stopping the reinitialization solve when steady-state is achieved according to the criteria defined by Olsson et al. (2007).
[Mesh]
type = GeneratedMesh
dim = 2
xmin = -1
xmax = 1
ymin = -1
ymax = 1
nx = 32
ny = 32
uniform_refine = 2
[]
[Variables]
[./phi]
[../]
[]
[AuxVariables]
[./phi_0]
[../]
[./marker]
[../]
[]
[Kernels]
[./time]
type = TimeDerivative
variable = phi
[../]
[./reinit]
type = LevelSetOlssonReinitialization
variable = phi
phi_0 = phi_0
epsilon = 0.03
[../]
[]
[Problem]
type = LevelSetReinitializationProblem
[]
[UserObjects]
[./arnold]
type = LevelSetOlssonTerminator
tol = 1
min_steps = 3
[../]
[]
[Executioner]
type = Transient
solve_type = NEWTON
start_time = 0
num_steps = 100
nl_abs_tol = 1e-14
scheme = crank-nicolson
line_search = none
petsc_options_iname = '-pc_type -pc_sub_type'
petsc_options_value = 'asm ilu'
dt = 0.003
[]
[Outputs]
[]
Figure 3 shows the results of the bubble problem with reinitialization, the result looks similar to the SUPG result. However, if you consider the area conservation discussed in the Area Comparison section, the reinitialization scheme yields the superior solution for this problem.
Figure 1: Results from solving the rotating circle problem with the level set equation alone.
Figure 2: Results from solving the rotating circle problem with the level set equation using SUPG stabilization.
Figure 3: Results from solving the rotating circle problem with the level set equation with reinitialization.
Area Comparison
Figure 4 is a plot of the area of the circle during the three simulations. Note that in the unstabilized, un-reinitialized level set equation, both area conservation and stability issues are readily apparent. The instabilities are especially obvious in Figure 4, where the drastic area changes are due to numerical oscillations in the solution field. Adding SUPG stabilization helps ameliorate the stability concern but it suffers from loss of area conservation. The re-initialization scheme is both stable and area-conserving.
Figure 4: Comparison of area inside the bubble during simulations.
The re-initialization methods performs well, but it is computationally expensive and picking the pseudo timestep size , steady-state criteria , and interface thickness correctly for the re-initialization problem is a non-trivial difficulty. Nevertheless, the level set module provides a valuable starting point and proof-of-concept for other researchers interested in the method, and the existing algorithm can no doubt be tuned to the needs of specific applications in terms of conservation and computational cost requirements.
References
- Elin Olsson, Gunilla Kreiss, and Sara Zahedi.
A conservative level set method for two phase flow ii.
Journal of Computational Physics, 225(1):785–807, 2007.
URL: http://dx.doi.org/10.1016/j.jcp.2006.12.027.[BibTeX]