Anisotropic parabolic problems with slowly or~rapidly growing terms
Agnieszka \'Swierczewska-Gwiazda
{Banacha 2, 02-097 Warsaw, Poland}
We consider an abstract parabolic problem in a framework of maximal monotone graphs, possibly multi-valued with growth conditions formulated with help of an $x-$dependent $N-$function. The main novelty of the paper consists in the lack of any growth restrictions on the ${ N}$--function combined with its anisotropic character, namely we allow the dependence on all the directions of the gradient, not only on its absolute value. This leads us to use the notion of modular convergence and studying in detail the question of density of compactly supported smooth functions with respect to the modular convergence. \medskip \noindent {AMS 2000 Classification:} 35K55, 35K20 \noindent {Keywords}: Musielak -- Orlicz spaces, modular convergence, nonlinear parabolic inclusion, maximal monotone graph % As a separate fact, the case of lower regularity with respect to $x$ is also considered under additional assumption on the growth of an ${\cal N}$--function (the so-called $\Delta_2-$condition for a conjugate function). %In the present paper we study the existence of weak solutions to an abstract parabolic initial-boundary value %problem. On the operator appearing in the equation we assume the coercivity conditions given by an %${\cal N}$--function (i.e. convex function satisfying conditions specified in the paper). The main novelty of the paper consists in the lack of any growth restrictions on the ${\cal N}$--function combined with an anisotropic character of the ${\cal N}$--function, namely we allow the dependence on all the directions of the gradient, not only on its absolute value. \\[2ex]
swierczewska_revised.tex
\documentclass[preprint,12pt]{article}
%\documentclass[12pt]{article}
\usepackage{amssymb}
\usepackage{amsfonts}
%\usepackage[notref,notcite]{showkeys}
\usepackage{color}
\usepackage{graphics}
\usepackage{epsfig,psfrag,graphicx}
%%
\usepackage{amssymb}
\usepackage{eepic,epic}

