Commit efbd417a authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

no double division by c in concentration eq

parent 8d611dca
Pipeline #5059 passed with stage
in 27 seconds
......@@ -53,20 +53,16 @@ void fluid_problem_advance_concentration(FluidProblem *problem, double dt) {
double phi[N_SF];
shape_functions(QP[iqp],phi);
double u[D] = {0};
double c = 0;
double cold = 0;
double a = 0;
for (int i = 0; i < N_SF; ++i) {
a += phi[i]*problem->concentration[iel*N_N+i];
c += phi[i]*problem->porosity[el[i]];
cold += phi[i]*problem->oldporosity[el[i]];
for (int id = 0; id < D; ++id) {
u[id] += phi[i]*velocity[el[i]*D+id];
}
}
double r = 0;
for (int id = 0; id < D; ++id) {
r -= u[id]*da[id]/c;
r -= u[id]*da[id];
}
const double jw = QW[iqp]*detj;
for (int i = 0; i < N_SF; ++i) {
......@@ -94,18 +90,16 @@ void fluid_problem_advance_concentration(FluidProblem *problem, double dt) {
double a0 = 0, a1 = 0, un = 0;
double phil[N_LQP];
l_shape_functions(LQP[iqp],phil);
double c = 0;
for (int i = 0; i < N_LSF; ++i) {
a0 += problem->concentration[iel0*N_N+cl0[i]]*phil[i];
a1 += problem->concentration[iel1*N_N+cl1[i]]*phil[i];
c += problem->porosity[el0[cl0[i]]]*phil[i];
for (int id = 0; id < D; ++id) {
un += velocity[el0[cl0[i]]*D+id]*phil[i]*n[id];
}
}
const double a = un > 0 ? a0 : a1;
const double r0 = -un*(a-a0)/c;
const double r1 = un*(a-a1)/c;
const double r0 = -un*(a-a0);
const double r1 = un*(a-a1);
const double jw = LQW[iqp]*detbnd;
for (int i = 0; i < N_LSF; ++i){
rhs[iel0*N_N+cl0[i]] += phil[i]*r0*jw;
......@@ -141,10 +135,8 @@ void fluid_problem_advance_concentration(FluidProblem *problem, double dt) {
double phil[N_LQP];
l_shape_functions(LQP[iqp],phil);
double a = 0;
double c = 0;
for (int i = 0; i < N_LSF; ++i) {
a += problem->concentration[iel*N_N+cl[i]]*phil[i];
c += problem->porosity[el[cl[i]]]*phil[i];
}
double un = 0;
for (int id = 0; id< D; ++id) {
......
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