Commit e06c3c58 authored by Michel Henry's avatar Michel Henry
Browse files

periodic singular matrix

parent 0991333a
Pipeline #8319 passed with stages
in 4 minutes and 24 seconds
......@@ -1914,6 +1914,10 @@ double *fluid_problem_element_size(const FluidProblem *p)
{
return p->element_size;
}
int *fluid_problem_parent_tag(const FluidProblem *p)
{
return p->mesh->parent_tag;
}
void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements, int n_boundaries, int *boundaries, int *boundary_tags,int n_physicals, char **physicals)
{
......
......@@ -2,22 +2,22 @@
* MigFlow - Copyright (C) <2010-2018>
* <Universite catholique de Louvain (UCL), Belgium
* Universite de Montpellier, France>
*
*
* List of the contributors to the development of MigFlow: see AUTHORS file.
* Description and complete License: see LICENSE file.
*
* This program (MigFlow) is free software:
* you can redistribute it and/or modify it under the terms of the GNU Lesser General
*
* This program (MigFlow) is free software:
* you can redistribute it and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation, either version
* 3 of the License, or (at your option) any later version.
*
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program (see COPYING and COPYING.LESSER files). If not,
* along with this program (see COPYING and COPYING.LESSER files). If not,
* see <http://www.gnu.org/licenses/>.
*/
......@@ -122,10 +122,10 @@ double *fluid_problem_old_porosity(const FluidProblem *p);
double *fluid_problem_porosity(const FluidProblem *p);
double *fluid_problem_concentration_dg(const FluidProblem *p);
double *fluid_problem_element_size(const FluidProblem *p);
void fluid_problem_advance_concentration(FluidProblem *problem, double dt);
// mesh data (for i/o)
int fluid_problem_n_nodes(const FluidProblem *p);
double *fluid_problem_coordinates(const FluidProblem *p);
......@@ -135,6 +135,7 @@ void fluid_problem_set_elements(FluidProblem *p, int n_nodes, double *x, int n_e
int fluid_problem_n_mesh_boundaries(const FluidProblem *p);
void fluid_problem_mesh_boundary_info(const FluidProblem *p, int bid, char **bname, int *bsize);
void fluid_problem_mesh_boundary_interfaces(const FluidProblem *p, int bid, int *binterfaces);
int *fluid_problem_parent_tag(const FluidProblem *p);
void fluid_problem_set_stab_param(FluidProblem *p, double stab_param);
......
......@@ -344,22 +344,22 @@ class FluidProblem :
rhs = np.zeros_like(solution_old)
self._lib.fluid_problem_reset_boundary_force(self._fp)
self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt))
# Create mapping
map = self.
periodic_map = self._get_matrix("parent_tag", self.n_nodes(), 1, c_int).flatten()
if self.sys is None :
self.sys = LinearSystem(map[self.elements()],self.n_fields(),self.solver_options)
#self.sys = LinearSystem(self.elements(),self.n_fields(),self.solver_options
)
self.sys = LinearSystem(periodic_map[self.elements()],self.n_fields(),self.solver_options)
# self.sys = LinearSystem(self.elements(),self.n_fields(),self.solver_options
# )
sys = self.sys
localv = np.ndarray([sys.localsize*self.n_elements()])
localm = np.ndarray([sys.localsize**2*self.n_elements()])
self._lib.fluid_problem_assemble_system(self._fp,_np2c(rhs),_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
sys.local_to_global(localv,localm,rhs)
# Periodic Conditions
if check_residual_norm > 0 :
norm0 = np.linalg.norm(rhs)
self.solution()[:] -= sys.solve(rhs)[map]
#self.solution()[:] -= sys.solve(rhs)
self.solution()[:] -= sys.solve(rhs)[periodic_map] # Matrice singuliere : logique car on mets plusieurs fois la meme ligne
# self.solution()[:] -= sys.solve(rhs)
if check_residual_norm > 0 :
sys.zero_matrix()
rhs[:] = 0
......
......@@ -4,16 +4,16 @@ y = 0;
lc = 0.05;
Point(1) = {0,0,0,lc};
Point(2) = {1,0,0,lc};
Point(3) = {1,1,0,lc};
Point(4) = {0,1,0,lc};
Point(2) = {L,0,0,lc};
Point(3) = {L,H,0,lc};
Point(4) = {0,H,0,lc};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {4,3};
Line(4) = {1,4};
Line Loop(1) = {1,2,-3,-4};
Plane Surface(1) = {1};
Periodic Curve {1,2} = {3,4};
Periodic Curve {1} = {3};
Physical Line("Left") = {4};
Physical Line("Right") = {2};
......
......@@ -48,11 +48,11 @@ g = np.array([0,0]) # gravity
rho = 1000 # fluid density
nu = 1e-3 # kinematic viscosity
mu = nu*rho # dynamic viscosity
tEnd = 10000 # final time
tEnd = 100 # final time
#numerical parameters
lcmin = .1 # mesh size
dt = 10 # time step
dt = 0.1 # time step
alpha = 1e-4 # stabilization coefficient
epsilon = alpha*lcmin**2 /nu # stabilization parametre
print('epsilon',epsilon)
......@@ -64,10 +64,10 @@ outf = 1 # number of iterations between ou
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
fluid = mbfluid.FluidProblem(2,g,nu*rho,rho)
fluid.load_msh("mesh.msh")
fluid.set_open_boundary("Left",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
fluid.set_wall_boundary("Bottom",velocity=[0,0], pressure=0)
fluid.set_wall_boundary("Top",velocity=[0,0], pressure=0)
fluid.set_open_boundary("Right",velocity=[lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]),0])
fluid.set_open_boundary("Right",pressure = 0)
fluid.set_open_boundary("Left", pressure = 1)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Top")
ii = 0
......
Supports Markdown
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