Skip to content
Snippets Groups Projects
Commit 36c559df authored by Fanny CHOPOT's avatar Fanny CHOPOT
Browse files

modifs de CL et k pour situation tube de choc avec diffusion nulle au debut

parent c239090d
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/main.cpp
+ 15
− 12
View file @ 36c559df
...@@ -138,7 +138,7 @@ int main(int argc, char *argv[]) ...@@ -138,7 +138,7 @@ int main(int argc, char *argv[])
const Kokkos::View<const double*> Vj = mesh_data.Vj(); const Kokkos::View<const double*> Vj = mesh_data.Vj();
const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr(); const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
const double tmax=0.7; const double tmax=1.5;
double t=0; double t=0;
int itermax=std::numeric_limits<int>::max(); int itermax=std::numeric_limits<int>::max();
...@@ -170,7 +170,7 @@ int main(int argc, char *argv[]) ...@@ -170,7 +170,7 @@ int main(int argc, char *argv[])
// ETAPE 2 DU SPLITTING - DIFFUSION // ETAPE 2 DU SPLITTING - DIFFUSION
double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj); double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj);
if (dt_euler <= dt_diff) { if (dt_euler <= dt_diff) {
dt_diff = dt_euler; dt_diff = dt_euler;
...@@ -179,7 +179,7 @@ int main(int argc, char *argv[]) ...@@ -179,7 +179,7 @@ int main(int argc, char *argv[])
double t_diff = t-dt_euler; double t_diff = t-dt_euler;
while (t > t_diff) { while (t > t_diff) {
finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns); finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns);
dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj); dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj);
if (t_diff+dt_diff > t) { if (t_diff+dt_diff > t) {
dt_diff = t-t_diff; dt_diff = t-t_diff;
} }
...@@ -208,7 +208,7 @@ int main(int argc, char *argv[]) ...@@ -208,7 +208,7 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n"; << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
/*
double error1 = 0.; double error1 = 0.;
error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax); error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax);
...@@ -226,7 +226,7 @@ int main(int argc, char *argv[]) ...@@ -226,7 +226,7 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset
<< ": " << rang::fgB::green << error << rang::fg::reset << " \n"; << ": " << rang::fgB::green << error << rang::fg::reset << " \n";
*/
/* /*
double error3 = 0.; double error3 = 0.;
error3 = finite_volumes_diffusion.error_L2_E(unknowns); error3 = finite_volumes_diffusion.error_L2_E(unknowns);
...@@ -251,11 +251,13 @@ int main(int argc, char *argv[]) ...@@ -251,11 +251,13 @@ int main(int argc, char *argv[])
{ // gnuplot output for density { // gnuplot output for density
const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> rhoj = unknowns.rhoj(); const Kokkos::View<const double*> rhoj = unknowns.rhoj();
double h = std::sqrt(1. - (tmax*tmax)/(50./9.)); //double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
std::ofstream fout("rho"); std::ofstream fout("rho1.5");
fout.precision(15); fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) { for (size_t j=0; j<mesh.numberOfCells(); ++j) {
fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n'; fout << xj[j][0] << ' ' << rhoj[j] << '\n';
// Kidder
// fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n';
} }
} }
...@@ -268,8 +270,8 @@ int main(int argc, char *argv[]) ...@@ -268,8 +270,8 @@ int main(int argc, char *argv[])
for (size_t j=0; j<mesh.numberOfCells(); ++j) { for (size_t j=0; j<mesh.numberOfCells(); ++j) {
//fout << xj[j][0] << ' ' << uj[j][0] << ' ' << std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2) <<'\n'; //cas k constant //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2) <<'\n'; //cas k constant
//fout << xj[j][0] << ' ' << uj[j][0] << ' ' << std::sin(pi*xj[j][0])*std::exp(-0.2) <<'\n'; // cas k non constant //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << std::sin(pi*xj[j][0])*std::exp(-0.2) <<'\n'; // cas k non constant
//fout << xj[j][0] << ' ' << uj[j][0] << '\n'; fout << xj[j][0] << ' ' << uj[j][0] << '\n';
fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n';
} }
} }
...@@ -277,13 +279,14 @@ int main(int argc, char *argv[]) ...@@ -277,13 +279,14 @@ int main(int argc, char *argv[])
const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> Ej = unknowns.Ej(); const Kokkos::View<const double*> Ej = unknowns.Ej();
//double pi = 4.*std::atan(1.); //double pi = 4.*std::atan(1.);
double h = std::sqrt(1. - (tmax*tmax)/(50./9.)); //double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
std::ofstream fout("E"); std::ofstream fout("E");
fout.precision(15); fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) { for (size_t j=0; j<mesh.numberOfCells(); ++j) {
//fout << xj[j][0] << ' ' << Ej[j] << ' ' << (-(std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0]))+(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])))*0.5*(std::exp(-4.*pi*pi*0.2)-1.) + 2. <<'\n'; // cas k constant //fout << xj[j][0] << ' ' << Ej[j] << ' ' << (-(std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0]))+(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])))*0.5*(std::exp(-4.*pi*pi*0.2)-1.) + 2. <<'\n'; // cas k constant
//fout << xj[j][0] << ' ' << Ej[j] << ' ' << ((xj[j][0]*pi*pi*0.5)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(xj[j][0]*pi)*std::cos(pi*xj[j][0])) - pi*0.5*std::sin(pi*xj[j][0])*std::cos(pi*xj[j][0]))*(std::exp(-2.*0.2)-1.) + 2. <<'\n' ; // cas k non constant //fout << xj[j][0] << ' ' << Ej[j] << ' ' << ((xj[j][0]*pi*pi*0.5)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(xj[j][0]*pi)*std::cos(pi*xj[j][0])) - pi*0.5*std::sin(pi*xj[j][0])*std::cos(pi*xj[j][0]))*(std::exp(-2.*0.2)-1.) + 2. <<'\n' ; // cas k non constant
fout << xj[j][0] << ' ' << Ej[j] << ' ' << (std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h) + (-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*(-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*0.5 << '\n'; //fout << xj[j][0] << ' ' << Ej[j] << ' ' << (std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h) + (-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*(-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*0.5 << '\n'; // kidder
fout << xj[j][0] << ' ' << Ej[j] << '\n';
} }
} }
......
src/mesh/Mesh.hpp
+ 2
− 2
View file @ 36c559df
...@@ -80,8 +80,8 @@ public: ...@@ -80,8 +80,8 @@ public:
m_x0("x0", 1), m_x0("x0", 1),
m_xmax("xmax", 1) m_xmax("xmax", 1)
{ {
double a = 0.; double a = 0.1;
double b = 10.; double b = 1.;
m_x0[0][0] = a; m_x0[0][0] = a;
m_xmax[0][0] = b; m_xmax[0][0] = b;
const double delta_x = (b-a)/connectivity.numberOfCells(); const double delta_x = (b-a)/connectivity.numberOfCells();
......
src/scheme/AcousticSolver.hpp
+ 9
− 4
View file @ 36c559df
...@@ -181,7 +181,7 @@ private: ...@@ -181,7 +181,7 @@ private:
} }
*/ */
Kokkos::View<Rd*> Kokkos::View<Rd*> // calcule u_r (vitesse au sommet du maillage et flux de vitesse)
computeUr(const Kokkos::View<const Rdd*>& Ar, computeUr(const Kokkos::View<const Rdd*>& Ar,
const Kokkos::View<const Rd*>& br, const Kokkos::View<const Rd*>& br,
const double& t) { const double& t) {
...@@ -197,14 +197,19 @@ private: ...@@ -197,14 +197,19 @@ private:
m_ur[r]=invAr(r)*br(r); m_ur[r]=invAr(r)*br(r);
}); });
m_ur[0]=zero;
m_ur[m_mesh.numberOfNodes()-1]=zero;
// Kidder
// Conditions aux limites dependant du temps // Conditions aux limites dependant du temps
//m_ur[0]=(-t/((50./9.)-t*t))*xr[0]; //m_ur[0]=(-t/((50./9.)-t*t))*xr[0];
//m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1]; //m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1];
// R(t) = x*h(t) a la place de x(t) // R(t) = x*h(t) a la place de x(t)
double h = std::sqrt(1. - (t*t)/(50./9.)); //double h = std::sqrt(1. - (t*t)/(50./9.));
m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0]; //m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0];
m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0]; //m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0];
return m_ur; return m_ur;
} }
......
src/scheme/FiniteVolumesDiffusion.hpp
+ 20
− 17
View file @ 36c559df
...@@ -115,7 +115,7 @@ private: ...@@ -115,7 +115,7 @@ private:
}); });
// Conditions aux bords // Conditions aux bords
/*
int cell_here = face_cells(0,0); int cell_here = face_cells(0,0);
int local_face_number_in_cell = face_cell_local_face(0,0); int local_face_number_in_cell = face_cell_local_face(0,0);
m_Fl(0) = -(kL(0) + kj(cell_here))*(1./(2*Vl(0)))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell))); m_Fl(0) = -(kL(0) + kj(cell_here))*(1./(2*Vl(0)))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell)));
...@@ -124,7 +124,7 @@ private: ...@@ -124,7 +124,7 @@ private:
local_face_number_in_cell = face_cell_local_face(m_mesh.numberOfFaces()-1,0); local_face_number_in_cell = face_cell_local_face(m_mesh.numberOfFaces()-1,0);
m_Fl(m_mesh.numberOfFaces()-1) = -(kR(0) + kj(cell_here))*(1/(2.*Vl(m_mesh.numberOfFaces()-1)))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell))); m_Fl(m_mesh.numberOfFaces()-1) = -(kR(0) + kj(cell_here))*(1/(2.*Vl(m_mesh.numberOfFaces()-1)))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
//m_Fl(m_mesh.numberOfFaces()-1) = -xr[m_mesh.numberOfNodes()-1][0]*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell))); //m_Fl(m_mesh.numberOfFaces()-1) = -xr[m_mesh.numberOfNodes()-1][0]*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
*/
// k = 0.5 // k = 0.5
//m_Fl(0,0) = -(t/((50./9.)-t*t))*0.5; //m_Fl(0,0) = -(t/((50./9.)-t*t))*0.5;
...@@ -134,9 +134,9 @@ private: ...@@ -134,9 +134,9 @@ private:
//m_Fl(0,0) = -(t/((50./9.)-t*t))*xr[0][0]; //m_Fl(0,0) = -(t/((50./9.)-t*t))*xr[0][0];
//m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*xr[m_mesh.numberOfFaces()-1][0]; //m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*xr[m_mesh.numberOfFaces()-1][0];
double h = std::sqrt(1. - (t*t)/(50./9.)); //double h = std::sqrt(1. - (t*t)/(50./9.));
m_Fl(0,0) = -(t/((50./9.)-t*t))*h*x0[0][0]; //m_Fl(0,0) = -(t/((50./9.)-t*t))*h*x0[0][0];
m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*h*xmax[0][0]; //m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*h*xmax[0][0];
return m_Fl ; return m_Fl ;
} }
...@@ -176,15 +176,15 @@ private: ...@@ -176,15 +176,15 @@ private:
}); });
// Conditions aux bords // Conditions aux bords
// m_Gl(0) = Fl(0)*uL(0); m_Gl(0) = Fl(0)*uL(0);
//m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0); m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
//m_Gl(0) = -(t/((50./9.)-t*t))*Fl(0,0)*xr(0); //m_Gl(0) = -(t/((50./9.)-t*t))*Fl(0,0)*xr(0);
//m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*Fl(m_mesh.numberOfFaces()-1,0)*xr(m_mesh.numberOfFaces()-1); //m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*Fl(m_mesh.numberOfFaces()-1,0)*xr(m_mesh.numberOfFaces()-1);
double h = std::sqrt(1. - (t*t)/(50./9.)); //double h = std::sqrt(1. - (t*t)/(50./9.));
m_Gl(0) = -(t/((50./9.)-t*t))*h*Fl(0,0)*x0(0); //m_Gl(0) = -(t/((50./9.)-t*t))*h*Fl(0,0)*x0(0);
m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*h*Fl(m_mesh.numberOfFaces()-1,0)*xmax(0); //m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*h*Fl(m_mesh.numberOfFaces()-1,0)*xmax(0);
return m_Gl ; return m_Gl ;
...@@ -244,12 +244,11 @@ public: ...@@ -244,12 +244,11 @@ public:
// Utilise la reduction definie dans la structure ReduceMin. // Utilise la reduction definie dans la structure ReduceMin.
KOKKOS_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
double diffusion_dt(const Kokkos::View<const double*>& rhoj, double diffusion_dt(const Kokkos::View<const double*>& rhoj,
const Kokkos::View<const double*>& kj) const { const Kokkos::View<const double*>& kj,
const Kokkos::View<const double*>& cj) const {
Kokkos::View<double*> dt_j("dt_j", m_mesh.numberOfCells()); Kokkos::View<double*> dt_j("dt_j", m_mesh.numberOfCells());
const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
const Kokkos::View<const double*>& Vl = m_mesh_data.Vl(); const Kokkos::View<const double*>& Vl = m_mesh_data.Vl();
const Kokkos::View<const double*>& Vj = m_mesh_data.Vj(); const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
...@@ -276,8 +275,12 @@ public: ...@@ -276,8 +275,12 @@ public:
sum += kj(cell_nodes(j,m)); sum += kj(cell_nodes(j,m));
} }
if (sum == 0.) {
dt_j[j] = Vj[j]/cj[j];
} else {
// dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./(2.*kj(j) + sum))*minVl; // dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./(2.*kj(j) + sum))*minVl;
dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./sum)*minVl; dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./sum)*minVl;
}
}); });
...@@ -343,8 +346,8 @@ public: ...@@ -343,8 +346,8 @@ public:
//Ej[j] -= (dt*inv_mj[j])*Vj(j)*((0.5*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((0.5*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
// ajout second membre pour kidder (k = x) // ajout second membre pour kidder (k = x)
uj[j][0] += (dt*inv_mj[j])*Vj(j)*(t/((50./9.)-t*t)); //uj[j][0] += (dt*inv_mj[j])*Vj(j)*(t/((50./9.)-t*t));
Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
}); });
......
src/scheme/FiniteVolumesEulerUnknowns.hpp
+ 16
− 12
View file @ 36c559df
...@@ -191,24 +191,23 @@ public: ...@@ -191,24 +191,23 @@ public:
const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
/*
if (xj[j][0]<0.5) { if (xj[j][0]<0.5) {
m_rhoj[j]=1; m_rhoj[j]=1;
} else { } else {
m_rhoj[j]=0.125; m_rhoj[j]=0.125;
} }
*/ //Kidder
m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/100.); //m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/100.);
}); });
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
/*if (xj[j][0]<0.5) { if (xj[j][0]<0.5) {
m_pj[j]=1; m_pj[j]=1;
} else { } else {
m_pj[j]=0.1; m_pj[j]=0.1;
} }
*/
m_pj[j] = 2.*std::pow(m_rhoj[j],3); //m_pj[j] = 2.*std::pow(m_rhoj[j],3);
}); });
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
...@@ -216,8 +215,8 @@ public: ...@@ -216,8 +215,8 @@ public:
}); });
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
//m_gammaj[j] = 1.4; m_gammaj[j] = 1.4;
m_gammaj[j] = 3.; //m_gammaj[j] = 3.;
}); });
BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj); BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
...@@ -237,16 +236,21 @@ public: ...@@ -237,16 +236,21 @@ public:
}); });
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_kj[j] = xj[j][0]; //m_kj[j] = xj[j][0];
//m_kj[j] = 0.5; //m_kj[j] = 0.5;
if (xj[j][0]<0.7) {
m_kj[j]=0.;
} else {
m_kj[j]=0.5;
}
}); });
// Conditions aux bords de Dirichlet sur u et k // Conditions aux bords de Dirichlet sur u et k
m_uL[0] = zero; m_uL[0] = zero;
m_uR[0] = zero; m_uR[0] = zero;
m_kL[0] = xj[0][0]; m_kL[0] = 0.01;
m_kR[0] = xj[m_mesh.numberOfCells()-1][0]; m_kR[0] = 0.5;
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment