The TS-GitLab will have to shut down towards the end of the year — please think about migrating your projects to GitLab.com or GitHub.
(This is still a very early message, meant to keep you informed. There will be more urgent ones in the future.)

Commit e8a76d10 authored by Lukas Riedel's avatar Lukas Riedel

fixes to jacobian_skeleton. Works really badly.

* multiply 'numFlux' to derivative of 'relCond'
* upwind position queries are always on '_s' side
parent e9fda0a2
......@@ -63,14 +63,16 @@ void Simulation<Traits>::operator_setup ()
go1 = std::make_unique<GO1>(*gfs,*cc,*gfs,*cc,*tlop,mbe_tlop);
igo = std::make_unique<IGO>(*go0,*go1);
slop->setTime(0.0);
// print jacobian
std::ofstream file;
file.open("jacobian.txt");
typename GO1::Traits::Jacobian jac(*go1);
jac = 0.0;
go1->jacobian(*this->unew,jac);
auto jac_native = Dune::PDELab::Backend::native(jac);
Dune::printmatrix(file,jac_native,"jacobian","");
// std::ofstream file;
// file.open("jacobian.txt");
// typename GO0::Traits::Jacobian jac(*go0);
// jac = 0.0;
// go0->jacobian(*this->unew,jac);
// auto jac_native = Dune::PDELab::Backend::native(jac);
// Dune::printmatrix(file,jac_native,"jacobian","");
// --- Solvers ---
ls = std::make_unique<LS>(5000,0);
......
......@@ -56,7 +56,7 @@ template<typename Traits, typename Parameter, typename Boundary, typename Source
class RichardsDGSpatialOperator
:
// 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::NumericalJacobianSkeleton<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,
......@@ -129,7 +129,7 @@ public:
int intorderadd_ = 2, int quadrature_factor_ = 2)
:
// 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::NumericalJacobianSkeleton<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_),
......@@ -553,12 +553,15 @@ public:
rel_cond_upwind_deriv_n = 0.0;
}
else {
// relCond depends on s
// relCond depends on n
rel_cond_upwind_deriv_s = 0.0;
rel_cond_upwind_deriv_n = param.cond_factor_deriv(saturation_upwind, p_global_n);
rel_cond_upwind_deriv_n *= param.saturation_deriv(head_upwind,p_global_n);
rel_cond_upwind_deriv_n = param.cond_factor_deriv(saturation_upwind, p_global_s);
rel_cond_upwind_deriv_n *= param.saturation_deriv(head_upwind,p_global_s);
}
// std::cout << "rel_cond_upwind_deriv_s: " << rel_cond_upwind_deriv_s << std::endl;
// std::cout << "rel_cond_upwind_deriv_n: " << rel_cond_upwind_deriv_n << std::endl;
// update jacobian (flux term)
// diffusion term
// consistency term
......@@ -569,7 +572,7 @@ public:
for (unsigned int j = 0; j<lfsu_s.size(); ++j){
// differentiate flux term
mat_ss.accumulate(lfsv_s,i,lfsu_s,j,
rel_cond_upwind_deriv_s * phi_s[j] * phi_s[i] * factor );
rel_cond_upwind_deriv_s * phi_s[j] * numFlux * phi_s[i] * factor );
mat_ss.accumulate(lfsv_s,i,lfsu_s,j,
relCond_upwind
* ( omega_s * satCond_s
......@@ -590,7 +593,7 @@ public:
for (unsigned int j = 0; j<lfsu_n.size(); ++j){
// differentiate flux term
mat_nn.accumulate(lfsv_n,i,lfsu_n,j,
- rel_cond_upwind_deriv_n * phi_n[j] * phi_n[i] * factor );
- rel_cond_upwind_deriv_n * phi_n[j] * numFlux * phi_n[i] * factor );
mat_nn.accumulate(lfsv_n,i,lfsu_n,j,
- relCond_upwind
* ( omega_n * satCond_n
......@@ -611,7 +614,7 @@ public:
for (unsigned int j = 0; j<lfsu_n.size(); ++j){
// differentiate flux term
mat_sn.accumulate(lfsv_s,i,lfsu_n,j,
rel_cond_upwind_deriv_n * phi_n[j] * phi_s[i] * factor );
rel_cond_upwind_deriv_n * phi_n[j] * numFlux * phi_s[i] * factor );
mat_sn.accumulate(lfsv_s,i,lfsu_n,j,
relCond_upwind
* ( omega_n * satCond_n
......@@ -632,7 +635,7 @@ public:
for (unsigned int j = 0; j<lfsu_s.size(); ++j){
// differentiate flux term
mat_ns.accumulate(lfsv_n,i,lfsu_s,j,
- rel_cond_upwind_deriv_s * phi_s[j] * phi_n[i] * factor );
- rel_cond_upwind_deriv_s * phi_s[j] * numFlux * phi_n[i] * factor );
mat_ns.accumulate(lfsv_n,i,lfsu_s,j,
- relCond_upwind
* ( omega_s * satCond_s
......@@ -648,28 +651,11 @@ public:
);
}
}
/*
// update residual (flux term)
// diffusion term
// consistency term
// + penalty term
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, relCond_upwind * numFlux * phi_s[i] * factor);
for (unsigned int i = 0; i<lfsv_n.size(); i++)
r_n.accumulate(lfsv_n,i, - relCond_upwind * numFlux * phi_n[i] * factor);
// update residual (symmetry term)
// (non-)symmetric IP term
// symmetry term
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, theta * omega_s * (tgradphi_s[i] * normal) * relCond_upwind * satCond_s * jump * factor);
for (unsigned int i = 0; i<lfsv_n.size(); i++)
r_n.accumulate(lfsv_n,i, theta * omega_n * (tgradphi_n[i] * normal) * relCond_upwind * satCond_n * jump * factor);
*/
} // quadrature rule
} // jacobian_skeleton
/**
* @brief Boundary integral depending on test and ansatz functions
* \todo there is additional code for the limited influx BC, which is currently commented out
......
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