  We are concerned with the simulation and control of a two phase flow model
governed by a coupled Cahn-Hilliard Navier-Stokes system involving a nonsmooth
energy potential. We establish the existence of optimal solutions and present
two distinct approaches to derive suitable stationarity conditions for the
bilevel problem, namely C- and strong stationarity. Moreover, we demonstrate
the numerical realization of these concepts at the hands of two adaptive
solution algorithms relying on a specifically developed goal-oriented error
estimator. In addition, we present a model order reduction approach using
proper orthogonal decomposition (POD-MOR) in order to replace high-fidelity
models by low order surrogates. In particular, we combine POD with
space-adapted snapshots and address the challenges which are the consideration
of snapshots with different spatial resolutions and the conservation of a
solenoidal property.