%\usepackage{natbib} % more flexible citations
\usepackage{epsfig} % postscript figures
\usepackage{amsmath} % ams math
\usepackage{amssymb} % ams symbols
\usepackage{amsbsy} % ams bold symbols
%%
\newcommand{\cblack}{\color{black}}
\newcommand{\cred}{\color{red}}
%\textheight22cm \textwidth15cm \voffset=-1.5cm \hoffset=-1cm
%\date{}
%%%%%%%%%%%%%%%%%%%%%%%%%
\newtheorem{definition}{Definition}[section]
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
%%%
%%%%%%%%%%%%%%
\renewcommand{\theequation}{\thesection.\arabic{equation}}
\newcommand{\nc}{\newcommand}
\newcommand{\mA}{\mathcal{A}}
\nc{\V}{{\cal V}} \nc{\M}{{\cal M}} \nc{\T}{{\cal T}}
\def\u{u}
\nc{\D}{{\cal D}} \nc{\W}{\mathbb W} \nc{\ti}{\tilde}
\nc{\wti}{\widetilde} \nc{\vep}{\varepsilon}
%\nc{\vp}{\varphi}
\nc{\R}{{\mathbb R}} \nc{\N}{{\mathbb N}}
%\renewcommand{\b}[1]{\mathbf{#1}}
\renewcommand{\S}{{\mathbb S}}
\nc{\di}{\displaystyle}
\renewcommand{\b}{\mathbf}
\nc{\pa}{\partial} \nc{\lra}{\longrightarrow}
% \nc{\qed}{\mbox{
%}\nolinebreak\hfill \rule{2mm}{2mm}}
 \nc{\weak}{\rightharpoonup}
\nc{\weakstar}{\stackrel{\ast}{\rightharpoonup}} \nc{\proof}{{\bf
Proof: }} \nc{\en}{^{\varepsilon,n}} \nc{\Rd}{{{\mathbb R}^{d}}}
\nc{\biting}{\stackrel{\,\,b}{\rightarrow}}
\newcommand{\ud}{\,\mathrm d}
\renewcommand{\div}{{\mathrm{div}}\,}
\nc{\bB}{B}
\nc{\bS}{A}
\newcommand{\vrho}{\varrho}
\nc{\modular}[1]{{\stackrel{ #1}{\longrightarrow\,}}}
\newcommand{\Ref}[1]{{\rm(\ref{#1})}}
\nc{\vd}{\bar{v}} \nc{\zd}{\bar{z}}
\def\bbbone{{\mathchoice {\rm 1\mskip-4mu l}
{\rm 1\mskip-4mu l} {\rm 1\mskip-4.5mu l} {\rm 1\mskip-5mu l}}}
\title{Anisotropic  parabolic problems with slowly or~rapidly growing terms  }
\author{Agnieszka \'Swierczewska-Gwiazda\\[4ex]
{\small Institute of Applied Mathematics and Mechanics}, \\{\small University of Warsaw},\\
{\small Banacha 2, 02-097 Warsaw, Poland}\\ {\small \texttt{aswiercz@mimuw.edu.pl}}}
%}
%\address{Institute of Applied Mathematics and Mechanics, \\University of Warsaw,\\Banacha 2, 02-097 Warsaw, Poland}
\begin{document}
\date{}
%\begin{frontmatter}
\maketitle
\begin{abstract}
We consider an abstract parabolic problem in a framework of maximal monotone graphs, possibly multi-valued with growth conditions formulated with help of an $x-$dependent $N-$function.
The main novelty of the paper consists in the lack of any growth restrictions on the ${ N}$--function combined with its anisotropic character, namely we allow the dependence on all the directions of the gradient, not only on its absolute value. This leads us to use the notion of modular convergence and studying in detail the question of density of compactly supported smooth functions with respect to the modular convergence. 

\medskip

\noindent
{\bf AMS 2000 Classification:} 35K55, 35K20 

\noindent
{\bf Keywords}: Musielak -- Orlicz spaces, modular
convergence, nonlinear parabolic inclusion, maximal monotone graph
% As a separate fact, the case of lower regularity with respect to $x$ is also considered under additional  assumption on the growth of an ${\cal N}$--function (the so-called $\Delta_2-$condition for a conjugate function).
%In the present paper we study the existence of weak solutions to an abstract parabolic initial-boundary value
%problem. On the operator appearing in the equation we assume the coercivity conditions given by an 
%${\cal N}$--function (i.e. convex function satisfying conditions specified in the paper). The main novelty of the paper consists in the lack of any growth restrictions on the ${\cal N}$--function combined with an anisotropic character of the ${\cal N}$--function, namely we allow the dependence on all the directions of the gradient, not only on its absolute value. \\[2ex]
\end{abstract}
%\begin{keyword}
%{\bf AMS 2000 Classification:} 35K55, 35K20 
%Musielak -- Orlicz spaces, modular
%convergence, nonlinear parabolic inclusion, maximal monotone graph
%\end{keyword}
%\end{frontmatter}

\section{Introduction}
 
Our interest is directed to the phenomenon of anisotropic behaviour in a parabolic problem. 
The proposed approach allows for capturing very general form of growth conditions of a nonlinear term.  
We concentrate on an abstract parabolic problem.
Let $\Omega\subset\R^d$ be an open, bounded set with a ${\cal C}^2$~boun\-da\-ry~$\partial \Omega$,  $(0,T)$ be the time interval with $T<\infty$, $Q:=(0,T)\times\Omega$ and 
${\mathcal A}$ be a maximal monotone graph satisfying the assumptions (A1)--(A5) formulated below. 
Given $f$  and $u_0$ we want to find 
 $u:Q\to\R$ and $ A:Q\to\R^d$ such that 
\begin{align}\label{P1a}
u_t-\div A = f\quad&{\rm in}\ Q,\\
\label{P1aa}
(\nabla u,A)\in{\mathcal A}(t,x)\quad&{\rm in}\ Q,\\
\label{P2a}
u(0,x)=u_0\quad&{\rm in}\ \Omega,\\
\label{P3a}
u(t,x)=0\quad&{\rm on}\ (0,T)\times\partial\Omega.
\end{align}

The main objective of the present paper is to obtain existence result for the widest possible class of maximal monotone graphs. Hence various non-standard possibilities are considered including anisotropic growth conditions, $x-$dependent growth conditions and also relations given by maximal monotone graph. The last ones provide the possibility of generalization of discontinuous relations, namely considering $A$ as a discontinuous function of $\nabla u$, where the jumps of $A$  are filled by intervals creating vertical parts of the graph $\mA$.
Most of these generalities shall arise in a function that will prescribe the growth/coercivity conditions. Contrarty to the usual case of Leray-Lions type operators, where the polynomial growth is assumed,
e.g. $|A(\xi)|\le c(1+|\xi|)^{p-1},\,A(\xi)\cdot\xi\ge C|\xi|^p$ for some nonnegative constants $c, C$ and $p>1$ we shall 
work with $N-$functions.  By an $N-$function we mean that  $M:\bar \Omega\times\R^d\to\R_+$, $M(x,a)$
  is measurable w.r.t. $x$ for all $a\in\R^d$ and  continuous w.r.t. $a$ for a.a. $x\in\bar \Omega$,  convex in $a$, has superlinear growth, $M(x,a)=0$ iff $a=0$ and  
$$\lim_{|a|\to\infty}\inf_{x\in \Omega}\frac{M(x,a)}{|a|}=\infty.$$
 Moreover the conjugate function $M^*$ is defined as 
$$M^*(x,b)=\sup_{a\in{\mathbb R}^d}(b\cdot a-M(x,a)).$$
%The last important property is the formulation 
%with help of a maximal monotone graph, not necessarily being a graph of $A$ as a function of $\nabla u$. 
%
The graph is expected to satisfy for almost all $(t,x)\in Q$ the following set of assumptions:
\begin{enumerate}
\item[{ (A1)}] $\mA$  comes through the origin.%$(\b0,\b0)\in\mA$.
\item[{ (A2)}] $\mA$  is a monotone graph, namely
$$
(A_1-A_2)\cdot (\xi_1-\xi_2) \ge 0 \quad \textrm{ for all } (\xi_1, A_1),(\xi_2,A_2)\in \mA(t,x)\,.
$$
\item[{ (A3)}] $\mA$  is a maximal monotone graph. Let $(\xi_2, A_2)\in \R^d \times \R^d$.
\begin{equation*}\begin{split}
&\textrm{If } ( {A_1} - A_2)\cdot( {\xi_1} - \xi_2) \geq 0 \quad \textrm{ for all }
({\xi_1}, A_1) \in \mA(t,x)\\&
 \textrm{ then } (\xi_2, A_2) \in \mA(t,x).
\end{split}\end{equation*}
\item[{ (A4)}]  $\mA$  is an  {\it $M-$ graph.} There are non-negative $k\in
L^1(Q)$,   $c_*>0$ and $N$-function $M$ such that
\begin{equation*} %\label{growthA}
A \cdot \xi \geq -k(t,x) +c_*(M(x,\xi) + M^*(x,A)) 
\end{equation*}
 for all $ (\xi,A)\in\mA(t,x).$ 
\item[{ (A5)}]  The existence of a measurable selection.  Either there is $\tilde A:Q\times \mathbb{R}^{d} \to
\mathbb{R}^{d}$ such that $(\xi,
\tilde A(t,x,\xi)) \in \mA(t,x)$ for all $\xi \in \Rd$ and
$\tilde A$ is measurable, 
 or there is $\tilde \xi:Q\times \mathbb{R}^{d} \to \mathbb{R}^{d}$ such
that $(\tilde\xi(t,x,A), A) \in \mA(t,x)$ for all $A \in \R^{d}$ and $\tilde\xi$ is measurable.
\end{enumerate}

Let us shortly refer again to the classical Leray-Lions operators. Within the setting  presented 
above  we would use the $N-$function
$M(a)=|a|^p$ with the conjugate function $M^*(a)=|a|^{p'}$, with $1/p+1/p'=1$. 

As we allow also for $x$ dependence, then the presented framework caputres also  the case of growth conditions in variable exponent case, namely
$M(a)=|a|^{p(x)}$. The further generalization is the anisotropic character and  functions different than only polynomials, hence the following example is acceptable 
$M(x,a)=a_1^{p_1(x)}\ln(|a|+1)+e^{a_2^{p_2(x)}}-1$ for $a=(a_1,a_2)\in\R^2$.
All the functions having a growth essentially different than polynomial (e.g. close to linear or exponential)
yield additional analytical difficulties and significantly constrain good properties of corresponding function spaces (like separability or reflexivity, or density of compactly supported smooth functions).  
We shall now discuss this issue in more detail. For this reason let us recall some definitions.
By the generalized Musielak-Orlicz class ${\mathcal L}_M(Q)$ we mean 
 the set of all measurable functions
 $\xi:Q\to\R^{d}$ for which the modular  $$\rho_{M,Q}(\xi)=\int_Q M(x,\xi(t,x)) \,dx\,dt$$ is finite. 
 By $L_M(Q)$ we mean the generalized Orlicz space which is the set of
all measurable functions
 $\xi:Q\to\R^{d}$ for  which $\rho_{M,Q}(\alpha\xi)\to0$ as $\alpha\to
 0.$
 %  $$\int_Q M(x,\lambda \xi(x))dxdt\to0\quad {\rm as}\,\lambda\to
 %0.$$
This is a Banach space with respect to the
 norm
$$\|\xi\|_M=\sup\left\{\int_Q \eta\cdot\xi dxdt : \eta\in L_{M^\ast}(Q),\int_Q M^\ast(x,\eta)
\,dx\,dt\le1\right\}.$$ 
%or the equivalent Luxemburg norm
 %$$\|\xi\|=\inf\left\{\lambda>0:\quad\int_Q
 %M\left(x,\frac{\xi}{\lambda}\right)dxdt\le 1\right\}.$$
 All over in the above definitions we used the notion of generalized Musielak-Orlicz spaces. Contrary to the classical Orlicz
 spaces we capture the case of $x-$dependent $N-$functions as well as functions dependent on the whole vector, not
 only on its absolute value (i.e. anisotropic).
Moreover, By $E_M(Q)$ we mean the closure of  bounded functions in
$L_M(Q)$. The space $L_{M^\ast}(Q)$ is the dual space of
$E_M(Q)$.
 A sequence $z^j$ is said to converge modularly to $z$ in $L_M(Q)$ if
there exists $\lambda>0$ such that
$$\rho_{M,Q}\left(\frac{z^j-z}{\lambda}\right)\to0$$
which is  denoted   by $z^j\modular{M} z$.
The basic estimates which we will frequently use in a sequel are 
the  H\"older inequality
\begin{equation}\label{hoelder}
\int_Q \xi \eta \,dx\,dt\le c\|\xi\|_M\|\eta\|_{M^*}
\end{equation}
and the  Fenchel-Young inequality
\begin{equation}\label{F-Y}
|\xi\cdot\eta|\le M(x,\xi)+M^*(x,\eta). 
\end{equation}
The essence of our considerations is  the lack of the assumption of $\Delta_2-$con\-di\-tion. We say 
that $M$ satisfies $\Delta_2-$condition if there exists a constant $c>0$ and a summable function $h$
such that
\begin{equation}\label{delta2}
M(x,2a)\le c M(x,a)+h(x)
\end{equation}
for all $a\in\R^d$. If $M$ satisfies \eqref{delta2} then $L_M(Q)$ is separable and compactly supported smooth functions are dense in strong topology. If additionally $M^*$ satisfies \eqref{delta2} then $L_M(Q)$ is reflexive. Notice that none of these assumptions is made in the present paper. For this reason the notion of modular topology and the issue of density of compactly supported smooth functions with respect to the modular topology are of crucial meaning. 
The basic properties which are mentioned above of anisotropic Musielak-Orlicz spaces were discussed and proved 
in~\cite{GwMiSw2013}.

As the density arguments become an essential tool, then the dependence of an $N-$function on $x$ 
becomes the significant constraint. The problem arises when we try to estimate uniformly the convolution operator. To handle this obstacle, we need some regularity with respect to the space variable. More precisely,  
we will assume that the function  $M$ satisfies  the following properties:
\begin{itemize}
% \item[(M1)] the function $M^*$ satisfies $\Delta_2$ condition in $\Omega$, i.e., 
% \begin{equation}\label{delta}
% M^*(t,x,2\xi)\le c M^*(t,x,\xi)+h(t,x) 
% \end{equation}
% for all $\xi\in\R^d$, a.a. $(t,x)\in(0,T)\times \Omega$ and some $h\in L^1((0,T)\times \Omega)$
%and  there exists a constant $H_1>0$ such that for all $ t,s\in[0,T],  |t-s|\le\frac{1}{2}$
%\begin{equation}\label{log-t}
%\frac{M(x,\xi)}{M(s,x,\xi)}\le |\xi|^\frac{H_1}{\log\frac{1}{|t-s|}}
%\end{equation}
\item[(M)] there exists a constant $H>0$ such that for all $x,y\in \Omega, |x-y|\le\frac{1}{2}$
and for all $\xi\in\R^d$ such that $|\xi|\ge 1$
\begin{equation}\label{log}
\frac{M(x,\xi)}{M(y,\xi)}\le |\xi|^\frac{H}{\log\frac{1}{|x-y|}}.
\end{equation}
Moreover, 
for every  bounded measurable set $G$  and every $z\in\R^d$
\begin{equation}\label{int}
\int_G M(x,z)< \infty.
\end{equation}
\end{itemize}
%
%We consider the problem of existence of weak solutions to the initial boundary value problem
%(\ref{P1})--\Ref{P3}. In the present paper we assume a full anisotropy of the ${\cal N}$--function, namely it depends on 
%the vector-valued argument and we do not put any assumptions on the growth of an 
%${\cal N}-$ function. Indeed we do not assume that  neither $M$ nor $M^*$ satisfies the so-called 
%$\Delta_2$ condition\footnote{We say that an ${\cal N}$--function $M$ satisfies $\Delta_2$--condition if for some constant $C>0$ it holds $M(2\xi)\le CM(\xi)\quad\mbox{for all }\xi\in\R^{d}.$
%}. The following pair of functions can serve as an  example   in $d=2$
%$$M(\xi)=M(\xi_1,\xi_2)=e^{|\xi_1|}-|\xi_1|-1+(1+|\xi_2|)\ln(1+|\xi_2|)-|\xi_2|,$$
%$$M^*(\eta)=M^*(\eta_1,\eta_2)=(1+|\eta_1|)\ln(1+|\eta_1|)-|\eta_1|+e^{|\xi_2|}-|\eta_2|-1.$$
%The result in anisotropic case is not a straightforward extension of the isotropic one. The new difficulty arising here concerns the density of compactly supported smooth functions w.r.t. the modular topology of the gradients. The detailed analysis of this issue appears in Section \ref{closures}. 
%
%We can serve some nontrivial examples of the operator  $A$:
%\begin{itemize}
%\item $A(t,x, \xi)=a(t,x)\xi \exp(|\xi|^2)$ where $0<c_1\le a(t,x)\le c_2<\infty$
%\item  $A(t,x, \xi)=a(t,x)\xi \ln(1+|\xi|)$ where $0<c_1\le a(t,x)\le c_2<\infty$
%\item $A(t,x, \xi_1,\xi_2)=a_1(t,x)\xi_1 \exp(|\xi_1|^2)+a_2(t,x)\xi_2 \ln(1+|\xi_2|)
%$ where $0<c_1\le a_i(t,x)\le c_2<\infty$ for $i=1,2$
%\end{itemize}
% In \cite{Donaldson} the operator $A$ was assumed  to be
%an elliptic second-order operator in divergence form and monotone.
%The growth and coercivity conditions were more general than the
%standard growth conditions in $L^p$, namely the ${\cal N}$--function
%formulation was stated\footnote{Assumption {\bf (A2)} is in
%fact a generalization of the growth and coercivity conditions
%assumed in \cite{Donaldson} for the case of $M$ dependent on
%vector-valued argument.}. Under the assumptions
%  on the ${\cal N}$--function $M$: $\xi^2<<M(|\xi|)$ (i.e., $\xi^2$
%grows essentially less rapidly than $M(|\xi|)$) and $M^\ast$
%satisfies a $\Delta_2$--condition, the existence results to
%\Ref{P1}--\Ref{P3} was established.
% The restrictions  on the
%growth of $M$ were abandoned in \cite{ElmahiMeskine}, but still  $M$ had an isotropic character.  
%
%%The elliptic problems were studied e.g. in \cite{Gossez1}.
%
%The review paper \cite{MustonenTienari} summarizes the
%monotone-like mappings techniques in Orlicz and Orlicz--Sobolev
%spaces\footnote{$W^mL_M$ is the Orlicz--Sobolev space of
%functions in $L_M$ with all distributional derivatives up
%to the order $m$ in $L_M$.}. 
%%The authors need the
%%essential modifications of the notions like monotonicity,
%%pseu\-do\-mo\-no\-to\-ni\-ci\-ty, operators of type $(M)$,
%%$(S_+)$, et al. The reason is  that Orlicz--Sobolev spaces are not
%%reflexive in general. Moreover, the nonlinear differential
%%operators in divergence form with standard growth conditions are
%%neither bounded nor everywhere defined.
%
%, which is reasonable from the physical point of
%view, since the electric field is a vector.
%The spaces with $M$ dependent on a vector-valued argument were
%extensively studied in \cite{Turett} and also in
%\cite{HudzikKaminska, Wisla1}.
%Assumptions (M1)--(M2) show that for those $x$ for which the conjugate function $M^*$ has moderate growth the existence of solutions can be easily shown. In those points where the growth starts to be faster (i.e. lack of $\Delta_2$ condition) we need to prescribe the regularity w.r.t. $t$ and $x$, which is the condition (M1). 
%The basic energy estimates  provide  a solution is in a generalized Musielak-Orlicz space. By generalized we mean that the $N-$function depends on the vector-valued argument
%(not  
%only on $|\xi|$, but on the whole vector $\xi$).
%
%Before collecting all definitions of the spaces we discuss some examples of where this framework can be used and related results. 
%
%wstawic
%Contrary to \cite{Musielak} we consider
% the ${\cal N}$--function $M$  not  dependent
%only on $|\xi|$, but on the whole vector $\xi$.

Below we formulate the definition and then state the existence theorem which is the main result of the present paper.  We shall use the following notation: by ${\cal C}_c^\infty(\Omega)$ we denote the space
of infinitely differentiable compactly supported functions in $\Omega$. Let $p\le 1\le \infty$ and 
$k\in\N$, then we denote by $(L^p(\Omega),\|\cdot\|_{L^p(\Omega)})$ the Lebesgue spaces and by 
$(W^{k,p}(\Omega), \|\cdot\|_{W^{k,p}(\Omega)})$ the Sobolev spaces.  By $W^{k,p}_0(\Omega)$  we
mean the closure of ${\cal C}_c^\infty(\Omega)$  with respect to the norm $ \|\cdot\|_{W^{k,p}(\Omega)}$ and $W^{-k,p'}(\Omega)$ with $1/p+1/p'=1$ denotes its dual space.  
Moreover we use the notation
${\cal C}_{weak}(0,T;L^2(\Omega))$ for the space of all $\varphi\in L^\infty(0,T;L^2(\Omega))$ which satisfy $(\varphi(t),v)\in{\cal C}([0,T])$ for all $v\in{\cal C}(\bar\Omega)$. 



\begin{definition}\label{d:1}
Assume that
$
u_0\in L^{2}(\Omega)
$ and $f\in L^\infty(Q).$ 
We say that $(u,A)$ is weak solution to \eqref{P1a}-\eqref{P3a}  if
\begin{equation}\label{LM}
u\in L^\infty(0,T;L^2(\Omega)),  \nabla u \in L_M(Q),\ 
A\in  L_{M^\ast}(Q)
\end{equation}
%&\lim_{t \to 0_+}\|u(t)-u_0\|_2^2 =0,\\
and \begin{equation}\label{Cweak}
u \in  {\cal C}_{\textrm{weak}}(0,T; L^{2}(\Omega)).
\end{equation}
Moreover, the following identity 
\begin{equation}
\begin{split}
\int_Q \left(-u \varphi_t +A
\cdot \nabla \varphi \right) \,dx \,dt+\int_\Omega u_0(x)\varphi(0,x) \,dx
=\int_Qf\varphi \,dx\,dt,
 %&\qquad \textrm{ for all } \varphi\in{\cal D}(\Omega) %\bw\in W^{1,1}  \textrm{ such that } \nabla \bw\in L^{\infty}(\Omega) 
% \textrm{ and a.a. } t\in (0,T),
\end{split} 
\end{equation}
is satisfied  for all  $\varphi\in{\cal C}_c^\infty((-\infty,T)\times\Omega)$
and
\begin{equation}
\left ( \nabla u ((t,x)), A(t,x) \right ) \in \mA(t,x)  \textrm{ for a.a. }
(t,x) \in Q.
\end{equation}
\end{definition}

%The space $Z^M_0$ appearing in the theorem 
%is defined in Section~\ref{closures}. 
%For sake of the existence result under assumption (M1) we observe that \eqref{delta} implies there exists some $q>1, c_1, c_2>0$ such that\marginpar{jakas funkcja jeszcze od t,x}
%\begin{equation}\label{est-M*}
%M^*(t,x,\xi)\le c_1|\xi|^q
%\end{equation} 
%and hence 
%\begin{equation}\label{est-M}
%M(x,\xi)\ge c_2|\xi|^{q'}.
%\end{equation} 
%This exponent $q$ will be used to conclude information on the boundedness of $u.$ Indeed,  $W^{1,q'}(\Omega)\subset L^r(\Omega),$ where 
%\begin{equation}\label{s}
%r=\frac{dq'}{d-q'}.
%\end{equation}
%%\marginpar{dopracowac lematem Lionsa Aubin, ale to raczej dalej.}

%\begin{theorem}\label{main}
%Let $M$ be an ${\cal N}$--function satisfying  (M1) and let  $A$ satisfy
%conditions (A1)--(A5). Given $f\in L^\infty(Q)$ %{\frac{d+q}{d}}(0,T;L^{r'}(\Omega)) $, where $1/r+1/r'=1$ and $r$ is given by \eqref{s}
%%in the form $f=\div F$ with $F\in\R_{sym}^{d\times d}$, $f\in
%%L^1(0,T;V_s^\ast)$  (LUKA- zal na F nie na f)
%and
% $u_0\in
%L^2(\Omega)$ 
%%and $f\in E_B$, where an $N-$function $B$ is given by \eqref{B} 
%there exists $u\in Z^M_0$ such that 
%\begin{equation}
%\int_Q \left(-u \varphi_t +A(t,x,\nabla u)
%\cdot \nabla \varphi \right) dx dt+\int_\Omega u_0(x)\varphi(0,x) dx=\int_Qf\varphi dxdt
%\end{equation}
%holds for all $\varphi\in{\cal C}^\infty_c((-\infty,T)\times\Omega).$
%\end{theorem}

\begin{theorem}\label{main2}
Let $M$ be an ${ N}$--function satisfying (M) and let  $A$ satisfy
conditions (A1)--(A5). Given $f\in L^\infty(Q) $
%in the form $f=\div F$ with $F\in\R_{sym}^{d\times d}$, $f\in
%L^1(0,T;V_s^\ast)$  (LUKA- zal na F nie na f)
and
 $u_0\in
L^2(\Omega)$ 
%and $f\in E_B$, where an $N-$function $B$ is given by \eqref{B} 
there exists  a weak solution to \eqref{P1a}-\eqref{P3a}.
\end{theorem}
%
%The assumption on the function $f$ in Theorem~\ref{main2} is rather strong. Note that to handle the right-hand side in an appropriate Musielak-Orlicz space we would need to be able to claim that $u$ is bounded in some Musielak-Orlicz space. For this purpose the results on Sobolev embeddings are needed. These kind of results, even in anisotropic case are available for Orlicz spaces (the ${\cal N}-$function independent of $t,x$) as well as in generalized Lebesgue spaces $L^{p(x)}$, see e.g.~\cite{Cianchi,Di2004, EdRa2000}. In this paper we shall not concentrate on the issue of Sobolev embeddings, hence the case of the simplest
%right-hand side is considered. 
%
The current paper provides complementary studies to the results presented in~\cite{Sw2013}. Here we also consider the problem of existence of weak solutions to the parabolic problem including multivalued terms. However, the essential differerence consists in the properties of an 
$N-$function describing the growth conditions of graph $\mA$. In~\cite{Sw2013} we concentrated on 
the case with time-dependent $N-$function. This required more delicate approximation theorem 
and excluded the possibility of anisotropic functions. 
%Moreover, the regularity of the boundary was higher. 
The studies presented here do not extend the results of the previous paper, but are parallel to them. We decided to omit here the dependence on time of an $N-$function, but added the possibility of  anisotropic behaviour. 


The anisotropic parabolic problems were conisdered also  in~\cite{GwSw2010}. 
This was however much simpler situation, namely the studies concerned an equation and 
the $N-$function was assumed to be homogeneous in space. The anisotropic and space-inhomogeneous problems, however in 
slightly different setting, namely in the case of systems describing flow of non-Newtonian fluids were considered in~\cite{GwSw2008, GwSw2008a, GwSwWr2010, Wr2010}. The authors assumed $\Delta_2-$condition on the conjugate $N-$function. The simplified problem, namely the generalized Stokes equation, in the case omitting the $\Delta_2-$condition on the conjugate $N-$function was considered in~\cite{GwSwWr2012}.

The approach of maximal monotone graphs also to problems arising in fluid mechanics was undertaken in  \cite{ BuGwMaSw2009, 
GwMaSw2007} for the $L^p$ setting and in  \cite{BuGwMaRaSw2012,BuGwMaSw2012} for the setting in Orlicz spaces. The latter ones however were restricted to classical Orlicz spaces with the assumption that $\Delta_2-$condition was satisfied.

%\cite{GwZa2005, GwZa2005a,GwZa2007}. Another novelty lies in the function space of  solutions. 

Most of the earlier results on existence of solutions to parabolic problems in non-standard setting concern
the case of classical Orlicz spaces, see e.g.~\cite{Donaldson}
and  later studies  of Benkirane, Elmahi and Meskine, cf.~\cite{BeEl1999, ElmahiMeskine, ElMe2005}.
%Landes~\cite{}, Landes and Mustonen~\cite{}.
 All of them concern the case of  an $N-$function dependent only on $|\xi|$ without the dependence on $x$.  

The paper is organized as follows: Section 2 contains the proof of Theorem~\ref{main2}, Section 3 
is devoted to the problems of density of compactly supported smooth functions with respect to 
the modular convergence. In the appendix we include some facts, which are used in the sequel and we refer 
to their proofs. 


\section{Existence of solutions }

The current section contains a proof of Theorem~\ref{main2}. The construction of an approximate problem follows in two steps.
By (A5) there exists a measurable selection $\tilde A:Q\times \R^{d}\to\R^{d}$ of the graph $\mA$.  Obviously, each such a selection $\tilde A$ defined
on $\R^{d}$, is monotone and due to (A4) satisfies
\begin{equation}\label{316}
 \tilde A(t,x,\xi)\cdot\xi\ge -k(x,t)+ c_*(M(x,\xi)+M^*(x,\tilde A(t,x,\xi))\mbox{ for all }\xi \in \R^{d}.
\end{equation} We mollify $\tilde A$  with a smoothing kernel
and then construct the finite-dimensional problem by means of Galerkin method.
 Indeed, let 
 \begin{equation}\label{S}
 S\in
{\cal C}^\infty_c(\R^{d}), \  \int_{\R^{d}} S(y)\;dy=1, \, S(y)=S(-y), \   S_\varepsilon(y):=1/\varepsilon^{d}S(y/\varepsilon)
\end{equation}
  with ${\rm supp }\, S$ in a unit ball $B(0,1)\subset\R^d$ and define
 \begin{equation}\label{Te}
 A^\varepsilon (t,x,\xi):=(\tilde A* S_\varepsilon )(t,x,\xi)=\int_{\R^d}\tilde A(t,x,\zeta)
 S_\varepsilon (\xi-\zeta)\, d\zeta.
%=\int_{\Rd}\tilde A(t,x,\bzeta)d\mu^\varepsilon _{\xi}(\bzeta).
\end{equation}
 Using the convexity of $M$ and $M^*$ and the Jensen inequality allows to conclude that the approximation $ A^\varepsilon$ satisfies a condition analogous to  \eqref{316}, namely
\begin{equation}\label{319}
A^\vep \cdot \nabla u \geq -k(t,x) +c_*(M(x,\nabla u) + M^*(x,A^\vep)). %\qquad \textrm{ for all } (\nabla u,A^\vep)\in\mA(t,x)\,.
\end{equation}
For the proof of analogous estimate for the approximation in case of polynomial conditions see~\cite{GwMaSw2007} and also~\cite{GwZa2007}. 

The assumption (A5) included either the possibility of  existence of a selection  $\tilde A$,
as was presented above, or  existence of a selection $\tilde\xi:Q\times\R^d\to\R^d$, such that 
$(\tilde \xi(t,x,A),A)$ is for all $A\in\R^d$ in the graph $\mA$. In the second case we would define
\begin{equation}
\xi^\vep(A):=(\tilde\xi*S_\vep)(t,x,A)+\vep A.
\end{equation}
Such a definition provides that the function $A\mapsto\xi^\vep(A)$ is invertible. Note that  since $\vep A\cdot A\ge 0$
one can show that for the pair $(\xi^\vep(A), A)$ an analogue of \eqref{319} holds, and consequently also for 
$(\xi, (\xi^\vep)^{-1}(\xi))$.
 Thus we may  define $A^\vep$ as follows
\begin{equation}
A^\vep:=(\tilde\xi*S_\vep+\vep\,  I\!d)^{-1}. 
\end{equation}
One  proceeds further analogously to the previous situation. In the sequel we present the proof for the case when there exists a selection
$\tilde\xi$ and $A^\vep$ is given by \eqref{Te}.

Consider now the basis consisting of eigenvectors of the Laplace operator  with Dirichlet boundary condtion and let  $u^{ \varepsilon,n}$ be the solution to the   finite dimensional problem with function~$A^\varepsilon$, namely 
$u^{ \varepsilon,n}(t,x):=\sum_{i=1}^nc_i^{ \varepsilon,n}(t)\omega_i(x)$ which solves the following system
\begin{equation}\label{aproksymacja}
\begin{split}
(u^{ \varepsilon,n}_t, \omega_i)+(A^\varepsilon(t,x,\nabla u^{ \varepsilon,n}),\nabla \omega_i)
=\langle f,\omega_i\rangle, \qquad i=1,\ldots,n, \\
u^{ \varepsilon,n}(0)=P^nu_0
\end{split}\end{equation}
where $P^n$ is the orthogonal projection of $L^2(\Omega)$ on the ${\rm span}\,\{\omega_1,
\ldots,\omega_n\}.$
 Let $Q^s:=(0,s)\times \Omega$ with 
$0<s<T$.
%For further estimates observe that in our particular case the Fenchel-Young inequality holds
%with an equality sign, namely
%\begin{equation}
%D\psi(\xi)\cdot\xi=\psi(\xi)+\psi^*(D\psi(\xi))
%\end{equation}
%for all $\xi\in\R^d$. Hence  
%In the standard manner 
%we conclude that  for $0<s<T$
%\begin{equation}\label{Galerkin}\begin{split}
%\frac{1}{2}\|u^{ \varepsilon,n}(s)\|^2_{2}+\int_{Q^s} A^\varepsilon(t,x,\nabla u^{ \varepsilon,n})\cdot \nabla u^{ \varepsilon,n}\, dx\,dt\\=\frac{1}{2}
%\|u^{ \varepsilon,n}(0)\|^2_{2}+\int_0^s\langle f,u^{ \varepsilon,n}\rangle dt
%\end{split}\end{equation}
%holds. 
%We estimate the right-hand side as follows %(using  H\"older's inequality and Young's inequality)
%\begin{equation}\begin{split}
%|\int_0^s\langle f,u^{ \varepsilon,n}\rangle dt|&\le \int_0^s\|f\|_{\infty}\|\u^{ \varepsilon,n}
%\|_{1}dt\le
%c \int_0^s\|f\|_{\infty}\|\nabla u^{ \varepsilon,n}\|_{1} dt\\
%&\le c\|f\|_{\infty}\int_0^s\|\nabla u^{ \varepsilon,n}\|_{1}dt \\
%&\le \frac{c_*}{2}\int_0^s\|\nabla  u^{ \varepsilon,n}\|_{1} dt+K(c_*,Q)\|f\|_{\infty}\\
%&\le  \frac{c_*}{2} \int_Q M(x,\nabla u^{ \varepsilon,n})dxdt+K\|f\|_{\infty}.\\
%%&\le \frac{c_*}{2}\int_Q M(t,x,\nabla u)dxdt+C\|f\|_{s'}^q.
%\end{split}\end{equation}
%\begin{equation}\begin{split}
%|\int_0^s\langle f,u^{\varepsilon,n}\rangle dt|&\le \int_0^s\|f\|_{r'}\|\u\|_{r}dt\le
%c \int_0^s\|f\|_{r'}\|\nabla u\|_{q'} dt\\
%&\le c\left(\int_0^s\|\nabla u\|_{q'}^{q'}dt\right)^{1/q'} \left(\int_0^s\|f\|_{r'}^q dt\right)^{1/q}\\
%&\le \frac{c_*}{2c_2}\int_0^s\|\nabla  u\|_{q'}^{q'} dt+K(c_*,c_2,q,Q)\int_0^s\|f\|_{r'}^q dt\\
%&\le  \frac{c_*}{2} \int_Q M(t,x,\nabla u)dxdt+K\int_0^T \|f\|_{r'}^q dt.\\
%%&\le \frac{c_*}{2}\int_Q M(t,x,\nabla u)dxdt+C\|f\|_{s'}^q.
%\end{split}\end{equation}
%Simple computations yield that  $q/r'=(d+q)/d,$ shall lead us to the uniform estimates. 
Using \eqref{319} allows to conclude
\begin{equation}\begin{split}
\sup\limits_{s\in(0,T)}\|u^{ \varepsilon,n}(s)\|_{L^2(\Omega)}^2&+c_*\int_QM(x,\nabla u^{ \varepsilon,n})
+M^*(x,A^\varepsilon(t,x,\nabla u^{ \varepsilon,n})) \,dx\,dt
\\&
\le c(\|u_0\|_{L^2(\Omega)}^2+\|f\|_{{L^\infty(Q)}}+\int_Q k \,dx\,dt
).\end{split}\label{energy-eps}
\end{equation}
%From \eqref{aproksymacja} and \eqref{energy-eps} it follows that 
%$c_i^{\vep,n}(t) $ is bounded in $L^\infty([0,T])$ and $\frac{d}{dt}c_i^{\vep,n}(t)$ is bounded in $L_{M^*}([0,T])$, hence uniformly integrable in $L^1([0,T])$ and there exists a monotone, continuous $L:\R_+\to\R_+$, with $L(0)=0$ such that
%$$\left|\int_{s_1}^{s_2}\frac{d}{dt}c_i^{\vep,n}(t)\,dt\right|\le L(|s_1-s_2|)$$
%and thus 
%$$|c_i^{\vep,n}(s_1)-c_i^{\vep,n}(s_2)|\le L(|s_1-s_2|).$$
%Hence by the Arzel\`a-Ascoli theorem we are able to conclude an existence of uniformly convergent subsequence 
%$\{c_i^{\vep_k,n}\}.$
%The limit passage with $\varepsilon\to 0$ is done on the level of finite-dimensional problem and follows the similar lines as in~\cite{BuGwMaSw2012}. Nevertheless, we include the details for completeness. 
 In a consequence of \eqref{energy-eps}  there exists a subsequence (labelled the same) such that
 % \marginpar{zastanowic sie nad tymi zbieznosciami}
   \begin{equation}\label{ep-to-zero}\begin{split}
  \nabla u^{ \varepsilon,n}\weakstar  \nabla u^{ n}\quad &\quad\mbox{weakly-star  in}\quad L_{M}(Q),\\
A^\varepsilon (\cdot,\cdot,\nabla u^{ \varepsilon,n})\weakstar A^{ n} \quad &\quad\mbox{ weakly-star  in}\quad L_{M^*}(Q).
\end{split}\end{equation}
Moreover, from relation \eqref{aproksymacja} we concude the boundedness of the sequence 
$u^{ \varepsilon,n}_t$ in $L_{M^*}(Q)$ and hence up to the subsequence we have 
\begin{equation}
u^{ \varepsilon,n}_t\weakstar u^{ n}_t\quad\quad\mbox{ weakly-star  in}\quad L_{M^*}(Q).
\end{equation}

Further we observe that  \eqref{aproksymacja}  implies  that 
$\frac{d}{dt}c_i^{\vep,n}(t)$ is bounded in the space $L_{M^*}([0,T])$, what implies the  uniform integrability in $L^1([0,T])$. Consequently  there exists a monotone, continuous $L:\R_+\to\R_+$, with $L(0)=0$ such that for all $s_1,s_2\in (0,T)$
$$\left|\int_{s_1}^{s_2}\frac{d}{dt}c_i^{\vep,n}(t)\,dt\right|\le L(|s_1-s_2|)$$
and thus the sequence $c_i^{\vep,n}$ is uniformly equicontinuous 
$$|c_i^{\vep,n}(s_1)-c_i^{\vep,n}(s_2)|\le L(|s_1-s_2|).$$
From   \eqref{energy-eps} we conclude that  $c_i^{\vep,n}(t) $ is bounded in $L^\infty([0,T])$ 
and hence by the Arzel\`a-Ascoli theorem there  exists a uniformly convergent subsequence 
$\{c_i^{\vep_k,n}\}$ in ${\cal C}([0,T])$  and taking into account the regularity of the basis 
$\{\omega_i\}_{i=1}^n$ we conclude

\begin{equation}\label{mocna}
 u^{ \varepsilon,n}\to  u^{ n}\quad \quad\mbox{strongly in}\quad {\cal C}([0,T];{\cal C}^1(\overline\Omega)).
\end{equation}
The limit passage with $\varepsilon\to 0$ is done on the level of finite-dimensional problem. It follows the similar lines as in~\cite{BuGwMaSw2012}, however we shall recall the main steps. 
%Note that the strong convergence follows directly from the Arzel\`a-Ascoli theorem and appropriate estimates 
%obtained with help of  equation
%\eqref{aproksymacja}.
%After passing to the limit we conclude the following limit identity
%\begin{equation}\label{limit}
%-\int_Q u \varphi_t dx dt +\int_Q{\chi} \cdot \nabla\varphi dx dt+\int_\Omega u_0(x)\varphi(0,x) dx=0
%\end{equation}
Using \eqref{ep-to-zero}-\eqref{mocna} we obtain the  following limit problem 
\begin{equation}\label{Galerkin}
\begin{split}
(u^{ n}_t, \omega_i)+(A^{ n},\nabla \omega_i)
=\langle f,\omega_i\rangle, \qquad i=1,\ldots,n, \\
u^{ n}(0)=P^nu_0.
\end{split}\end{equation}
To complete the limit passage we need to provide that  
\begin{equation}\label{in-graph}
( \nabla u^{ n}, A^{ n})\in{\cal A}.
\end{equation}
Following \cite{BuGwMaSw2012} and also \cite{Sw2013}, with simple algebraic tricks and estimates which are not included in the present paper,
we conclude that for all $B\in\Rd$ and for a.a. $(t,x)\in Q$
\begin{equation}
(A^{ n}-\tilde A(t,x,B)) \cdot (\nabla u^{ n}-B)\ge0 \,.
\end{equation}
Hence, using  the equivalence
  of $(i)$ and $(ii)$ in  Lemma~\ref{LS*}, we arrive to \eqref{in-graph}.
% Since $\tilde A$ is the selection of the graph, according to Lemma~\ref{LS*} (a2)  for all $\zeta, B\in\Rd$ and a.a. $(t,x)\in Q$ it holds
%\begin{equation}\label{mon}
%(\tilde A(t,x,\zeta)-\tilde A(t,x,B))\cdot(\zeta-B)\ge0.
%\end{equation}
%We shall  add and subtract the term $(\tilde A(t,x,\zeta)-\tilde A(t,x,B))\cdot \nabla u^{ \varepsilon,n}$ and
%then integrate  with respect to the probability measure, which  has the
%density $K^\varepsilon (\nabla u^{ \varepsilon,n}-\zeta)$ and obtain  that
%\begin{equation}\begin{split}\label{mon1}
%\int_{\Rd}&(\tilde A(t,x,\zeta)-\tilde A(t,x,B))\cdot (\nabla u^{ \varepsilon,n}-B)K^\varepsilon (\nabla u^{ \varepsilon,n}-\zeta)d\zeta\\
%&\ge \int_{\Rd}(\tilde A(t,x,\zeta)-\tilde A(t,x,B))\cdot (\nabla u^{ \varepsilon,n}-\zeta)K^\varepsilon (\nabla u^{ \varepsilon,n}-\zeta)d\zeta.
%\end{split}\end{equation}
%Consider $|\zeta|\le \|\nabla u^{ \varepsilon,\eta}\|_\infty
%+ \vep$. The difference $(\tilde A(t,x,\zeta)-\tilde A(t,x,B))$ can be estimated  by a constant
%dependent only  on $\bB$. Hence from  \eqref{mon1} we conclude that
%\begin{equation}\label{concl_mon}\begin{split}
%\left(\int_{\Rd}\tilde A(t,x,\zeta)K^\varepsilon (\nabla u^{ \varepsilon,n}-\zeta)d\zeta-\tilde A(t,x,B)\right)\cdot(\nabla u^{ \varepsilon,n}-B)\ge\\
%-C_n(B) \int_{\Rd}|\nabla u^{ \varepsilon,n}-\zeta|K^\varepsilon (\nabla u^{ \varepsilon,n}-\zeta)d\zeta.
%\end{split}\end{equation}
%Using the strong convergence \eqref{ep-to-zero}$_1$ we see that the right hand side of \eqref{concl_mon} vanishes as $\vep \to 0_+$ and  we get
%\begin{equation}\begin{split}
%\liminf\limits_{\vep\to0_+}\left(A^\vep(t,x,\nabla u^{ \varepsilon,n})-\tilde A(t,x,B)\right)\cdot(\nabla u^{ \varepsilon,n}-B)\ge0 \qquad \textrm{ for a.a. } (t,x)\in Q.
%\end{split}\end{equation}
%The strong convergence of  $\nabla u^{ \varepsilon,n}$ and weak-star  convergence of  
%$A^\vep(t,x,
%\nabla u^{ \varepsilon,n})$
%yields that for all $B\in\Rd$ and for a.a. $(t,x)\in Q$
%\begin{equation}
%(A^{ n}-\tilde A(t,x,B)) \cdot (\nabla u^{ n}-B)\ge0 \,.
%\end{equation}
%Thus, Lemma \ref{LS*} yields that 
%\begin{equation*}
%(\nabla u^{ n}(t,x),A^{ n}(t,x))\in\mA(t,x) \quad \textrm{ for a.a. } (t,x)\in Q\,.
%\end{equation*}
%%The remaining part of the proof shall concern the limit passage with $n\to\infty$ and is an essence 
%%of the limit procedure. 
Before passing to the limit with $n\to\infty$ we notice that 
in the same manner as before we obtain the estimates, which are uniform with respect to $n$, namely 
\begin{equation}\begin{split}
\sup\limits_{s\in(0,T)}\|u^{ n}(s)\|_{L^2(\Omega)}^2+\int_QM(x,\nabla u^{ n})+M^*(x,A^{ n}) \,dx\,dt\\
\le c(\|u_0\|_{L^2(\Omega)}^2+\|f\|_{L^\infty(Q)}+\|k\|_{L^1(Q)}). 
\end{split}\label{esty2}\end{equation}
Consequently  there exists a subsequence, labelled the same, such that 
  \begin{equation}\begin{split}
 \nabla u^{ n}\weakstar \nabla u\quad &\quad\mbox{weakly-star in}\quad L_M(Q),\\
  %\nabla u^{ n}\weakstar \nabla u\quad &\quad\mbox{weakly-star in}\quad L_M(Q),\\
  u^{ n}\weak  u\quad &\quad\mbox{weakly in}\quad L^1(0,T;W^{1,1}(\Omega)),\\
A^{ n}\weakstar A\quad &\quad\mbox{ weakly-star  in}\quad L_{M^*}(Q),\\
u^{ n}\weakstar u\quad&\quad\,\mbox{weakly-star in}\,\,L^\infty(0,T;L^2(\Omega)).\\
u^{ n}_t\weakstar u_t\quad&\quad\mbox{ weakly-star  in}\quad W^{-1,\infty}(0,T;L^2(\Omega)).\end{split}\label{228}\end{equation}
Using \eqref{228} we let $n\to\infty$ and conclude from \eqref{Galerkin} that the following  identity
\begin{equation}\label{limit}
 u_t  -\div {A} =  f
\end{equation}
holds in a distributional sense. 
Again, to complete the limiting procedure, we need to show that $(\nabla u, A)\in\mA(t,x)$. 
This case however requires more attention, contrary to the previous limit passage on the level of
fixed finite dimension $n$. The essence of this step is using the maximal monotonicity of the graph
$\mA$, in particular the property formulated in Lemma~\ref{Minty2}. As the  assumptions 
 \eqref{11}-\eqref{33} are obviously satisfied, then our attention shall be directed to \eqref{44}. 
For this aim we need to establish a strong energy inequality. Since testing \eqref{limit} with a solution is not possible, we first approximate it with respect to the space variable. 
By  Theorem~\ref{main1} there exists a sequence $v^j\in L^\infty(0,T;{\cal C}_c^\infty(\Omega
))$ such that 
\begin{equation}\label{zb_w_M}
\nabla  v^j\modular{M} \nabla u\ \mbox{ modularly in}\  L_M(Q) 
\ \mbox{and}\  v^j\to u \ \mbox{strongly in }\ 
L^2(Q).
\end{equation}
%and we can write
%\begin{equation}
%\frac{\partial v^j}{\partial t}=\div v_A^j+v_f^j\quad\mbox{and}\quad
%\frac{\partial u}{\partial t} =\div v_A+v_f
%\end{equation}
%with 
%\begin{equation}\label{zb_w_M*}
% v^j_A\modular{M^*}  v_A \ \mbox{ modularly in }\ L_{M^*}(Q)  
% \quad\mbox{and}\quad v^j_f\to v_f \mbox{ strongly in } L^2(Q).
% \end{equation}
%We are aiming to choose a test function
%\begin{equation}\label{test-function}
%{\cal K}_\varepsilon u=K_\vep*[(K_\vep*z)\cdot\bbbone_{(s_0,s)}]
%\end{equation}
% where $*$ is a convolution with respect to time.  Since ${\cal K}_\varepsilon u$ is not smooth with respect to $x$ we consider 
% an approximation $\{{\cal K}_\varepsilon u^j\}_{j\in\N}$ such that ${\cal K}_\varepsilon u^j$ converges weakly to 
% ${\cal K}_\varepsilon u$ in $L^1(0,T;W^{1,1}_{0}(\Omega))$ as $j\to\infty$. For passing to the limit in all terms of \eqref{limit}
% we introduce  the truncation operator
%$$\T_m(\eta)=\left\{\begin{array}{cc}
%\eta&|\eta|\le m,\\
%m \frac{\eta}{|\eta|}&|\eta|> m.
%\end{array}\right.$$
%%The sequence of functions $\{u^{\vep,j}\}$ is weakly convergent to
%%$u^\vep$ 
%Observe the following identity
%\begin{equation}\label{test-function_bis}\begin{split}
%\int_{s_0}^s \int_\Omega u_t \cdot {\cal K}_\varepsilon u^j\,dx\,dt&= \int_{s_0}^s \int_\Omega
%(\T_m(A)-A) \cdot \nabla{\cal K}_\varepsilon u^j\,dx\,dt\\
%&-\int_{s_0}^s \int_\Omega
%\T_m(A) \cdot \nabla {\cal K}_\varepsilon u^j \,dx\,dt\\
%%& +\int_{s_0}^s \int_\Omega
%%(\T_m(f)-f) \cdot {\cal K}_\varepsilon u^j\,dx\,dt\\
%&-\int_{s_0}^s \int_\Omega f \cdot {\cal K}_\varepsilon u^j dxdt.
%\end{split}
%\end{equation}
% We fix $k \in
%\N$ and using the convexity of $M$ and that $M^\ast$ satisfies
%$\Delta_2-$condition we estimate the integral
%\begin{equation}\begin{split}
%\int_{s_0}^s \int_\Omega |(\T_m(A)-A) \cdot \nabla {\cal K}_\varepsilon u^j|
%\,dx\,dt& \le \int_{s_0}^s \int_\Omega M^\ast(t,x,2^r(\T_m(A)-A))
%\,dx\,dt\\
%&+ \int_{s_0}^s \int_\Omega M(t,x,\frac{1}{2^r}\nabla {\cal K}_\varepsilon u^j) \,dx\,dt \\
%&\le C^r \int_{s_0}^s \int_\Omega M^\ast(t,x,\T_m(A)-A)
%\,dx\,dt\\
%&+ r\int_{s_0}^s \int_\Omega h(x) \bbbone_{\{|A(t,x)|>m\}} \,dx\,dt\\
%&+\frac{1}{2^r} \int_{s_0}^s \int_\Omega M(t,x,\nabla {\cal K}_\varepsilon u^j) \,dx\,dt.
%\end{split}\end{equation}
%Note that for each $\vep>0$
%\begin{equation}
%\lim \limits_{r \to \infty}
% \frac{1}{2^r}\sup\limits_{j \in \N}
%\int_{s_0}^s \int_{\Omega} M(t,x,\nabla {\cal K}_\varepsilon u^j) \,dx\,dt=0.
%\end{equation}
%Due to the convexity of $M^\ast$, symmetry and that
%$M^\ast(t,x,0)=0$ a.e. it holds
%$$M^\ast(t,x,\T_m(A)-A)\le M^\ast(t,x,A).$$
%Hence $A$ is an element of ${\cal L}_{M^\ast}(Q_T)$ and this
%yields by the Lebesgue convergence theorem   that
%$\int_{Q_T}M^\ast (t,x,\T_m(A)-A) dxdt$ converges to zero as
%$m\to\infty$. Hence
%\begin{equation}
%\lim \limits_{r \to \infty}\lim\limits_{m \to \infty} \int_{s_0}^s
%\int_{\Omega} C^rM^\ast(t,x,\T_m(A)-A) + r h(x)
%\bbbone_{\{|A(t,x)|>m\}} dxdt =0.
%\end{equation}
%
And hence we shall test with a  function of the form 
\begin{equation}\label{testfunkcja}
u^{j,\epsilon}= K^\epsilon\ast (K^\epsilon\ast v^j \bbbone_{(s_0,s)})
\end{equation}
   with $ K\in
{\cal C}_c^\infty(\R),$  $ K(\tau)= K(-\tau),\ 
\int_\R K(\tau)d\tau=1$ and  defining
$ K^\epsilon(t)=\frac{1}{\epsilon} K(t/\epsilon)$, $\epsilon<\min\{ s_0,T-s\}. $ Thus
\begin{equation}\label{335}
\int_{s_0}^{s}\int_\Omega (u*K^\epsilon)\cdot \partial_t(v^j*K^\epsilon)\,dx\,dt
=\int_Q A\cdot\nabla u^{j,\epsilon}\,dx\,dt-\int_Q fu^{j,\epsilon}\,dx\,dt.
\end{equation}
Because of \eqref{zb_w_M} we easily pass to the limit with $j\to\infty$.
Indeed, the left-hand side of \eqref{335} can be easily handled since this term can be reformulated
to $\int_Q((\partial_t K^\epsilon)*K^\epsilon*u) v^j\,dx\,dt$ and hence the limit passage is obvious.
Note that for all $0<s_0<s<T$ it follows
	\begin{equation}
	\begin{split}
 	 \int_{s_0}^{s}\int_\Omega
	(K^\epsilon\ast u)\cdot\partial_t(K^\epsilon\ast u) \,dx\,dt
	 &
= 
%\int_{s_0}^{s}\int_\Omega
	%( K^\epsilon\ast  u)\cdot\partial_t( K^\epsilon\ast  u) dxd t
	\int_{s_0}^{s} \frac{1}{2}\frac{d}{d t}\| K^\epsilon\ast
	u\|^2_{L^2(\Omega)}\,d t\\
	&=\frac{1}{2} \| K^\epsilon\ast u(s)\|^2_{L^2(\Omega)}-
	\frac{1}{2}\| K^\epsilon\ast u(s_0)\|^2_{L^2(\Omega)}.
	\end{split}
	\end{equation}
Passing  to the limit with $\epsilon\to0$ yields  for almost
all $s_0, s$, namely for all Lebesgue points of the function $ u(t)$ that the following identity
	\begin{equation}\label{vt}
	\lim\limits_{\epsilon\to0}
	\int_{s_0}^{s}\int_\Omega( u\ast  K^\epsilon)\cdot \partial_t ( u\ast K^\epsilon)=
	\frac{1}{2}\| u(s)\|_{L^2(\Omega)}^2-\frac{1}{2}\| u(s_0)\|^2_{L^2(\Omega)}
	\end{equation}
holds.
Observe now the term
	$$\int_0^T\int_\Omega  A \cdot( K^\epsilon\ast(( K^\epsilon\ast \nabla u)\ 	\bbbone_{(s_0,s)}))dxdt
	=\int_{s_0}^{s}\int_\Omega
	( K^\epsilon\ast A )\cdot( K^\epsilon\ast \nabla u) dxdt.$$
Both of the sequences $\{ K^\epsilon\ast A \}$ and
$\{ K^\epsilon\ast \nabla u\}$ converge in measure in $Q$ by Proposition~\ref{ae}.
 %due to Proposition~\ref{ae}. 
 Moreover%, since  $\nabla u\in L$
%and $ A \in {\mathcal{L}}_{M^*}(Q;\R^{n\times n}_{\rm sym})$ provide that 
	$$\int_Q (M(x,\nabla u)+ M^\ast(x, A ))\,dx\,dt< \infty	.
	$$ 
 Hence by Proposition~\ref{sup} we conclude 
that the sequences
$\{M^*\left( x,K^\epsilon\ast A \right)\}$ and 
$\{M\left( x,K^\epsilon\ast \nabla u\right)\}$
are uniformly integrable and with help of 
Lemma~\ref{modular-conv} we have
	\begin{equation*}
	\begin{split}
	 K^\epsilon\ast \nabla u
	&\modular{M}\nabla u\quad\mbox{modularly in }L_{M}(Q),\\
	 K^\epsilon\ast A 
	&\modular{M^\ast} A \quad\quad\mbox{modularly in }L_{M^*}(Q).
	\end{split}
	\end{equation*}
Applying Proposition~\ref{product} allows to conclude
	\begin{equation}\label{chi}
	\lim\limits_{\epsilon\to0}\int_{s_0}^{s}\int_\Omega
	( K^\epsilon\ast A )\cdot( K^\epsilon\ast \nabla u) dxdt
	=\int_{s_0}^{s}\int_\Omega  A \cdot \nabla u dxdt.
	\end{equation}
Passing to the limit with $\epsilon\to0_+$ in the right-hand side  is obvious. 
Hence for the moment we are able to claim that the following holds
\begin{equation}
\frac{1}{2}\|u(s)\|_2^2-\frac{1}{2}\|u(s_0)\|_2^2+\int_{Q^s}A\cdot \nabla u \,dx\,dt=\int_{Q_s} fu \,dx\,dt
\end{equation}
 for almost all $0<s_0<s<T$. For further considerations we need to know that the same holds for $s_0=0$, hence let us
 pass to the limit with $s_0\to0$. Thus, we need to establish that \eqref{Cweak} holds.  
  We shall observe that using the approximate equation we estimate 
the sequence $\{\frac{d\u^n}{dt}\}$ uniformly (with respect to $n$) in the space $L^1(0,T;W^{-r,2}(\Omega))$, where
 $r>\frac{d}{2}+1$.
Consider 
 $\varphi\in L^\infty(0,T;W^{r,2}_0(\Omega))$, $\|\varphi\|_{L^\infty(0,T;W^{r,2}_0)}\leq 1$ and observe that 

	\begin{equation*}
	\left\langle\frac{du^n}{dt}, \varphi \right \rangle=
	\left\langle\frac{du^n}{dt}, P^n\varphi\right \rangle
	= -\int_\Omega
	 A^n\cdot \nabla(P^n\varphi)\,dx\\
	 +\int_\Omega f\cdot P^n\varphi\,dx.
	\end{equation*}
Since the orthogonal projection is continuous in  ${W^{r,2}_0(\Omega)}$ and 
$W^{r-1,2}(\Omega)\subset L^\infty(\Omega)$ we estimate as follows 
	\begin{equation}\label{dt1}
	\begin{split}
	&\Big{|}\int_0^T\int_\Omega  A^n\cdot \nabla(P^n\varphi)
	dxdt\Big{|}\le \int_0^T\|
	 A^n\|_{L^1(\Omega)}\|\nabla(P^n\varphi)\|_{L^\infty(\Omega)}dt\\
	&\le c\int_0^T\| A^n\|_{L^1(\Omega)}\|P^n\varphi\|_{W^{r,2}_0}dt
	%\le c \int_0^T\|A(t,\cdot,\nabla u^n)\|_{L^1(\Omega)}\|\varphi\|_{W^{r,2}_0}dt\\ 
	\le c
	\| A^n\|_{L^1(Q)}\|\varphi\|_{L^\infty(0,T;W^{r,2}_0)}.
\end{split}	\end{equation}
From \eqref{esty2} and Lemma~\ref{uni-int} we conclude there exists   
a monotone, continuous function  $L:\R_+\to\R_+$, with $L(0)=0$, 
 independent of $n$, such that 
	$$\int_{s_1}^{s_2}\| A^n\|_{L^1(\Omega)}
	\le L(|s_1-s_2|)$$
for all $s_1,s_2\in[0,T]$. Conseqently,  \Ref{dt1} gives us
	$$\left|\int_{s_1}^{s_2}\left\langle\frac{d u^n}{dt},\varphi\right\rangle dt\right|	\le
	L(|s_1-s_2|)$$
for all $\varphi$ with 
${\rm supp}\ \varphi\subset(s_1,s_2)\subset[0,T]$ and 
$\|\varphi\|_{L^\infty(0,T;W^{r,2}_0)}\le 1$. 
Since 
	\begin{equation}\begin{split}
	\|u^n(s_1)-u^n(s_2)\|_{W^{-r,2}}
%	=
%	\sup\limits_{\stackrel{\psi\in V_r}{\|\psi\|_{V_r}\le1}}\left|\left\langle
%	u^n(s_1)-u^n(s_2),\psi\right\rangle\right |
	=
	\sup\limits_{\|\psi\|_{W^{r,2}_0}\le1}\left|\left\langle
	\int_{s_1}^{s_2}\frac{du^n(t)}{dt},\psi\right\rangle\right |
%	\\ \le
%	\sup\left\{\int_0^T\left|\left\langle\frac{du^n(\tau)}{d\tau},\varphi\right\rangle
%	\right|dt\ :\ \|\varphi\|_{L^\infty(0,T;V_r)}\le 1,\ {\rm supp }\ 
%	\varphi\subset(s_1,s_2)\right\}
	\end{split}\end{equation}
then
	\begin{equation}
	\label{equi}
	\sup\limits_{n\in\N}\|u^n(s_1)-u^n(s_2)\|_{W^{-r,2}}\le L(|s_1-s_2|),
	\end{equation}
namely the family of functions $\u^n:[0,T]\to W^{-r,2}(\Omega)$ is equicontinuous. 
Moreover, it is  uniformly bounded in $L^\infty(0,T;L^2(\Omega))$ and hence the sequence $\{\u^n\}$ is relatively compact 
in  ${\cal C}([0,T];W^{-r,2}(\Omega))$ and the limit $\u\in {\cal C}([0,T];W^{-r,2}(\Omega))$.  Thus there exists  a sequence 
$\{s_0^i \}_i$, 
$s_0^i \to 0^+$ as $i\to\infty$ such that 
%
%
%Since $\u\in L^{\infty}(0,T;L^2(\Omega))$, we can choose a sequence 
%$\{s_0^i \}_i$, 
%$s_0^i \to 0^+$ as $i\to\infty$. Thus $\{\u(s^i_0)\}_i$ is weakly convergent in 
%$L^2(\Omega)$. The estimate (\ref{equi}) provides 
%that the family of functions $\u^n:[0,T]\to V_r^*$ is equicontinuous. 
%Uniform bound in $L^\infty(0,T;L^2(\Omega))$ and $L^2\subset\subset W^{-r,2}(\Omega)$  
%the Arzel\`a-Ascoli theorem imply that the sequence $\{\u^n\}$ is relatively compact 
%in  $C([0,T];W^{-r,2}(\Omega))$ and $\u\in C([0,T];W^{-r,2}(\Omega))$. Consequently we obtain that
	\begin{equation}
	\u(s_0^i){{\stackrel{i\to \infty}{\longrightarrow\,}}}\u(0)\quad\mbox{in }W^{-r,2}(\Omega).
	\end{equation}
The limit above coincides with the weak limit of $\{\u(s^i_0)\}$ in $L^2(\Omega)$
what allows to claim that 
%	
%	
%Then use the fact that strong convergence in weaker space and boundetnes in 
%$L^\infty(0,T;\Omega)$ imply weak convergence in $L^2$. 
%Hence 
	\begin{equation}\label{liminfu0}
	\liminf\limits_{i\to\infty}\|\u(s_0)\|_{L^2(\Omega)}\geq\|\u_0\|_{L^2(\Omega)}.
	\end{equation} 
We obtain from \eqref{Galerkin} for any   Lebesgue point $s$  of $\u$ that 
	\begin{equation}\label{777}
	\begin{split}
	\limsup\limits_{n\to\infty}\int_{Q_s} &A^n\cdot \nabla\u^n \,dx\,dt\\
	& = 
	 \frac{1}{2} \|\u_0\|^2_{2} -
	\liminf\limits_{k\to\infty} \frac{1}{2}  \|\u^n(s)\|^2_{2}
	+\lim\limits_{n\to\infty}\int_{Q_s} f u^n\,dx\,dt\\
	& \leq
	 \frac{1}{2} \|\u_0\|^2_{2} -
	 \frac{1}{2}  \|\u(s)\|^2_{2}+\int_{Q_s} f u\,dx\,dt\\
	 &\leq
	 \liminf\limits_{i\to \infty}\left( \frac{1}{2} \|\u(s^i_0)\|^2_{2} -
	 \frac{1}{2}  \|\u(s)\|^2_{2}\right)+\int_{Q_s} f u\,dx\,dt\\
	  &{=} \lim\limits_{i\to\infty} 
	  \int_{s^i_0}^s\int_\Omega A\cdot\nabla\u \,dx\,dt \\%+\int_{Q_s} f u \,dx\,dt\\
	  &= \int_{0}^s\int_\Omega A\cdot\nabla\u \,dx\,dt
	 \end{split}
	\end{equation}
	what provides that  \eqref{44} is satisfied and  Lemma~\ref{Minty2} allows to complete the proof. 
	
\section{Approximation}
In this section we shall concentrate on the issue of density of compactly supported smooth functions with respect to the modular topology. The fundamental studies in this direction are due to Gossez for the case of classical Orlicz spaces and elliptic equations \cite{Gossez1, Gossez2}.
The similar considerations for isotropic $x-$dependent $N-$functions are due to Benkirane et al. cf.~\cite{Be2011}, see also \cite{GwMiWr2012} for anisotropic case with an application to elliptic problems. 
 Note that the main idea is analogous to \cite{GwMiWr2012}. However, Gwiazda et al. approxi\-ma\-te the truncated functions which are appropriate test functions in the con\-si\-dered elliptic equation. This is not the case of parabolic problems. Hence the presented approximation theorem is under weaker assumptions and the dependence on time is taken into account. 
Since  this  result is  essential  for proving existence of weak solutions, then we include the details 
for completeness. 

\begin{theorem}\label{main1}
If $u\in  L^\infty(0,T;L^2(\Omega))\cap L^1(0,T;W^{1,1}_0(\Omega)), \nabla u\in L_M(Q)$ 
%and $u_t\in .....$, 
then there exists a sequence $v^j \in L^\infty(0,T;{\cal C}_c^\infty(\Omega))$
satisfying 
\begin{equation}\label{zb_w_M1}
\nabla v^j\modular{M} \nabla u\ \mbox{ modularly in}\  L_M(Q) 
\ \mbox{and}\  v^j\to u \ \mbox{strongly in }\ 
L^2(Q).
\end{equation}
\end{theorem}

\begin{proof}
%The current section provides the proof of Theorem~\ref{main1}.
Already for Lipschitz domain $\Omega$ there exists a finite family of star-shaped Lipschitz domains $\{\Omega_i\}$ such that 
$$\Omega=\bigcup\limits_{i\in J}\Omega_i,$$
cf.~\cite{Novotny}. We introduce the partition of unity $\theta_i$ with $0\le\theta_i\le1,\, \theta_i\in {\cal C}^\infty_c(\Omega_i), \,
{\rm supp} \,\theta_i=\Omega_i, \sum_{i\in J}\theta_i(x)=1$ for $x\in\Omega$
%For  $u\in Z^M_0$ one easily observes that $u\in L^1(0,T;W^{1,1}_0(\Omega))$. 
%The proof of  (\ref{inc}) we split also into two parts. 
%For readability of the proof we introduce the notation $$V_M:=\{u \in L^1(0,T;W^{1,1}_0(\Omega))\cap L^\infty(0,T;L^2(\Omega))\,\,\,|\,\,\, \nabla u \in L_M(Q)\}.$$
and  define the truncation operator $T_\ell(u)$  as follows 
\begin{equation}
T_\ell(u)=\left\{
\begin{array}{rcl}
u&{\rm if}&|u|\le \ell,\\
\ell&{\rm if}& u>\ell,\\
-\ell&{\rm if}& u< -\ell.
\end{array}\right.\end{equation}
Define $Q_i:=(0,T)\times \Omega_i$. Obviously $$T_\ell(u)\in L^\infty(0,T;L^2(\Omega))\cap L^1(0,T;W^{1,1}_0(\Omega)), \nabla T_\ell u\in L_M(Q)$$ and for each 
$i\in J$
$$\theta_i\cdot T_\ell(u)\in L^\infty(Q_i)\cap L^1(0,T;W^{1,1}_0(\Omega_i))\cap L^\infty(0,T;L^2(\Omega_i)) .$$
 Introducing the truncation of $u$ was necessary to provide that 
$$\nabla T_\ell(u)\cdot \theta_i + T_\ell(u)\cdot \nabla \theta_i= \nabla (T_\ell(u)\cdot \theta_i) \in
 L_M(Q_i).$$
%hence\marginpar{zbior gwiazdzisty wzgl kuli, wiec $\vep_\lambda$ zalezy od promienia kuli}
%$\theta_i\cdot T_\ell(u)\in  Z^M_0$ and consequently it is also bounded  in $L^1(0,T;W^{1,1}_0(\Omega))$. 
%Hence using the introduced partition of unity we may now without loss of generality assume that $\Omega$ is a star-shaped domain. 
Without loss of generality assume that all $\Omega_i$ are  star-shaped domains with respect to a ball of radius R, i.e. $B(0,R)$.
% then for each $\lambda\in(0,1)$
%$$(1-\lambda)x+\lambda y\in \Omega\quad\mbox{for each }x\in \Omega, y\in B(0,R).
%$$
%Hence 
%for $\delta<R$ we may choose $\lambda=\delta/R$ and conclude that 
%$$\left(1-\frac{\delta}{R}\right)\Omega+\delta B(0,1)\subset \Omega.$$
%
%As the  set $(0,T)\times \Omega$ is also star-shaped, then
%
%In the next step  
% we will show that 
%$$(V_M\cap L^\infty(Q)) \subset Y^M_0,$$
%which together with \Ref{domk}  and \Ref{inc1} will prove \Ref{inc}.
%If $u\in V_M\cap  L^\infty(Q)$ then
%$$u\cdot \theta_i\in  L^1(0,T;W^{1,1}_0(\Omega_i))\cap L^\infty(0,T;L^2(\Omega_i)) \cap  L^\infty((0,T)\times\Omega_i))$$
%and $$\nabla u\cdot \theta_i + u\cdot \nabla \theta_i= \nabla (u\cdot \theta_i) \in
% L_M((0,T)\times\Omega_i),$$ where $\Omega_i={\rm supp}\ \theta_i$. % and $\nabla (u\cdot \theta_i) =
% Now we follow the case of star-shaped domains to complete the proof. \\
%The proof of an opposite inclusion we will split into two steps. In the first step we will assume 
% that 
%$\Omega$ is a star-shaped domain, whereas in the second step we  extend the idea for arbitrary Lipschitz domains.  \\
%\noindent
% {\it Step 1. Star-shaped domains.} 
% We are aiming to show that 
% \begin{equation}\label{inc}
% Z^M_0 \subset Y^M_0.
% \end{equation}
% The starting point is extending the function $u$ by zero outside of  $\Omega$ to the whole $\R^d$ to mollify it. 
% To provide the extension is properly defined we split (\ref{inc}) as follows
% 
% \begin{equation}\label{inc1}
% Z^M_0 \subset V_M
% \end{equation}
%  and
%  \begin{equation}\label{inc2}
%V_M\subset Y^M_0.
% \end{equation}
% 
% Inclusion (\ref{inc1}) is obvious and to prove (\ref{inc2})
 %Also  $Z^M_0 \subset \{u \in L^1(0,T;W^{1,1}_0(\Omega))\,\,\,|\,\,\, \nabla u \in L_M(Q)\}$.
%It remains  to prove that 
%\begin{equation}\label{W11}
%\{u \in L^1(0,T;W^{1,1}_0(\Omega))\,\,\,|\,\,\, \nabla u \in L_M(Q)\} \subset W^M_0.
%\end{equation}
We define for $(t,x)\in(0,T)\times\Omega$
%$$T_\ell^\delta(u)(t,x):=T_\ell(u)(t+(1-\delta/R)x)$$
%and
\begin{equation}\label{Sdelta2}
\begin{split}
 {\cal S}_\delta (\theta_iT_\ell(u))(t,x)
:=\frac{1}{\left(1-{\delta}/{R}\right)}\int_Q
  S_\delta(x-y)\theta_i T_\ell(u)\left(t,\left(1-{\delta}/{R}\right)y\right)  \,dy.
\end{split}\end{equation}
%$$T_\ell^\lambda(u)(t,x):=T_\ell(u)(\lambda(t-T/2)+T/2,\lambda (x-x_0)+x_0)$$  where $x_0$ is a vantage point of $\Omega$.
%Let $\varepsilon_\lambda=\frac{1}{2}{\rm dist}\,(Q,\partial(\lambda Q))$ where $\lambda Q:=\{(s,y)=(\lambda(t-T/2)+T/2,\lambda  (x-x_0)+x_0)\,\,|\,\,(t,x)\in Q\}$. 
%Define then\marginpar{tutaj nie definiowac tego tak, tylko tak jak w lemacie poprzednio}  $$[K_\varepsilon T_\ell^{\lambda}(u)](t,x):= K_\varepsilon*T_\ell^\lambda(u)(t,x)$$
%where $ K_\varepsilon$ is a standard mollifier,  the convolution is done w.r.t. $x$ and $t$ and $\varepsilon<\varepsilon_\lambda$.
%This approximation has the property that if  $u\in V_M$, then
%$u^{\lambda,\varepsilon} \in V_M$
%Note that for $\varepsilon<\varepsilon_\lambda$ the approximation 
%$$[K_\varepsilon T_\ell^{\lambda}(u)]\in C^\infty([0,T]\times\overline\Omega)\cap Z_0^M\cap L^1(0,T;W^{1,1}_0(\Omega)).$$
Our aim is to show that there exists a constant $\lambda>0$ such that
\begin{equation}\label{216}
\lim\limits_{l\to\infty}\lim\limits_{\delta\to 0_+}\varrho_{M,Q_i}\left( \frac{\nabla u -\nabla {\cal S}_\delta(\theta_i T_\ell(u))}{\lambda} \right)
=0.
\end{equation}
For this purpose we introduce a sequence of simple functions
$$\xi_n(t,x):=\sum_{j=1}^n \alpha_j^n \bbbone_{G_j}(t,x),\quad  \alpha_j^n \in\R,\ \bigcup_{j\in\{1,\ldots,n\}} G_j=Q$$ 
which converges to $ \nabla (\theta_i\cdot T_\ell(u)) $ modularly in $L_M(Q)$. 
Moreover, let $\lambda_0, \lambda_1, \lambda_2$ be some appropriate constants which we specify later such that the following estimate holds
\begin{equation}
\begin{split}
\varrho_{M,Q_i}&\left( \frac{\nabla u -\nabla {\cal S}_\delta(\theta_i T_\ell(u))}{\lambda} \right)\\
&\le \frac{\lambda_0}{\lambda}
\rho_{M,Q_i}\left(\frac{ {\cal S}_\delta  \nabla(\theta_iT_\ell( u))- {\cal S}_\delta \xi_n}{\lambda_0}\right)
+ \frac{\lambda_0}{\lambda}
\rho_{M,Q_i}\left(\frac{\nabla (\theta_iT_\ell( u))-\xi_n}{\lambda_0}\right)\\
&\quad+\frac{\lambda_1}{\lambda}\rho_{M,Q_i}\left(\frac{ {\cal S}_\delta \xi_n-\xi_n}{\lambda_1}\right)
+\frac{\lambda_2}{\lambda}\varrho_{M,Q_i}\left( \frac{ \nabla u -\nabla (T_\ell(u)\theta_i)}{\lambda_2}\right)\\
&=I_1+I_2+I_3+I_4.
%&\le \frac{\lambda_0(1+c)}{\lambda}
%\rho_M\left(\frac{\nabla (\theta_iT_\ell( u))- \xi_n}{\lambda_0}\right)
%+\rho_M\left(\frac{ {\cal S}_\delta \xi_n-\xi_n}{\lambda_1}\right)
%+\rho\left(\frac{\xi-\xi_n}{\lambda_0}\right)
%\right]
\end{split}\end{equation}
Consider first $I_3$. 
The existence of a sequence $\xi_n$ is provided by Lemma \ref{lem-dense}.
%Note that the sequence  $\xi_n$ does not have to be in a gradient form.  
Let 
$B_\delta:=\{y\in \Omega\ :\ |y|<\delta\}$.
Then 
\begin{equation}
\begin{split}
{\cal S}_\delta \xi_n-\xi_n
=\int_{B_\delta}  S_\delta(y)\sum\limits_{j=1}^n\left(\alpha_j^n\bbbone_{G_j}
(t, (1-\delta/R)(x-y))-
\alpha_j^n\bbbone_{G_j}(t, x)\right) \,dy
\end{split}\end{equation}
and the  Jensen inequality and Fubini theorem yield
\begin{equation}\label{proste}\begin{split}
&\rho_{M,Q_i}\left(\frac{ {\cal S}_\delta \xi_n(t,x)-\xi_n}{\lambda_1}\right)\\
&=\int_Q M(x,\frac{1}{\lambda_1}\int_{B_1}  S(y)\sum\limits_{j=1}^n(\alpha_j^n\bbbone_{G_j}
(t, (1-\delta/R)(x-\delta y))\\&\qquad-
\alpha_j^n\bbbone_{G_j}(t, x)) \,dy) \,dt\,dx\\
&\le \int _{B_1}  S(y)( \int_Q  M(x,\frac{1}{\lambda_1} \sum\limits_{j=1}^n\alpha_j^n(\bbbone_{G_j}
(t, (1-\delta/R)(x-\delta y))\\&\qquad-
\bbbone_{G_j}(t, x))) \,dt\,dx) \,dy.\\
%&\le \int _{B_1}  S(y)( \int_Q  M(x, \frac{1}{\lambda_1} \sum\limits_{j=1}^n\alpha_j^n(\bbbone_{G_j}
%(t, (1-\delta/R)(x-\delta y))\\&\qquad-
%\bbbone_{G_j}(t, x)) \,dt\,dx) \,dy.
\end{split}\end{equation}
Note that  $\{\frac{1}{\lambda_1} \sum\limits_{j=1}^n\alpha_j^n\left(\bbbone_{G_j}
(t, (1-\delta/R)(x-\delta y))-
\bbbone_{G_j}(t, x)\right) \,dt\,dx) \}_{\delta>0}$ converges a.e. in $Q$ to zero as $\delta\to0_+$ and 
\begin{equation}\begin{split}\label{estimate2}
M(x, \frac{1}{\lambda_1} \sum\limits_{j=1}^n\alpha_j^n\left(\bbbone_{G_j}
(t, (1-\delta/R)(x-\delta y))-
\bbbone_{G_j}(t, x)\right)\\
\le \sup\limits_{|z|=1}M(x,\frac{1}{\lambda_1} \sum\limits_{j=1}^n\alpha_j^n z).
\end{split}\end{equation}
Assumption \eqref{int} provides that  the right-hand side of \eqref{estimate2} is integrable, hence the Lebesgue dominated convergence theorem allows to conclude that $I_3$ vanishes as $\delta\to0_+$. 
%\begin{equation}\label{zb-prostych}
%\lim\limits_{\delta\to0_+}\rho_M\left(\frac{ {\cal S}_\delta \xi_n(t,x)-\xi_n}{\lambda_1}\right)=0.
%\end{equation}
%According to ~Lemma~\ref{lem-dense} there exists $\lambda_0>0$ such that  
%\begin{equation}\label{dense}
%\lim\limits_{n\to\infty}\rho_M\left(\frac{\xi_n-\nabla (\theta_i T_\ell( u))}{\lambda_0}\right).
%\end{equation}
Lemma~\ref{modular-topology}  allows to estimate $I_1$ on each $\Omega_i$ as follows
\begin{equation}
I_1=\frac{\lambda_0}{\lambda}\rho_{M,Q_i}\left(\frac{ {\cal S}_\delta  (\nabla(\theta_iT_\ell( u))- \xi_n)}{\lambda_0}\right)
\le c \rho_{M,Q_i}\left(\frac{ \nabla(\theta_iT_\ell( u))- \xi_n}{\lambda_0}\right)
\end{equation}
and hence by Lemma~\ref{lem-dense} there exists a constant $\lambda_0$ such that 
$$\lim\limits_{n\to\infty} (I_1+I_2)=0.$$
%By \eqref{zb-prostych} and \eqref{dense}, passing first with $\delta\to0_+$ and then with 
%$n\to \infty$  we conclude that $\lim\limits_{\delta\to0_+}\rho_M\left(\frac{ {\cal S}_\delta
%\nabla (\theta_i T_\ell( u))-\nabla (\theta_iT_\ell( u))}{\lambda}\right)=0.$
%%Then we pass with $\lambda\to1$ and obtain that $T_\ell^{\lambda}(u)\xrightarrow{\lambda\to1}T_\ell( u)$ in $L^1(0,T;W^{1,1}_0(\Omega))$ and 
%%$\nabla T_\ell(u)^{\lambda}\xrightarrow{\lambda\to1} \nabla T_\ell(u)$ modularly in $L_M(Q)$.
Moreover, as $\ell\to\infty$ we observe the following convergence
$$T_\ell(u)\to u\ \ \  {\rm strongly \,in}\ \ \   L^1(0,T;W^{1,1}_0(\Omega))$$
and hence also, at least for a subsequence, almost everywhere.
To find a uniform estimate we observe  that 
 $M(x,\nabla T_\ell(u(t,x)))\le M(x,\nabla u(t,x))$ a.e. in Q. Indeed,  $T_\ell(u)$ and $u$ coincide for $|u|\le \ell $ and on the remaining two sets, where $T_\ell (u)$
 is equal to $\ell$ or $-\ell$ we have that
 $ T_\ell(u)\in L^1(0,T;W^{1,1}_0(\Omega))$, then $\nabla T_\ell(u)$  is almost everywhere equal to zero.
  Consequently $M(x,\nabla T_\ell(u(t,x)))$ is uniformly 
integrable, which combined with pointwise convergence provides
$$\nabla T_\ell(u)\to \nabla u\ \ \  {\rm modularly\  in}\ \ \  L_M(Q)$$
as $\ell\to\infty$, hence there exists a constant $\lambda_2$ such that $\lim_{\ell\to\infty} I_4=0$.
Finally, choosing $\lambda>\max\{3 \lambda_0, 3\lambda_1, 3\lambda_2\}$, passing first with $\delta\to0_+$, then 
$n\to\infty$ and $\ell\to\infty$ we arrive to \eqref{216}.

The strong convergence in $L^2$ is straighforward, since an $N-$function $M(x,a)=|a|^2$ satisfies 
$\Delta_2-$condition and the strong and modular convergence coincide. 

%In the second step we shall show that for this approximation the conditions on the time derivative are valid.  Let us 
%write $\frac{\partial u}{\partial t} =\div v_A+v_f$ with $v_A\in L_{M^*}(Q)$ and $v_f\in L^\infty(Q)$. We will show there exists 
%$\lambda_2$ and sequences $v_A^j, v_f^j$ such that 
%\begin{equation}
%\lim\limits_{j\to\infty}\rho_{M_*}\left(\frac{v_A^j-v_A}{\lambda_2}\right)=0
%\end{equation}
%and 
%\begin{equation}
%v_f^j\to v_f\quad \mbox{strongly in}\ L^2(Q).
%\end{equation}
%Observe that 
%\begin{equation}
%\begin{split}
%\frac{\partial S_\delta(\theta_i\beta(u))}{\partial t}&=
%S_\delta\theta_i\frac{\partial\beta(u)}{\partial t}
%=
%S_\delta\theta_i(\div v_A+v_f)\\&=S_\delta \div(v_A\theta_i)
%-S_\delta\nabla \theta_iv_A+S_\delta\theta_iv_f\\&=\div(S_\delta(\theta_iv_A))
%+S_\delta\theta_iv_f-S_\delta\nabla \theta_iv_A
%\end{split}\end{equation}
%\qed\\
%\end{proof}
%\begin{corollary}\label{Cor}
%Let an $N-$function satisfy ... and let the assumptions of Lemma~\ref{closures} be satisfied. Then 
%for any $\varphi\in ...$ there exists a sequence  $\{\varphi^j\}_{j=1}^{\infty} \subset
%C^\infty_c((-\infty,T)\times\Omega)$ such that there exists $\lambda>0$ and $\rho_M\left(\frac{1}{\lambda}(\varphi^j-\varphi)\right)\to0$ as $j\to\infty$.
%\end{corollary}


%Arbitrary Lipschitz domains - Step 2
%\begin{itemize}
%\item
%$T_n(u)$ -- truncation of the function $u$
%\item
%If $$u\in W^{1,1}_0(\Omega), \nabla u \in L_M (\Omega)$$ then
%$$T_n(u)\in  W^{1,1}_0(\Omega), \nabla T_n(u)\in L_M (\Omega).$$
%\item
%$$T_n(u)\to u\ \ \  {\rm in}\ \ \  W^{1,1}_0(\Omega) $$
%\item
%$$\nabla T_n(u)\to \nabla u\ \ \  {\rm in}\ \ \  L_M(\Omega) $$

%\end{itemize}
%\noindent
%{\bf Remark} Note that in the classical case $M=M(|\cdot |)$ (i.e. the modular does not depend on the direction)  the proof is simpler.  The problem which arises here is the  lack of proper Poincar\'e inequality, cf. \cite{Cianchi}, which in the isotropic case has the following form  
%$$\int_QM(|\nabla u(t,x)|)dxdt\ge c \int_Q M(|u|) dxdt. $$

% Integrating by parts
%{\bf  Classical case}
%$H$-- Hilbert space, $X$-- Banach space, 
%$X\hookrightarrow H\simeq
%H^\ast\hookrightarrow X^\ast$, $u\in L^p(0,T;X), \frac{du}{dt}\in
%L^{p'}(0,T;X^\ast)$. 
%Then for all $s_0,s_1\in (0,T)$
%$$\int_{s_0}^{s_1}\left\langle \frac{du(t)}{dt},u(t)\right\rangle_Xdt=
%\frac{1}{2}\|u(s_1)\|_H^2-\frac{1}{2}\|u(s_0)\|^2_H.$$ 

%We know
%that \begin{equation*} \begin{split} u&\in L^\infty(0,T;L^2(\Omega)),\\
%Du&\in L_M(Q), \,\,\, u_t,\chi\in L_{M^\ast}(Q)
%\end{split}\end{equation*}
%{\bf Problem:}
% $$L_M(Q)\neq L_M(0,T;L_M(\Omega))$$

\end{proof}




%\begin{lemma}\label{kapelusz}
%Let $\varphi\in{\cal C}_c^\infty(\R^{d+1})$ with $0\le\varphi\le 1,\ \varphi(t,x)=1$ for $|(t,x)|\le1$ and 
%$\varphi(t,x)=0$ for  $|(t,x)|\ge 2$. Let $\varphi_r(t,x)=\varphi(t/r, x/r)$ and let $z_r=\varphi_r z$. If 
%$z\in{\cal L}_M(Q)$, then \marginpar{zamiescic dowod}
%\begin{equation}
%\lim\limits_{r\to\infty}\rho_M(z_r-z)=0
%\end{equation}
%
%\end{lemma}


%\begin{lemma}\label{ciaglosc-1}
%Let  $K\in
%C^\infty_c(\R)$  and
%$\int_{\R} K(\tau)d\tau=1, K(t)=K(-t)$. We define
%$ { K}_\varepsilon(t):=1/\vep K(t/\vep).$ Consider the family of operators
%\begin{equation}
% {\cal K}_\varepsilon z(t,x)=K_\vep*[(K_\vep*z)\cdot\bbbone_{(s_0,s)}]
% %\int_0^T  K_\varepsilon(t-s)z(s,x)  \,dw.
%\end{equation}
%where $\vep<\frac{1}{2}\min\{s_0,T-s\}$ and the convolution is done with respect to  time.
%Moreover, let an $N-$function satisfy condition \eqref{log-t} and:
%\begin{itemize}
%\item there exists constant $c>0$ such that $\inf_{(t,x)\in Q, |\xi|=1}M(x,\xi)\ge c$
%\item %if $G$ is bounded measurable set, then 
%$\int_Q M(x,\xi)\,dx\,dt\le \infty$ for every $\xi\in\R^d$
%\item $M^*(t,x,1)\le C$ a.e.
%\end{itemize}
%Then there exists a constant $c>0$ (independent of $\vep$) such that 
%\begin{equation}\label{cont}
%\int_Q M(x, {\cal K}_\varepsilon z(t,x))\;dx\;dt\le c\int_G M(x,z(t,x))\;dx\;dt
%\end{equation}
%for every $z\in L_M(Q)$.
%\end{lemma}
%
%\begin{proof}
%Hence
%\begin{equation}
%\int_G M(y,[K_\varepsilon z](y))dy=\sum\limits_{k\in I}\int_{5B_{i_k}}M(y,[K_\varepsilon z](y))dy.
%\end{equation}
%Define 
%\begin{equation}
%m_k^\vep(\xi):=\inf_{y\in 5B_{i_k}\cap G}M(y,\xi)\le \inf_{y\in B_{i_k}\cap G}M(y,\xi)
%\end{equation}
%and 
%\begin{equation}
%\alpha_k(y,\vep):=\frac{M(y,[K_\varepsilon z](y))}{m_k^\vep([K_\varepsilon z](y))}.
%\end{equation}
%Then obviously
%\begin{equation}
%\int_G M(y,[K_\varepsilon z](y))dy=\sum\limits_{k\in I}\int_{5B_{i_k}}\alpha_i(y,\vep)m_k^\vep([K_\varepsilon z](y)) dy.
%\end{equation}
%We shall now concentrate on the uniform estimates of $\alpha_i(y,\vep)$ for sufficiently small $\vep$ and 
%$x\in B_{i_k}$.  Without loss of generality one can assume that $\|z\|_M\le 1$. By H\"older inequality \eqref{hoelder} we obtain for $\vep<R/2$
%\begin{equation}\label{22}\begin{split}
%|[K_\varepsilon z](y)|&\le\left| \frac{1}{\vep^n}\left(1-\frac{\vep}{R}\right)^{-1}\sup_{B(0,1)}|K(y)|\int_G\bbbone_{B(0,\vep)}(y)z((1-\frac{\vep}{R})y) \,dy\right|\\
%&\le \frac{2}{\vep^n}\sup_{B(0,1)}|K(y)|\|z\|_{1}\le \frac{c}{\vep^n}\|z\|_M
%\le \frac{c}{\vep^n}
%\end{split}\end{equation}
%Let now $y_k$ be the point where the infimum of $M(y,\xi)$ is obtained in the set $5B_{i_k}$. Then by log-Hoelder
%regularity we have
%\begin{equation}
%\alpha_k(y,\vep)=\frac{M(y,[K_\varepsilon z](y))}{M(y_k,[K_\varepsilon z](y))}\le |[K_\varepsilon z](y) |^\frac{H}{\ln\frac{1}{|y-y_k|}}
%\end{equation}
%And as $y\in B_{i_k}$ and $y_k\in 5B_{i_k}$ then   $|y-y_k|\le \frac{5}{2}\vep$. Hence for sufficiently small $\vep$, e.g. $\vep<\frac{1}{5}$ we have
%$$|[K_\varepsilon z](y) |^\frac{H}{\ln\frac{1}{|y-y_k|}}\le |[K_\varepsilon z](y) |^\frac{H}{\ln\frac{2}{5\vep}}.
%$$ 
%Further we use  \eqref{22} to estimate as follows
%\begin{equation}
% |[K_\varepsilon z](y) |^\frac{H}{\ln\frac{2}{5\vep}}\le(c\vep^{-n})^\frac{H}{\ln\frac{2}{5\vep}}
% \le c^\frac{H}{\ln2} \cdot \frac{2}{5}\left(e^{\ln\frac{2}{5\vep}}\right)^\frac{nH}{\ln\frac{2}{5\vep}}
% \le \frac{2}{5} c^\frac{H}{\ln2}e^{nH}:=C
%\end{equation}
%Using this uniform estimate we have
%\begin{equation}
%\int_G M(y,[K_\varepsilon z](y))dy\le C\sum\limits_{k\in I}\int_{5B_{i_k}}m_k^\vep([K_\varepsilon z](y)) dy.
%\end{equation}
%
%\end{proof}



\begin{lemma}\label{modular-topology} Let an $N-$function satisfy condition (M),
  $S$ and $S_\delta$ be given by \eqref{S}
%\in
%C^\infty_c(\R^{d})$  and
%$\int_{\R^{d}} S(y)\;dy\;d\tau=1, S(x)=S(-x)$. We define
%$ S_\delta(x):=1/\delta^{d}S(x/\delta).$ 
and  assume that $\Omega$ is a star-shaped domain  with respect to a ball centered at the origin $B(0,R)$ for some $R>0$.
We define  the family of operators
\begin{equation}\label{Sdelta}
 {\cal S}_\delta z(t,x):=\left(1-{\delta}/{R}\right)^{-1}\int_\Omega
  S_\delta(x-y)z\left(t,\left(1-{\delta}/{R}\right)y\right)  \,dy.
\end{equation}
 % and:
%\begin{itemize}
%\item there exists constant $c>0$ such that $\inf_{(t,x)\in Q, |\xi|=1}M(x,\xi)\ge c$
%\item  $\int_Q M(x,\xi)\le \infty$ for every $\xi\in\R^d$
%\item $M^*(t,x,1)\le C$ a.e.
%\end{itemize}
Then there exists a constant $c>0$ (independent of $\delta$) such that 
\begin{equation}\label{cont}
\int_Q M(x, {\cal S}_\delta z(t,x))\,dx\,dt\le c\int_Q M(x,z(t,x))\,dx\,dt
\end{equation}
holds for every $z\in L_M(Q)\cap L^\infty(0,T;L^1(\Omega))$.
\end{lemma}

\begin{proof}
Since $\Omega$ is
a star-shaped domain with respect to  $B(0,R)$, then for each $\lambda\in(0,1)$
$$(1-\lambda)x+\lambda y\in \Omega\quad\mbox{for each }x\in \Omega, y\in B(0,R).
$$
Hence 
for $\delta<R$ we may choose $\lambda=\delta/R$ and conclude that 
$$\left(1-\frac{\delta}{R}\right)\Omega+\delta B(0,1)\subset \Omega.$$

%Define the family of operators
%\begin{equation}
% {\cal S}_\delta z(t,x)=\left(1-\frac{\delta}{R}\right)^{-1}\int_Q  S_\delta(x-w)z\left(t,\left(1-\frac{\delta}{R}\right)w\right) \,dw.
%\end{equation}
Let ${\cal S}_\delta z(t,x)$ be defined by~\eqref{Sdelta}. 
Since $\overline{\left(1-\frac{\delta}{R}\right)\Omega+\delta B(0,1)}\subset \Omega,$ then %for each $t\in(0,T)$ 
it holds
$ {\cal S}_\delta z\in L^\infty(0,T;{\cal C}^\infty_c(\Omega))$. For every $\delta>0$ there exists $N=N(\delta)$ such that  a family of closed
cubes $\{D_{\delta,k}\}_{k=1}^N$ with disjoint interiors and the length of an edge equal to $\delta$ covers $\Omega$, i.e.
$\Omega\subset \bigcup_{k=1}^ND_{\delta,k}$. 
 Hence %for a.a. $t\in(0,T)$
%
% balls ${\cal F}=\{B_i\}_{i=1}^{N}$ that cover the set $G$ 
%such that for each $i\in\{1,\ldots,N(\vep)\}$
%$${\rm diam}\,B_i\le \vep.$$
%By the Vitali covering lemma there exists a disjoint subfamily ${\cal F}'$ of balls in ${\cal F}$, i.e. 
%${\cal F}'=\{B_{i_k}\}_{k\in I}$ such that
%$$G\subset\bigcup_{k\in I} 5B_{i_k}, \qquad |G\setminus \bigcup_{k\in I} 5B_{i_k}|=0 $$
%where $5B_i$ is the ball centered the same as $B_i$ and with a radius five times the radius of $B_i$.
%%$$G=\bigcup_{i=1}^{N}\gamma^\vep_i.$$
%and hence
%$$\int_\Omega \xi
\begin{equation}
\int_0^T\int_\Omega M(x,{\cal S}_\delta z(t,x))\,dx=%\sum\limits_{i=1}^{[T/\delta]}
\sum\limits_{k=1}^{N}
\int_0^T\int_{D_{\delta,k}\cap\Omega}M(x,{\cal S}_\delta z(t,x))\,dx\,dt.
\end{equation}
For each $\delta, k$ by $G_{\delta,k}$ we shall mean a cube with an edge of the length $2\delta$ and centered the same as the corresponding $D_{\delta,k}$. 
Note that if $x\in D_{\delta,k}$, then there exist $2^d$ cubes 
$G_{\delta,k}$ such that $x\in G_{\delta,k}$. 
%Then divide the interval $[0,T]$ for the 
%subintervals of the length $\delta$, which we will call $I_{i,\delta}$. Moreover by $J_{i,\delta}$ we shall mean the intervals of the length $2\delta$, namely $((i-3/2)\delta,(i+1/2)\delta)$. 
 Define 
\begin{equation}
m_{k}^\delta(\xi):=\inf_{(t,x)\in ((0,T)\times G_{\delta,k})\cap Q}M(x,\xi)\le 
\inf_{(t,x)\in ((0,T)\times D_{\delta,k})\cap Q}M(x,\xi)
\end{equation}
and 
\begin{equation}
\alpha_{k}(t,x,\delta):=\frac{M(x,{\cal S}_\delta z(t,x))}{m_{k}^\delta({\cal S}_\delta z(t,x))}.
\end{equation}
Then 
\begin{equation}
\int_0^T\int_\Omega M(x,{\cal S}_\delta z(t,x))\,dx\,dt=\sum\limits_{k= 1}^N
\int_0^T\int_{D_{\delta,k}\cap\Omega}\alpha_k(t,x,\delta)m_k^\delta
({\cal S}_\delta z(t,x)) \,dx\,dt.
\end{equation}
We are aiming to estimate the term 
 $\alpha_{k}(t,x,\delta)$ and the main tool here will be the regularity with respect to $x$, which is assumed on M, namely
 condition \eqref{log}. 
 For this purpose let now $(t_k,x_k)$ be the point where the infimum of $M(x,\xi)$ is obtained in the set $(0,T)\times G_{\delta,k}$. Then 
\begin{equation}
\alpha_{k}(t,x,\delta)=\frac{M(x,{\cal S}_\delta z(t,x))}{M(x_k,{\cal S}_\delta z(t,x))}\le |{\cal S}_\delta z(t,x) |^\frac{H}{\ln\frac{1}{|x-x_k|}}.
\end{equation}
%   for sufficiently small $\delta$ and 
%$(t,x)\in (0,T)\times D_{\delta,k}$.  
Without loss of generality one can assume that 
$\|z\|_{L^\infty(0,T;L^1(\Omega))}\le 1$. By H\"older inequality \eqref{hoelder} we obtain for $\delta<R$
\begin{equation}\label{22}\begin{split}
|{\cal S}_\delta z(t,x)|&\le\left| \frac{1}{\delta^{d}}\left(1-\frac{\delta}{R}\right)^{-1}\sup_{B(0,1)}|S(y)|\int_\Omega\bbbone_{B(0,\delta)}(y)z(t,(1-\frac{\delta}{R})y) \,dy\right|\\
&\le \frac{2}{\delta^{d}}\sup_{B(0,1)}|S(y)|\|z\|_{L^\infty(0,T;L^1(\Omega))}%\le \frac{c}{\delta^{d}}\|z\|_M
\le \frac{c}{\delta^{d}}.
\end{split}\end{equation}
Since $x\in D_{\delta,k}$ and $x_k\in G_{\delta,k}$ then   $|x-x_k|\le \delta\sqrt{d}$  and for sufficiently small $\delta$, e.g. $\delta<\frac{1}{2\sqrt{d}}$ 
%we have
%$$|{\cal S}_\delta z(t,x) |^\frac{H}{\ln\frac{1}{|x-x_k|}}\le |{\cal S}_\delta z(t,x) |^\frac{H}{\ln\frac{1}{\delta\sqrt{d}}}.
%$$ 
with  use of  \eqref{22} we obtain %to estimate as follows for $\delta<\frac{1}{2^d\sqrt{d}}$
\begin{equation}\begin{split}
 |{\cal S}_\delta z(t,x) |^\frac{H}{\ln\frac{1}{\delta\sqrt{d}}}&\le(c\delta^{-d})^\frac{H}{\ln\frac{1}
 {\delta\sqrt{d}}}
 \le c^\frac{H}{\ln2} \cdot d^\frac{dH}{\ln 4}\left(e^{\ln \delta\sqrt{d}}\right)^\frac{dH}{\ln\delta\sqrt{d}}
 \le d^\frac{dH}{\ln 4} c^\frac{H}{\ln2}e^{dH}:=C.
\end{split}\end{equation}
Consequently
\begin{equation}\label{osza}
|\alpha_{k}(t,x,\delta)|\le C.
\end{equation}
Define $\tilde M(x,\xi):=\max_{k}m_{k}^\delta(\xi)$ where the maximum is taken with respect to all the sets 
$(0,T)\times G_{\delta,k}$. Obviously  $\tilde M(x,\xi)\le M(x,\xi)$ for all $(t,x)\in Q.$
Using the uniform estimate \eqref{osza} and the Jensen inequality we have 
\begin{equation}\begin{split}
\int_Q &M(x,{\cal S}_\delta z(t,x))dxdy\le C\sum\limits_{k=1}^N
\int_0^T\int_{D_{\delta,k}}m_{k}^\delta({\cal S}_\delta z(t,x)) 
\,dx\,dt\\&
\le C\sum\limits_{k=1}^N\int_{B(0,\delta)}|S_\delta(y)|\,dy
\int_0^T\int_{(1-\frac{\delta}{R})G_{\delta,k}}
m_{k}^\delta( z(t,x))\,dx\,dt\\&
\le 2^{d}C\int_{Q}\tilde M(x,z(t,x)) \,dx\,dt
\le 2^{d}C\int_{Q} M(x,z(t,x)) \,dx\,dt
\end{split}\end{equation}
which completes the proof. 
\end{proof}
\appendix
\section{Auxilary facts}

%In this short section we collect some facts about ${\cal N}$--functions and Orlicz spaces, which are used in the proof of the main theorem. All the proofs to the collected lemmas and propositions can be found e.g in~\cite{gwiazda-swierczewska}.  
\begin{lemma}\label{lem-dense}
Let ${\mathbb S}$ be the set of all simple, integrable functions on $Q$ and let \eqref{int} hold. Then ${\mathbb S}$ is dense with respect to the modular topology in $L_M(Q)$. 
\end{lemma}
For the proof in isotropic case see \cite[Theorem 7.6]{Musielak}. The anisotropic case follows exactly the same lines.

Below we formulate some facts concerning convergence in generalized Musielak-Orlicz spaces.
For the proofs of these lemmas and propositions see~\cite{GwSw2008}. 

\begin{lemma}\label{modular-conv}
Let $z^j:Q\to\R^d$ be a measurable sequence. Then
$z^j\modular{M} z$ in $L_M(Q)$ modularly if and only if
$z^j\to z$ in measure and there exist some $\lambda>0$ such that
the sequence $\{M(x,\lambda z^j)\}$ is uniformly integrable in $L^1(Q)$,
i.e.,
$$\lim\limits_{R\to\infty}\left(\sup\limits_{j\in\N}\int_{\{(t,x):|M(x,\lambda z^j)|\ge
R\}}M(x,\lambda z^j)dxdt\right)=0.$$
%there exist a convergent sequence $\{f^j\}$ in $L^1(\Omega)$ and
%some $k>0$ such that $M(\cdot,kz^j)\le f^j$ a.e. in $\Omega$.
\end{lemma}
\begin{lemma}\label{uni-int}
Let $M$ be an  ${ N}$--function and for all $j\in\N$ let $\int_Q
M(x,z^j)\,dx\,dt\le c$.
% Moreover $\lim\limits_{|\xi|\to\infty}\inf\limits_{x\in\Omega}
%\frac{M(x,\xi)}{|\xi|}=\infty$.
Then the sequence $\{z^j\}$ is
uniformly integrable in $L^1(Q)$.
\end{lemma}
\begin{proposition}\label{product}
Let $M$ be an  ${ N}$--function and $M^\ast$ its complementary
function. Suppose that the sequences $\psi^j:Q\to\R^d$ and
$\phi^j:Q\to\R^d$ are uniformly bounded in $L_M(Q)$ and
$L_{M^\ast}(Q)$ respectively. Moreover $\psi^j\modular{M}\psi$
modularly in $L_M(Q)$ and $\phi^j\modular{M^\ast}\phi$  modularly
in $L_{M^\ast}(Q)$. Then $\psi^j\cdot\phi^j\to\psi\cdot\phi$
strongly in $L^1(Q)$.
\end{proposition}
\begin{proposition}\label{ae}
Let $ K^j$ be a standard mollifier, i.e., $ K\in
C^\infty(\R),$ $ K$ has a compact support and
$\int_\R K(\tau)d\tau=1,  K(t)= K(-t)$. We define
$ K^j(t)=j K(jt).$ Moreover let $\ast$ denote a
convolution in the variable $t$. Then for any  function
$\psi:Q\to\R^d$ such that $\psi\in L^1(Q)$ it holds
$$(\vrho^j\ast\psi)(t,x)\to\psi(t,x)\quad\mbox{in measure}.
$$
\end{proposition}


%However the situation with $t-$dependent modular in correspondence with convolution in $t$ requires a different proof. The simple argumentation by  convexity of $M$ used in \cite{gwiazda-swierczewska} is not sufficient. Therefore the proposition is proved here. 
\begin{proposition}\label{sup}
Let $ K^j$ be defined as in Proposition \ref{ae}. Given an
 ${ N}$--function $M$ and a   function $\psi:Q\to\R^d$ such that
$\psi\in{\mathcal L}_M(Q)$,
%, namely the modular $\rho_{M,Q}(\psi)$ is finite,  
the sequence $\{M(\vrho^j\ast\psi)\}$
is uniformly integrable.
\end{proposition}

%\begin{proof} We  start with an abstract fact concerning the uniform\marginpar{na razie wkopiowany nalezy go poprawic}
%integrability. Namely, the following two conditions are equivalent
%for any measurable sequence $\{z^j\}$
%\begin{enumerate}
%\item[(a)] $\forall \varepsilon>0\quad\exists \delta>0:\quad
%\sup\limits_{j\in\N}
%\sup\limits_{|A|\le\delta}\int_A|z^j(x)|dxdt\le\varepsilon$,
%$\lim\limits_{R\to\infty}\sup\limits_{j\in\N}
%\int_{\{x:|z^j(x)|\ge R\}}z^j(x)dx=0$
%\item[(b)] $\forall \varepsilon>0\quad\exists \delta>0:\quad
%\sup\limits_{j\in\N}
%\int_Q\left||z^j(x)|-\frac{1}{\sqrt{\delta}}\right|_+
%dxdt\le\varepsilon$,
%$\lim\limits_{R\to\infty}\sup\limits_{j\in\N}
%\int_{\{x:|z^j(x)|\ge R\}}(z^j(x)-R)dx=0$
%\end{enumerate}
% where we use the notation
%$$|\xi|_+=\max\{0,\xi\}.$$
%The implication $(a)\Rightarrow(b)$ is obvious. To show that also
%$(b)\Rightarrow(a)$ holds let us %choose
%$R=\frac{1}{\sqrt{|A_k|}},$ where $A_k$ is an arbitrary set and
%estimate
%\begin{equation*}\begin{split}
%\sup\limits_{j\in\N} \sup\limits_{|A|\le\delta} \int_{A}|z^j| dx
%dt&\le
%\sup\limits_{j\in\N}
%\sup\limits_{|A|\le\delta}
%|A|\cdot\frac{1}{\sqrt{\delta}} +
%\sup\limits_{j\in\N}\int_{Q}\left||z^j|-
%\frac{1}{\sqrt{\delta}}\right|_+
%dxdt\\
%&\le \sqrt{\delta}+ \sup\limits_{j\in\N}\int_{Q}\left||z^j|-
%\frac{1}{\sqrt{\delta}}\right|_+ dx dt.
%\end{split}\end{equation*}
%Notice that since $M$ is a convex function, then  the following
% inequality holds for
%all $\delta>0$
%\begin{equation}\label{R}
%\int_Q\left|M(x,\psi)-\frac{1}{\sqrt{\delta}}\right|_+dxdt\ge
%\int_Q\left|M(x,\vrho^j\ast\psi)-\frac{1}{\sqrt{\delta}}\right|_+dxdt.
% %$M(x,\psi)\in L^1(\Omega)$
%\end{equation}
% Finally, since $\psi\in{\mathcal L}(Q)$,
%then also $\int_Q|M(x,\psi)-\frac{1}{\sqrt{\delta}}|_+ dxdt$
%is finite and hence
%the integral on the right-hand side is finite.
%taking supremum over $j\in\N$ in \Ref{R} we prove the
%assertion.\end{proof}


The next  lemma  is the main tool for showing that  the limits of approximate sequences are in the graph ${\cal A}$ provided that the graph is maximal monotone. This lemma in such a form was formulated in 
\cite{BuGwMaRaSw2012}, see also \cite{Sw2013}. 
\begin{lemma}\label{Minty2}
Let $\mathcal{A}$ be maximal monotone $M$-graph.
 Assume that there are
sequences $\{A^n\}_{n=1}^\infty$ and
$\{\nabla u^n\}_{n=1}^{\infty}$ defined on $Q$ such that %for some $G'\subset G$
the following conditions hold:
\begin{align}\label{11}
(\nabla u^n(t,x),A^n(t,x))& \in \mathcal{A}(t,x) &&\textrm{ a.e. in } Q,\\\label{22}
\nabla u^n &\weakstar \nabla u &&\textrm{ weakly}^*  \textrm{in } L_M(Q),\\
A^n &\weakstar A &&\textrm{ weakly}^*  \textrm{in }  L_{M^*}(Q),\label{33}\\\label{44}
\limsup_{n\to \infty} \int_{Q}A^n \cdot \nabla u^n \,d x\,dt &\le \int_{Q} A \cdot \nabla u \, d x\,dt.
\end{align}
Then %for almost all $(t,x)\in Q'$ 
\begin{equation*}
(\nabla u(t,x),A(t,x))\in \mathcal{A}(t,x)\quad \textrm{ a.e. in } Q,\label{Minty2-2}
\end{equation*}
\end{lemma}


%\section{Selections}\label{selections}
%Let $\mA$ be a maximal monotone graph satisfying {\bf (A1)}--{\bf (A4)}. Consider mapping
%$\tilde A:\R^{d}\to\R^{d}$   assigning to each  $\bB\in\R^{d}$
%exactly one value $\tilde A(\bB)\in \R^{d}$ so that $(\bB,\tilde A(\bB))\in\mA$.
%Such a mapping $\tilde A$ is called selection of $\mA$. Obviously, each such a selection $\tilde A$ defined
%on $\R^{d\times d}_{\rm sym}$, is monotone and due to {\bf (A4)} satisfies
%\begin{enumerate}
%\item[{\bf (A4$^*$)}]  $\quad\tilde A(\nabla u)\cdot\nabla u\ge c_*(\psi(\nabla u)+\psi^*(\tilde A(\nabla u))-d_* $ for all $\nabla u\in
% \R^{d}.$
%\end{enumerate}
%It is also not difficult to observe  (see also \cite{AlAm}) that the following condition replaces {\bf (A3)}
%\begin{enumerate}
%\item[{\bf (A3$^*$)}] For $(\nabla u,\bS)\in\R^{d}$:
%$${\rm if}\  (\bS-\tilde A(\bB),\nabla u-\bB)\ge0 \textrm{ for all } \bB\in\R^{d},
%\textrm{ then } (\nabla u,\bS)\in\mA.$$
%\end{enumerate}
%In general, $\tilde A$ does not have to be a Borel function. On the other hand, there is a selection
%$\tilde A$ that is a Borel function, see e.g. \cite{AubinFrankowska}, and only such  a selection is considered in the remaining parts of the paper.
\noindent
Finally we summarize some properties of selections. 
\begin{lemma}\label{LS*}
Let $\mathcal{A}(t,x)$ be maximal monotone $M$-graph satisfying {(A1)}--{ (A5)} with measurable selection
$\tilde A:Q\times \mathbb{R}^{d} \to \mathbb{R}^{d}$. Then $\tilde A$
% such that
%Then, according to Francfort et al.
%\cite{francfort-m-t}, there exists a (possibly discontinuous)
%function $\tilde A\colon Q \times \Rd \to \Rd$ (a selection)
%such that
 %for all $\xi \in \Rd$ we have $(\xi, \tilde A(t,x,\xi))
%\in \mA(t,x)$ and $\tilde A$
satisfies the following conditions:
\begin{enumerate}
\item [{\it (a1)}]$\mathrm{Dom}\, \tilde A(t,x,\cdot) = \Rd$ a.e. in $Q$;
\item [{\it (a2)}]$\tilde A$ is  monotone, i.e. for every $\xi_1$,
$\xi_2 \in \Rd$ and a.a. $(t,x)\in Q$
\begin{equation} \label{monot}
(\tilde A(t,x,\xi_1) - \tilde A(t,x,\xi_2))\cdot( \xi_1 - \xi_2) \ge 0;
\end{equation}
\item [{\it (a3)}] There are non-negative $k\in
L^1(Q)$,   $c_*>0$ and $N$-function  $M$ such that for all $\nabla u\in\Rd$ the function
 $\tilde A$ satisfies
\begin{equation} \label{growthS*}
\tilde A \cdot \nabla u \geq -k(t,x) +c_*(M(x,\nabla u) +M^*(x,\tilde A))
\end{equation}
%\item [{\it a4.}]for all $\xi \in \Rd$ the function $\tilde A$
%satisfies  the following growth and coercivity conditions a.e. in
%$Q$
%\begin{equation} \label{growth}
%\begin{split}
%|\tilde A(x,\xi)| &\leq c_1|\xi|^{q-1} + k(x),\\
%\tilde A(x,\xi)\cdot\xi &\geq c_2|\xi|^{q} -m(x),
%\end{split}
%\end{equation}
%where $c_1$, $c_2>0$,
% $k\in L^{q'}(\Omega)$ and $m\in L^1(\Omega)$ are nonnegative
%functions.
\end{enumerate}
Moreover, let
$U$ be a dense set in $\Rd$ and
%$\tilde A(t,x,\xi)$ be measurable w.r.t. ${\mLl(Q)} \otimes {\mB(\Rd})$ $\sigma$--algebra, and
$(\bB,\tilde A(t,x,\bB)) \in \mA(t,x)$ for a.a. $(t,x)\in Q$ and for all $\bB\in U$. Let also
  $(\nabla u, \bS)\in \Rd \times \Rd$. Then the following conditions are equivalent:
\begin{equation}\begin{split}
\textrm{(i)} \quad &({\bS} - \tilde A(t,x,\bB))\cdot( {\nabla u} - \bB) \geq 0 \quad \textrm{ for all } \quad (\bB,\tilde A(t,x,\bB))\in\mA(t,x)\,,\\
\textrm{(ii)} \quad & (\nabla u, \bS) \in \mA(t,x).
\end{split} \label{maxmon}
\end{equation}
\end{lemma}
For the proof see~\cite{BuGwMaSw2012}.

\bigskip
\noindent
{\bf Acknowledgements }\\
The author was supported by the grant IdP2011/000661.

%\subsection{Convergence tools}
%
%In this subsection, we formulate simple criteria that allows us to verify that $\bS$ and $\nabla u$ fulfill the nonlinear constitutive equation \eqref{I2}, which is equivalent to $(\nabla u,\bS)\in\mA$. The first criterion concerns in particular the case when $\psi$ and $\psi^*$ fulfill {\bf (H1)} and {\bf (H2)}. The second criterion concerns the general (non-reflexive) case. We state the first lemma since it is sufficient for many  applications and its proof is simpler than the proof for the  general case.
%
%\begin{lemma}\label{Minty}
%Let $\mathcal{A}$ be maximal monotone $\psi$-graph
%satisfying {\bf (A1)}--{\bf (A4)}. Assume that there are
%sequences $\{\bS^n\}_{n=1}^\infty$ and
%$\{\nabla u^n\}_{n=1}^{\infty}$ defined on a measurable set $G\subset \R^k$ such
%the following holds:
%\begin{align}
%(\nabla u^n,\bS^n)& \in \mathcal{A} &&\textrm{ a.e. in }G,\label{Ass1}\\
%\nabla u^n &\weak \nabla u &&\textrm{ weakly}\  \textrm{in } L^\psi(G),\label{1.16}\\
%\bS^n &\weak \bS &&\textrm{ weakly}\  \textrm{in }  L^{\psi^*}(G),\label{1.17}\\
%\limsup_{n\to \infty} \int_{G}\bS^n \cdot \nabla u^n \,d z &\le \int_{G} \bS \cdot \nabla u \, d z.\label{Ass2}\\
%%\end{align}
%\intertext{Then}
%(\nabla u,\bS)&\in \mathcal{A} &&\textrm{ a.e. in }G. \label{Minty2eq}
%\end{align}
%\end{lemma}
%\begin{proof}
%We first observe that \eqref{Ass1}--\eqref{Ass2} imply that
%\begin{equation}
%\limsup\limits_{n\to\infty}
%\int_{G}(\bS^n-\tilde A(\nabla u))\cdot (\nabla u^n-\nabla u)\,d z\le 0.\label{Minty3}
%\end{equation}
%Since the graph is monotone, \eqref{Minty3} is equivalent to
%\begin{equation}
%\limsup\limits_{n\to\infty}
%\int_{G}\left|(\bS^n-\tilde A(\nabla u))\cdot (\nabla u^n-\nabla u)\right|\,d z= 0.
%\end{equation}
%Therefore, $(\bS^n-\tilde A(\nabla u))\cdot (\nabla u^n-\nabla u)$ converges strongly in $L^1(G)$ and consequently weakly in $L^1(G)$, namely, we have  for all $\varphi\in L^\infty(G)$
%\begin{equation}\label{weakly}
%\lim\limits_{n\to\infty}
%\int_{G}(\bS^n-\tilde A(\nabla u))\cdot (\nabla u^n
%-\nabla u)\varphi\,d z=0.
%\end{equation}
%From \eqref{weakly}, \eqref{1.16} and \eqref{1.17} it follows that
%\begin{equation}\label{limit_with_phi}
%\lim\limits_{n\to\infty}\int_{G} \bS^n\cdot \nabla u^n\varphi\,d z=
%\lim\limits_{n\to\infty}\int_{G} \bS^n\cdot \nabla u\varphi \,d z=
%\int_{G}\bS\cdot \nabla u \varphi \,d z.
%\end{equation}
%Consequently, since the graph is monotone, we observe that for an arbitrary fixed matrix $\bB\in \mathbb{R}^{d\times d}_{sym}$ and all nonnegative $\varphi\in L^\infty(G)$
%\begin{equation*}
%0\le \lim\limits_{n\to\infty}\int_{G}(\bS^n-\tilde A(\bB))\cdot(\nabla u^n-\bB)\varphi \,d z =
%\int_{G}(\bS-\tilde A(\bB))\cdot(\nabla u-\bB)\varphi \,d z.
%\end{equation*}
%Since $\varphi\ge0$ is arbitrary, this implies that  %and a.a. $(t,x)\in Q'$
%\begin{equation}\label{mon-mon}
%(\bS-\tilde A(\bB))\cdot(\nabla u-\bB)\ge0\  \textrm{ for all } \bB\in\R^{d\times d}_{\rm sym} \textrm{ a.e. in } G,
%\end{equation}
%and the assertion follows  from {\bf (A3$^*$)}.
%\end{proof}
%
%
%{\bf Proof of Lemma \ref{Minty2}}
%\begin{proof}
%Since  the graph is monotone, we have
%\begin{equation}\label{43}
%\int_{G}(\tilde A(\bB)-\bS^n))\cdot(\bB-\nabla u^n)\,d z \geq 0 \quad \textrm{for all }\, \bB
%\in L^{\infty}(G).
%\end{equation}
%Note that if $\bB\in L^\infty(G)$ then  $\tilde A(\bB) \in L^{\infty}(G)$. Indeed, assume the opposite, i.e.,
% $\tilde A(\bB)$ is unbounded. Then, it follows from {\bf (A4)} and the fact that $\psi$ is nonnegative that
%$$
%|\bB|\ge c_*\frac{\psi^\ast(\tilde A(\bB))-d_*}{|\tilde A(\bB)|}.
%$$
%Since $\psi^*$ is an $N-$function, then the right-hand side tends to infinity, which contradicts the fact that $\bB \in L^{\infty}(G)$.
%
%Next, rearranging  \eqref{43} we have
%\begin{equation}\label{44}
%\int_{G} \bS^n\cdot \nabla u^n \,d z  \geq \int_{G}
%\bS^n\cdot\bB \,d z + \int_{G} \tilde A(\bB)\cdot(\nabla u^n-\bB)\,d z.
%\end{equation}
%Letting  $n\to\infty$ in \eqref{44}, we conclude from \eqref{1.26}--\eqref{Ass2-2} that
%\begin{equation}\label{45}
%\int_{G}\bS\cdot \nabla u \,d z \geq \int_{G} \bS\cdot\bB \,d z +
%\int_{G} \tilde A(\bB)\cdot(\nabla u-\bB)\,d z
%\end{equation}
%and consequently
%\begin{equation}\label{46}
%\int_{G} (\tilde A(\bB)-\bS)\cdot(\bB-\nabla u)\,d z \geq 0 \quad \textrm{for all }\, \bB
%\in L^{\infty}(G).
%\end{equation}
%Next, we use a generalization of  the Minty method to prove \eqref{Minty2-2}. For  any $j>0$ we define the set
%$G_{j} := \{z\in G;  |\nabla u(z)|\leq j\}
%$
%and use $ \bbbone_{G_{j}}$ to denote the characteristic function of $G_j$. Since $|\nabla u| \in L^1(G)$ due to \eqref{1.26} we observe that
%\begin{equation}
%|G\setminus G_j|\le \frac{C}{j}. \label{smallset}
%\end{equation}
%Next, let $i,j\in\N$
%be arbitrary such that $0<j<i$. Then we take $B$ in \eqref{46} of the form
%$$
%\bB:=\nabla u \bbbone_{G_{i}} + h\,W\bbbone_{G_{j}},
%$$
%where  $h>0$  and $W\in L^{\infty}(G)$ is arbitrary (note that $\bB$ is bounded). Doing so, we obtain
%\begin{equation}\label{13}
%\int_{G_{j}}(\tilde A(\nabla u+h W)-\bS)\cdot W \,d z \geq \frac{1}{h}\int_{G \setminus G_{i}} \left(\tilde A(\b0)\cdot \nabla u - \bS\cdot \nabla u\right) \,d z.
%\end{equation}
%Next, using \eqref{1.26}--\eqref{1.27}, and the H\"{o}lder inequality we obtain
%$$
%\int_{G} \left|\tilde A(\b0)\cdot \nabla u - \bS\cdot \nabla u\right|\,d z<\infty.
%$$
%Consequently, we appeal to the Lebesgue dominated convergence theorem and \eqref{smallset} to let $i\to \infty$ on the right-hand side of \eqref{13} and conclude that
%$$
%\lim\limits_{i\to\infty}\int_{G \setminus G_{i}}
%\left(\tilde A(\b0)\cdot \nabla u+\bS\cdot\nabla u\right) \,d z=0,
%$$
%which directly implies that
%\begin{equation}\label{MB1-1}
%\int_{G_{j}}(\tilde A(\nabla u+h W)-\bS)\cdot W \,d z\geq 0\quad \textrm{for all }j\in\N\,.
%\end{equation}
%Our goal is to let $h\to 0_+$ and to show \eqref{Minty2-2}. Using the definition of $G_j$ it is easy to see that for a subsequence
%\begin{align*}
%\tilde A(\nabla u+h W)&\rightharpoonup \bar{\bS} &&\textrm{weakly in }L^2(G_j),\\
%\nabla u+h W &\to \nabla u &&\textrm{strongly  in }L^2(G_j),\\
%(\nabla u+h W,\tilde A(\nabla u+h W))&\in {\mathcal A} &&\textrm{a.e. in }G_j.
%\end{align*}
%Therefore, we can use Lemma~\ref{Minty} (note that \eqref{Ass2} is direct consequence of above convergence results) with $L^\psi=L^2$ and we observe that
%\begin{equation}
%(\nabla u,\bar\bS)\in\mA \qquad \textrm{ a.e. in }G_j\,. \label{begr}
%\end{equation}
%Moreover, since $\nabla u$ is bounded in $G_j$ then also $\bar\bS$ is bounded in $G_j$ due to \eqref{begr}. Finally, taking the limit $h\to 0_+$ in \eqref{MB1-1}, we have
%$$
%\int_{G_j}(\bar \bS-\bS)\cdot W \,d z \ge 0 \qquad \textrm{for all } W\in L^{\infty}(G_j).
%$$
%Setting,  $ W:=-\frac{(\bar\bS-\bS)}{|\bar\bS-\bS|}\bbbone _{\{\bar\bS\neq\bS\}}$   yields
%$$
%\int_{G_j}|\bar\bS-\bS|\,d z\leq0
%$$
%and therefore \eqref{begr} implies that  $(\nabla u,\bS)\in\mA$ a.e. in $G_j$. But since $j$ was arbitrary, we use \eqref{smallset} and conclude that
% $(\nabla u,\bS)\in\mA$ a.e. in $G$.
%\end{proof}
%
%WPISANE TW O WLOZENIACH DO WYKORZYSTANIA W DRUGIEJ PRACY


%\section{Embeddings}
%\begin{theorem}[\cite{Cianchi}]
%Let $M:\R^d\to\R_+$ be an $N-$function such that for every $\lambda>0, \{\xi\in\R^d: M(\xi)\le\lambda\}$
%is a compact set whose interior contains zero.  Denote by $M_*:\R_+\to\R_+$ the left-continuous increasing function satisfying
%\begin{equation}
%|\{\xi\in\R^d:M(\xi)\le 1\}|=|\{\xi\in\R^d:M_*(|\xi|)\le t\}|
%\end{equation}
%for every $t\ge0$. Assume that 
%\begin{equation}
%\int_{\R_+}\left(\frac{t}{M_*(t)}\right)^{d'-1}dt<\infty
%\end{equation}
%where $d'=d/(d-1)$ and define $H:\R_+\to\R_+$ as follows
%\begin{equation}
%H(r)= \left(\int_0^r \left(\frac{t}{M_*(t)}\right)^{d'-1}dt \right)^{1/d'}
%\end{equation}
%for $r\ge0$. Moreover,  $B$ is given by 
%\begin{equation}\label{B}
%B=M_*\circ H^{-1}
%\end{equation}
%where $H^{-1}$ is the left-continuous inverse of $H$. Then there exists a constant $K$, depending only on $d$ such  that
%\begin{equation}
%\|u\|_{B, \R^d}\le K\|\nabla u\|_{M,\R^d}.
%\end{equation}
%\end{theorem}

\noindent
%{\bf Acknowledgments.}


\bibliographystyle{abbrv}
\bibliography{genorlicz}


\end{document}












swierczewska_revised.bbl
\def\ocirc#1{\ifmmode\setbox0=\hbox{$#1$}\dimen0=\ht0 \advance\dimen0
  by1pt\rlap{\hbox to\wd0{\hss\raise\dimen0
  \hbox{\hskip.2em$\scriptscriptstyle\circ$}\hss}}#1\else {\accent"17 #1}\fi}
\begin{thebibliography}{10}

\bibitem{Be2011}
A.~Benkirane, J.~Douieb, and M.~Ould Mohamedhen~Val.
\newblock An approximation theorem in {M}usielak-{O}rlicz-{S}obolev spaces.
\newblock {\em Comment. Math.}, 51(1):109--120, 2011.

\bibitem{BeEl1999}
A.~Benkirane and A.~Elmahi.
\newblock An existence theorem for a strongly nonlinear elliptic problem in
  {O}rlicz spaces.
\newblock {\em Nonlinear Anal.}, 36(1, Ser. A: Theory Methods):11--24, 1999.

\bibitem{BuGwMaRaSw2012}
M.~Bul{\'\i}\v{c}ek, P.~Gwiazda, J.~M{\'a}lek, K.~R. Rajagopal, and
  A.~{\'S}wierczewska-Gwiazda.
\newblock On flows of fluids described by an implicit constitutive equation
  characterized by a maximal monotone graph.
\newblock In {\em Mathematical Aspects of Fluid Mechanics}, volume 402 of {\em
  London Math. Soc. Lecture Note Ser.}, pages 23--51. Cambridge Univ. Press,
  2012.

\bibitem{BuGwMaSw2009}
M.~Bul{\'\i}\v{c}ek, P.~Gwiazda, J.~M{\'a}lek, and A.~{\'S}wierczewska-Gwiazda.
\newblock {On steady flows of incompressible fluids with implicit
  power-law-like rheology.}
\newblock {\em Adv. Calc. Var.}, 2(2):109--136, 2009.

\bibitem{BuGwMaSw2012}
M.~Bul{\'\i}\v{c}ek, P.~Gwiazda, J.~M{\'a}lek, and A.~{\'S}wierczewska-Gwiazda.
\newblock On unsteady flows of implicitly constituted incompressible fluids.
\newblock {\em SIAM J. Math. Anal.}, 44(4):2756--2801, 2012.

\bibitem{Donaldson}
T.~Donaldson.
\newblock Inhomogeneous {O}rlicz-{S}obolev spaces and nonlinear parabolic
  initial value problems.
\newblock {\em J. Differential Equations}, 16:201--256, 1974.

\bibitem{ElmahiMeskine}
A.~Elmahi and D.~Meskine.
\newblock Parabolic equations in {O}rlicz spaces.
\newblock {\em J. London Math. Soc. (2)}, 72(2):410--428, 2005.

\bibitem{ElMe2005}
A.~Elmahi and D.~Meskine.
\newblock {Strongly nonlinear parabolic equations with natural growth terms in
  Orlicz spaces.}
\newblock {\em Nonlinear Anal., Theory Methods Appl.}, 60(1):1--35, 2005.

\bibitem{Gossez1}
J.-P. Gossez.
\newblock Nonlinear elliptic boundary value problems for equations with rapidly
  (or slowly) increasing coefficients.
\newblock {\em Trans. Amer. Math. Soc.}, 190:163--205, 1974.

\bibitem{Gossez2}
J.-P. Gossez.
\newblock Some approximation properties in {O}rlicz-{S}obolev spaces.
\newblock {\em Studia Math.}, 74(1):17--24, 1982.

\bibitem{GwMaSw2007}
P.~Gwiazda, J.~M{\'a}lek, and A.~{\'S}wierczewska.
\newblock {On flows of an incompressible fluid with a discontinuous
  power-law-like rheology.}
\newblock {\em Comput. Math. Appl.}, 53(3-4):531--546, 2007.

\bibitem{GwMiSw2013}
P.~Gwiazda, P.~Minakowski, and A.~\'Swierczewska-Gwiazda.
\newblock On the anisotropic {O}rlicz spaces applied in the problems of
  continuum mechanics.
\newblock {\em Discrete Contin. Dyn. Syst. Ser. S}, 6(5):1291--1306, 2013.

\bibitem{GwMiWr2012}
P.~Gwiazda, P.~Minakowski, and A.~Wr{{\'o}}blewska-Kami{{\'n}}ska.
\newblock Elliptic problems in generalized {O}rlicz-{M}usielak spaces.
\newblock {\em Cent. Eur. J. Math.}, 10(6):2019--2032, 2012.

\bibitem{GwSw2008}
P.~Gwiazda and A.~{\'S}wierczewska-Gwiazda.
\newblock On non-{N}ewtonian fluids with a property of rapid thickening under
  different stimulus.
\newblock {\em Math. Models Methods Appl. Sci.}, 18(7):1073--1092, 2008.

\bibitem{GwSw2008a}
P.~Gwiazda and A.~{\'S}wierczewska-Gwiazda.
\newblock {On steady non-Newtonian fluids with growth conditions in generalized
  Orlicz spaces.}
\newblock {\em Topol. Methods Nonlinear Anal.}, 32(1):103--113, 2008.

\bibitem{GwSw2010}
P.~Gwiazda and A.~\'Swierczewska-Gwiazda.
\newblock Parabolic equations in anisotropic {O}rlicz spaces with general
  {N}-functions.
\newblock {\em special Birkhaeuser volume ''Parabolic Problems. The Herbert
  Amann Festschrift'', Progress in Nonlinear Differential Equations and Their
  Applications}, 60:301--311, 2010.

\bibitem{GwSwWr2010}
P.~Gwiazda, A.~{\'S}wierczewska-Gwiazda, and A.~Wr{{\'o}}blewska.
\newblock Monotonicity methods in generalized {O}rlicz spaces for a class of
  non-{N}ewtonian fluids.
\newblock {\em Math. Methods Appl. Sci.}, 33(2):125--137, 2010.

\bibitem{GwSwWr2012}
P.~Gwiazda, A.~{\'S}wierczewska-Gwiazda, and A.~Wr{{\'o}}blewska.
\newblock Generalized {S}tokes system in {O}rlicz spaces.
\newblock {\em Discrete Contin. Dyn. Syst.}, 32(6):2125--2146, 2012.

\bibitem{GwZa2007}
P.~Gwiazda and A.~Zatorska-Goldstein.
\newblock {On elliptic and parabolic systems with $x$-dependent multivalued
  graphs.}
\newblock {\em Math. Methods Appl. Sci.}, 30(2):213--236, 2007.

\bibitem{Musielak}
J.~Musielak.
\newblock {\em Orlicz spaces and modular spaces}, volume 1034 of {\em Lecture
  Notes in Mathematics}.
\newblock Springer-Verlag, Berlin, 1983.

\bibitem{Novotny}
A.~Novotn\'y and I.~Stra\v{s}kraba.
\newblock {\em Introduction to the mathematical theory of compressible flow.}
\newblock Oxford Lecture Series in Mathematics and its Applications 27. Oxford:
  Oxford University Press., 2004.

\bibitem{Sw2013}
A.~{\'S}wierczewska-Gwiazda.
\newblock Nonlinear parabolic problems in {M}usielak--{O}rlicz spaces.
\newblock {\em http://arxiv.org/abs/1306.2186}, 06 2013.

\bibitem{Wr2010}
A.~Wr{{\'o}}blewska.
\newblock Steady flow of non-{N}ewtonian fluids---monotonicity methods in
  generalized {O}rlicz spaces.
\newblock {\em Nonlinear Anal.}, 72(11):4136--4147, 2010.

\end{thebibliography}
