2 #ifndef DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONCCFV_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_CONVECTIONDIFFUSIONCCFV_HH
5 #include<dune/common/exceptions.hh>
6 #include<dune/common/fvector.hh>
7 #include<dune/common/typetraits.hh>
8 #include<dune/geometry/referenceelements.hh>
66 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
67 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
70 const auto& cell = eg.entity();
73 auto geo = eg.geometry();
74 auto ref_el = referenceElement(geo);
75 auto local_inside = ref_el.position(0,0);
78 auto c = param.c(cell,local_inside);
81 r.accumulate(lfsu,0,(c*x(lfsu,0))*geo.volume());
85 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
92 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
93 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
97 const auto& cell = eg.entity();
100 auto geo = eg.geometry();
101 auto ref_el = referenceElement(geo);
102 auto local_inside = ref_el.position(0,0);
105 auto c = param.c(cell,local_inside);
108 mat.accumulate(lfsu,0,lfsu,0,c*geo.volume());
113 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
115 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
116 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
117 R& r_s, R& r_n)
const
120 using RF =
typename LFSU::Traits::FiniteElementType::
121 Traits::LocalBasisType::Traits::RangeFieldType;
124 const auto dim = IG::Entity::dimension;
127 auto cell_inside =
ig.inside();
128 auto cell_outside =
ig.outside();
131 auto geo =
ig.geometry();
132 auto geo_inside = cell_inside.geometry();
133 auto geo_outside = cell_outside.geometry();
136 auto geo_in_inside =
ig.geometryInInside();
139 auto ref_el = referenceElement(geo);
140 auto face_local = ref_el.position(0,0);
143 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
146 auto ref_el_inside = referenceElement(geo_inside);
147 auto ref_el_outside = referenceElement(geo_outside);
148 auto local_inside = ref_el_inside.position(0,0);
149 auto local_outside = ref_el_outside.position(0,0);
152 auto tensor_inside = param.A(cell_inside,local_inside);
153 auto tensor_outside = param.A(cell_outside,local_outside);
154 auto n_F =
ig.centerUnitOuterNormal();
155 Dune::FieldVector<RF,dim> An_F;
156 tensor_inside.mv(n_F,An_F);
157 auto k_inside = n_F*An_F;
158 tensor_outside.mv(n_F,An_F);
159 auto k_outside = n_F*An_F;
160 auto k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
163 auto iplocal_s = geo_in_inside.global(face_local);
164 auto b = param.b(cell_inside,iplocal_s);
167 if (vn>=0) u_upwind = x_s(lfsu_s,0);
else u_upwind = x_n(lfsu_n,0);
170 auto global_inside = geo_inside.global(local_inside);
171 auto global_outside = geo_outside.global(local_outside);
174 global_inside -= global_outside;
175 auto distance = global_inside.two_norm();
178 r_s.accumulate(lfsu_s,0,(u_upwind*vn)*face_volume+k_avg*(x_s(lfsu_s,0)-x_n(lfsu_n,0))*face_volume/distance);
179 r_n.accumulate(lfsu_n,0,-(u_upwind*vn)*face_volume-k_avg*(x_s(lfsu_s,0)-x_n(lfsu_n,0))*face_volume/distance);
183 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
185 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
186 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
187 Y& y_s, Y& y_n)
const
192 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
194 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
195 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
196 M& mat_ss, M& mat_sn,
197 M& mat_ns, M& mat_nn)
const
200 using RF =
typename LFSU::Traits::FiniteElementType::
201 Traits::LocalBasisType::Traits::RangeFieldType;
204 const auto dim = IG::Entity::dimension;
207 auto cell_inside =
ig.inside();
208 auto cell_outside =
ig.outside();
211 auto geo =
ig.geometry();
212 auto geo_inside = cell_inside.geometry();
213 auto geo_outside = cell_outside.geometry();
216 auto geo_in_inside =
ig.geometryInInside();
219 auto ref_el = referenceElement(geo);
220 auto face_local = ref_el.position(0,0);
223 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
226 auto ref_el_inside = referenceElement(geo_inside);
227 auto ref_el_outside = referenceElement(geo_outside);
228 auto local_inside = ref_el_inside.position(0,0);
229 auto local_outside = ref_el_outside.position(0,0);
232 auto tensor_inside = param.A(cell_inside,local_inside);
233 auto tensor_outside = param.A(cell_outside,local_outside);
234 auto n_F =
ig.centerUnitOuterNormal();
235 Dune::FieldVector<RF,dim> An_F;
236 tensor_inside.mv(n_F,An_F);
237 auto k_inside = n_F*An_F;
238 tensor_outside.mv(n_F,An_F);
239 auto k_outside = n_F*An_F;
240 auto k_avg = 2.0/(1.0/(k_inside+1E-30) + 1.0/(k_outside+1E-30));
243 auto iplocal_s = geo_in_inside.global(face_local);
244 auto b = param.b(cell_inside,iplocal_s);
248 auto global_inside = geo_inside.global(local_inside);
249 auto global_outside = geo_outside.global(local_outside);
252 global_inside -= global_outside;
253 auto distance = global_inside.two_norm();
256 mat_ss.accumulate(lfsu_s,0,lfsu_s,0, k_avg*face_volume/distance );
257 mat_ns.accumulate(lfsu_n,0,lfsu_s,0, -k_avg*face_volume/distance );
258 mat_sn.accumulate(lfsu_s,0,lfsu_n,0, -k_avg*face_volume/distance );
259 mat_nn.accumulate(lfsu_n,0,lfsu_n,0, k_avg*face_volume/distance );
262 mat_ss.accumulate(lfsu_s,0,lfsu_s,0, vn*face_volume );
263 mat_ns.accumulate(lfsu_n,0,lfsu_s,0, -vn*face_volume );
267 mat_sn.accumulate(lfsu_s,0,lfsu_n,0, vn*face_volume );
268 mat_nn.accumulate(lfsu_n,0,lfsu_n,0, -vn*face_volume );
274 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
276 const LFSV& lfsv, R& r)
const
278 if (not first_stage)
return;
281 const auto& cell = eg.entity();
284 auto geo = eg.geometry();
285 auto ref_el = referenceElement(geo);
286 auto local_inside = ref_el.position(0,0);
289 auto cellcapacity = param.d(cell,local_inside)*geo.volume();
290 auto celldt = cellcapacity/(cellinflux+1E-30);
291 dtmin = std::min(dtmin,celldt);
297 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
299 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
303 using RF =
typename LFSU::Traits::FiniteElementType::
304 Traits::LocalBasisType::Traits::RangeFieldType;
307 const auto dim = IG::Entity::dimension;
310 auto cell_inside =
ig.inside();
313 auto geo =
ig.geometry();
314 auto geo_inside = cell_inside.geometry();
317 auto geo_in_inside =
ig.geometryInInside();
320 auto ref_el = referenceElement(geo);
321 auto face_local = ref_el.position(0,0);
324 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
327 auto ref_el_inside = referenceElement(geo_inside);
328 auto local_inside = ref_el_inside.position(0,0);
331 auto bctype = param.bctype(
ig.intersection(),face_local);
337 auto global_inside = geo_inside.global(local_inside);
338 auto global_outside = geo.global(face_local);
339 global_inside -= global_outside;
340 auto distance = global_inside.two_norm();
343 auto tensor_inside = param.A(cell_inside,local_inside);
344 auto n_F =
ig.centerUnitOuterNormal();
345 Dune::FieldVector<RF,dim> An_F;
346 tensor_inside.mv(n_F,An_F);
347 auto k_inside = n_F*An_F;
350 auto iplocal_s = geo_in_inside.global(face_local);
351 auto g = param.g(cell_inside,iplocal_s);
354 auto b = param.b(cell_inside,iplocal_s);
355 auto n =
ig.centerUnitOuterNormal();
358 r_s.accumulate(lfsu_s,0,(b*n)*g*face_volume + k_inside*(x_s(lfsu_s,0)-g)*face_volume/distance);
369 auto j = param.j(
ig.intersection(),face_local);
372 r_s.accumulate(lfsu_s,0,j*face_volume);
380 auto iplocal_s = geo_in_inside.global(face_local);
381 auto b = param.b(cell_inside,iplocal_s);
382 auto n =
ig.centerUnitOuterNormal();
385 auto o = param.o(
ig.intersection(),face_local);
388 r_s.accumulate(lfsu_s,0,((b*n)*x_s(lfsu_s,0) + o)*face_volume);
394 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
396 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
400 using RF =
typename LFSU::Traits::FiniteElementType::
401 Traits::LocalBasisType::Traits::RangeFieldType;
404 const auto dim = IG::Entity::dimension;
407 auto cell_inside =
ig.inside();
410 auto geo =
ig.geometry();
411 auto geo_inside = cell_inside.geometry();
414 auto geo_in_inside =
ig.geometryInInside();
417 auto ref_el = referenceElement(geo);
418 auto face_local = ref_el.position(0,0);
421 auto face_volume = geo.integrationElement(face_local) * ref_el.volume();
424 auto ref_el_inside = referenceElement(geo_inside);
425 auto local_inside = ref_el_inside.position(0,0);
428 auto bctype = param.bctype(
ig.intersection(),face_local);
434 auto global_inside = geo_inside.global(local_inside);
435 auto global_outside = geo.global(face_local);
436 global_inside -= global_outside;
437 auto distance = global_inside.two_norm();
440 auto tensor_inside = param.A(cell_inside,local_inside);
441 auto n_F =
ig.centerUnitOuterNormal();
442 Dune::FieldVector<RF,dim> An_F;
443 tensor_inside.mv(n_F,An_F);
444 auto k_inside = n_F*An_F;
447 mat_ss.accumulate(lfsu_s,0,lfsv_s,0, k_inside*face_volume/distance );
455 auto iplocal_s = geo_in_inside.global(face_local);
456 auto b = param.b(cell_inside,iplocal_s);
457 auto n =
ig.centerUnitOuterNormal();
460 mat_ss.accumulate(lfsu_s,0,lfsv_s,0, (b*n)*face_volume );
467 template<
typename EG,
typename LFSV,
typename R>
471 const auto& cell = eg.entity();
474 auto geo = eg.geometry();
475 auto ref_el = referenceElement(geo);
476 auto local_inside = ref_el.position(0,0);
479 auto f = param.f(cell,local_inside);
481 r.accumulate(lfsv,0,-f*geo.volume());
485 void setTime (
typename TP::Traits::RangeFieldType t)
491 void preStep (
typename TP::Traits::RangeFieldType time,
typename TP::Traits::RangeFieldType dt,
497 void preStage (
typename TP::Traits::RangeFieldType time,
int r)
504 else first_stage =
false;
513 typename TP::Traits::RangeFieldType
suggestTimestep (
typename TP::Traits::RangeFieldType dt)
const
521 mutable typename TP::Traits::RangeFieldType dtmin;
522 mutable typename TP::Traits::RangeFieldType cellinflux;
553 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
554 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
557 const auto& cell = eg.entity();
560 auto geo = eg.geometry();
561 auto ref_el = referenceElement(geo);
562 auto local_inside = ref_el.position(0,0);
565 auto capacity = param.d(cell,local_inside);
568 r.accumulate(lfsu,0,capacity*x(lfsu,0)*geo.volume());
572 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename Y>
579 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
584 const auto& cell = eg.entity();
587 auto geo = eg.geometry();
588 auto ref_el = referenceElement(geo);
589 auto local_inside = ref_el.position(0,0);
592 auto capacity = param.d(cell,local_inside);
595 mat.accumulate(lfsu,0,lfsu,0,capacity*geo.volume());
static const int dim
Definition: adaptivity.hh:84
const IG & ig
Definition: constraints.hh:149
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Definition: convectiondiffusionccfv.hh:44
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:67
void postStage()
to be called once at the end of each stage
Definition: convectiondiffusionccfv.hh:508
void jacobian_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, M &mat_ss) const
Definition: convectiondiffusionccfv.hh:395
@ doAlphaVolume
Definition: convectiondiffusionccfv.hh:53
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusionccfv.hh:93
void setTime(typename TP::Traits::RangeFieldType t)
set time in parameter class
Definition: convectiondiffusionccfv.hh:485
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:468
TP::Traits::RangeFieldType suggestTimestep(typename TP::Traits::RangeFieldType dt) const
to be called once before each stage
Definition: convectiondiffusionccfv.hh:513
@ doAlphaSkeleton
Definition: convectiondiffusionccfv.hh:54
void preStep(typename TP::Traits::RangeFieldType time, typename TP::Traits::RangeFieldType dt, int stages)
to be called once before each time step
Definition: convectiondiffusionccfv.hh:491
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Definition: convectiondiffusionccfv.hh:86
@ doAlphaBoundary
Definition: convectiondiffusionccfv.hh:55
@ doPatternVolume
Definition: convectiondiffusionccfv.hh:49
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusionccfv.hh:298
@ doPatternSkeleton
Definition: convectiondiffusionccfv.hh:50
@ doLambdaVolume
Definition: convectiondiffusionccfv.hh:56
void preStage(typename TP::Traits::RangeFieldType time, int r)
to be called once before each stage
Definition: convectiondiffusionccfv.hh:497
void alpha_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, R &r_s, R &r_n) const
Definition: convectiondiffusionccfv.hh:114
void alpha_volume_post_skeleton(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:275
ConvectionDiffusionCCFV(TP ¶m_)
Definition: convectiondiffusionccfv.hh:60
@ doLambdaSkeleton
Definition: convectiondiffusionccfv.hh:57
@ doLambdaBoundary
Definition: convectiondiffusionccfv.hh:58
void jacobian_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, M &mat_ss, M &mat_sn, M &mat_ns, M &mat_nn) const
Definition: convectiondiffusionccfv.hh:193
void jacobian_apply_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, Y &y_s, Y &y_n) const
Definition: convectiondiffusionccfv.hh:184
Definition: convectiondiffusionccfv.hh:539
void jacobian_apply_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, Y &y) const
Definition: convectiondiffusionccfv.hh:573
@ doAlphaVolume
Definition: convectiondiffusionccfv.hh:545
ConvectionDiffusionCCFVTemporalOperator(TP ¶m_)
Definition: convectiondiffusionccfv.hh:547
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionccfv.hh:554
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusionccfv.hh:580
@ doPatternVolume
Definition: convectiondiffusionccfv.hh:542
Type
Definition: convectiondiffusionparameter.hh:113
@ Neumann
Definition: convectiondiffusionparameter.hh:113
@ Outflow
Definition: convectiondiffusionparameter.hh:113
@ Dirichlet
Definition: convectiondiffusionparameter.hh:113
Default flags for all local operators.
Definition: flags.hh:19
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericaljacobianapply.hh:287
sparsity pattern generator
Definition: pattern.hh:14
sparsity pattern generator
Definition: pattern.hh:30