Commit 518c0906 authored by Lukas Riedel's avatar Lukas Riedel 📝

implemented jacobian_boundary. still not working

parent e8a76d10
......@@ -57,7 +57,7 @@ template<typename Traits, typename Parameter, typename Boundary, typename Source
:
// public Dune::PDELab::NumericalJacobianVolume<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint> >,
// public Dune::PDELab::NumericalJacobianSkeleton<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint> >,
public Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint> >,
// public Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint> >,
public Dune::PDELab::FullSkeletonPattern,
public Dune::PDELab::FullVolumePattern,
public Dune::PDELab::LocalOperatorDefaultFlags,
......@@ -130,7 +130,7 @@ public:
:
// Dune::PDELab::NumericalJacobianVolume<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
// Dune::PDELab::NumericalJacobianSkeleton<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
// Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
param(param_), boundary(boundary_), sourceTerm(sourceTerm_), method(method_), weights(weights_),
penalty_factor(config.get<RF>("dg.penaltyFactor")), mapper(view_),
intorderadd(intorderadd_), quadrature_factor(quadrature_factor_),
......@@ -576,14 +576,14 @@ public:
mat_ss.accumulate(lfsv_s,i,lfsu_s,j,
relCond_upwind
* ( omega_s * satCond_s
* ( - ( tgradphi_s[j] * normal ) + penalty * phi_s[j] )
* ( - ( tgradphi_s[j] * normal ) + penalty * phi_s[j] )
+ omega_n * satCond_n * penalty * phi_s[j] )
* phi_s[i] * factor
);
// differentiate symmetry term
mat_ss.accumulate(lfsv_s,i,lfsu_s,j,
theta * omega_s * (tgradphi_s[i] * normal) * satCond_s
* phi_s[j] * ( rel_cond_upwind_deriv_s + 1 )
theta * omega_s * satCond_s * (tgradphi_s[i] * normal)
* phi_s[j] * ( rel_cond_upwind_deriv_s + 1.0 )
* factor
);
}
......@@ -604,7 +604,7 @@ public:
// differentiate symmetry term
mat_nn.accumulate(lfsv_n,i,lfsu_n,j,
theta * omega_n * (tgradphi_n[i] * normal) * satCond_n
* phi_n[j] * ( rel_cond_upwind_deriv_n - 1 )
* phi_n[j] * ( rel_cond_upwind_deriv_n - 1.0 )
* factor
);
}
......@@ -625,7 +625,7 @@ public:
// differentiate symmetry term
mat_sn.accumulate(lfsv_s,i,lfsu_n,j,
theta * omega_s * (tgradphi_s[i] * normal) * satCond_s
* phi_n[j] * ( rel_cond_upwind_deriv_n - 1 )
* phi_n[j] * ( rel_cond_upwind_deriv_n - 1.0 )
* factor
);
}
......@@ -646,7 +646,7 @@ public:
// differentiate symmetry term
mat_ns.accumulate(lfsv_n,i,lfsu_s,j,
theta * omega_n * (tgradphi_n[i] * normal) * satCond_n
* phi_s[j] * ( rel_cond_upwind_deriv_s + 1 )
* phi_s[j] * ( rel_cond_upwind_deriv_s + 1.0 )
* factor
);
}
......@@ -899,6 +899,223 @@ public:
}
}
template<typename IG, typename LFSU, typename X, typename LFSV, typename LocalMatrix>
void jacobian_boundary (const IG& ig,
const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
LocalMatrix& mat_ss) const
{
// get polynomial degree
const int order_s = lfsu_s.finiteElement().localBasis().order();
const int degree = order_s;
const int intorder = intorderadd + quadrature_factor * degree;
// geometric factor of penalty
const RF h_F = ig.inside().geometry().volume()/ig.geometry().volume(); // Houston!
// get element geometry
auto gtface = ig.geometry();
// loop over quadrature points and integrate normal flux
std::size_t i = 0;
for (const auto& it : quadratureRule(gtface,intorder))
{
// calculate unique id
const auto it_id = create_unique_id(ig.intersection(),i);
i++;
// position of quadrature point in local coordinates of elements
const auto p_local_s = ig.geometryInInside().global(it.position());
// evaluate basis functions
const auto& phi_s = cache[order_s].evaluateFunction(p_local_s,lfsu_s.finiteElement().localBasis());
// integration factor
const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
// evaluate u
RF u_s = 0.;
for (unsigned int i = 0; i<lfsu_s.size(); i++)
u_s += x_s(lfsu_s,i) * phi_s[i];
// check if a BC type is cached
BoundaryCondition::Type bc;
const auto bc_cached = bc_type_cache.find(it_id);
if(bc_cached!=bc_type_cache.end())
bc = bc_cached->second;
else // query type of boundary condition
bc = boundary.bc(ig.intersection(),it.position(),time);
// set the internal bcType flag
auto bcType = BCType::none;
RF normal_flux;
if (BoundaryCondition::isNeumann(bc))
{
// evaluate flux boundary condition
normal_flux = boundary.j(ig.intersection(),it.position(),time);
bcType = BCType::vonNeumann;
}
else if (BoundaryCondition::isDirichlet(bc))
{
bcType = BCType::dirichlet;
}
else if (BoundaryCondition::isLimitedInFlux(bc))
{
// modified Neumann boundary condition, set flux to 0 where potential is positive
normal_flux = .0;
if(u_s < .0) normal_flux = boundary.j(ig.intersection(),it.position(),time);
bcType = BCType::vonNeumann;
}
else if (BoundaryCondition::isEvaporation(bc))
{
// estimate the dirichlet flux
// evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
const auto& gradphi_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
// transform gradients to real element
const auto jac = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
std::vector<Vector> tgradphi_s(lfsu_s.size());
for (unsigned int i = 0; i<lfsu_s.size(); i++)
jac.mv(gradphi_s[i][0],tgradphi_s[i]);
// compute gradient of u
Vector gradu_s(0.);
for (unsigned int i = 0; i<lfsu_s.size(); i++)
gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
// update with gravity vector
gradu_s -= param.gravity();
// face normal vector
const Domain normal = ig.unitOuterNormal(it.position());
// jump relative to Dirichlet value
const RF g = boundary.g(ig.intersection(),it.position(),time);
const RF jump = u_s - g;
const Domain p_global_s = ig.inside().geometry().global(p_local_s);
// conductivity at saturation
const RF satCond_s = param.condRefToMiller(param.K(p_global_s), p_global_s);
const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
// compute numerical flux estimation
const RF numFlux = satCond_s * ( - (gradu_s * normal) + penalty * jump);
// switches between Neumann and Dirichlet bc
// this choice is kept constant for one time step
normal_flux = boundary.j(ig.intersection(),it.position(),time);
if(numFlux < normal_flux) {
bcType = BCType::dirichlet;
bc_type_cache.emplace(it_id,BoundaryCondition::Dirichlet);
}
else {
bcType = BCType::vonNeumann;
bc_type_cache.emplace(it_id,BoundaryCondition::Neumann);
}
}
else if (BoundaryCondition::isOther(bc))
bcType = BCType::none;
// update residual according to byType flag
if (bcType == BCType::vonNeumann)
{
// update jacobian (flux term)
for (unsigned int i = 0; i<lfsv_s.size(); ++i) {
for (unsigned int j = 0; j<lfsu_s.size(); ++j) {
mat_ss.accumulate(lfsv_s,i,lfsu_s,j, 0.0);
}
}
continue;
}
else if (bcType == BCType::dirichlet)
{
// evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
const auto& gradphi_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
// transform gradients to real element
const auto jac = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
std::vector<Vector> tgradphi_s(lfsu_s.size());
for (unsigned int i = 0; i<lfsu_s.size(); i++)
jac.mv(gradphi_s[i][0],tgradphi_s[i]);
// compute gradient of u
Vector gradu_s(0.);
for (unsigned int i = 0; i<lfsu_s.size(); i++)
gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
// update with gravity vector
gradu_s -= param.gravity();
// face normal vector
const Domain normal = ig.unitOuterNormal(it.position());
// jump relative to Dirichlet value
const RF g = boundary.g(ig.intersection(),it.position(),time);
const RF jump = u_s - g;
const Domain p_global_s = ig.inside().geometry().global(p_local_s);
// conductivity at saturation
const RF satCond_s = param.condRefToMiller(param.K(p_global_s), p_global_s);
const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
// compute numerical flux
const RF numFlux = satCond_s * ( - (gradu_s * normal) + penalty * jump);
RF head_upwind;
if (numFlux > 0)
head_upwind = param.headMillerToRef(u_s, p_global_s);
else
head_upwind = g; // Miller here?!
// compute saturation from matrix head
const RF saturation_upwind = param.saturation(head_upwind, p_global_s);
// compute relative conductivity
const RF relCond_upwind = param.condFactor(saturation_upwind, p_global_s);
// derivatives of relCond to s
RF rel_cond_upwind_deriv_s;
if (numFlux > 0){
// relCond depends on s
rel_cond_upwind_deriv_s = param.cond_factor_deriv(saturation_upwind, p_global_s);
rel_cond_upwind_deriv_s *= param.saturation_deriv(head_upwind,p_global_s);
}
else {
// relCond depends on boundary condition
rel_cond_upwind_deriv_s = 0.0;
}
for (unsigned int i = 0; i<lfsv_s.size(); ++i) {
for (unsigned int j = 0; j<lfsu_s.size(); ++j) {
// update jacobian (flux term)
mat_ss.accumulate(lfsv_s,i,lfsu_s,j, satCond_s * ( - ( tgradphi_s[j] * normal ) + penalty * phi_s[j] ) * phi_s[i] * factor );
// update jacobian (symmetry term)
mat_ss.accumulate(lfsv_s,i,lfsu_s,j, theta * ( tgradphi_s[i] * normal ) * satCond_s * phi_s[j] * ( rel_cond_upwind_deriv_s + 1.0 ) * factor );
}
}
continue;
}
else if (bcType == BCType::none)
continue;
} // quadrature rule
} // jacobian_boundary
/**
* @brief Volume integral depending only on test functions
*/
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment