Skip to content

GitLab

Commit e90ed53f authored by Nathan Coppin's avatar Nathan Coppin
Browse files

Merge branch 'master' into nathan

parents 29e0f581 13b06185
Pipeline #7757 passed with stages
in 3 minutes and 13 seconds
Hide whitespace changes
Inline Side-by-side
.gitignore
View file @ e90ed53f
output
*.so
*.core
*.msh
*.pos
__pycache__
......@@ -8,3 +6,11 @@ __pycache__
*.aux
*.log
*.pdf
*.vtu
*.vtp
*.vtm
*.pvd
*.npy
build
.vscode
wip
.gitlab-ci.yml
View file @ e90ed53f
......@@ -19,13 +19,60 @@
# along with this program (see COPYING and COPYING.LESSER files). If not,
# see <http://www.gnu.org/licenses/>.
mbtests :
image : immc/migflow-valid:v0.8
mfbuild :
image : immc/migflow-build:v0.6
stage : build
script:
# linux
- mkdir build
- cd build/
- cmake .. -DENABLE_PETSC=1 -DPETSC_DIR=/usr/lib/petscdir/petsc3.9/ -DPETSC_ARCH=x86_64-linux-gnu-real
- PPATH=$(pwd)
- cmake ..
- make -j4
- cd ../validation
- PYTHONPATH=$PPATH:$PYTHONPATH python3 mbtests.py
- python3 setup.py bdist_wheel --plat-name=manylinux1_x86_64 -d ../dist
- cd ..
# windows
- mkdir build-mingw
- cd build-mingw
- cmake .. -DCMAKE_TOOLCHAIN_FILE=/cmake-mingw
- make -j4
- mv */*.dll */*.dll.a migflow
- python3 setup.py bdist_wheel --plat-name=win_amd64 -d ../dist
- cd ..
# OSX
- mkdir build-osxcross
- cd build-osxcross
- $(osxcross-conf) && cmake .. -DCMAKE_TOOLCHAIN_FILE=/osxcross/target/toolchain.cmake
- make -j4
- python3 setup.py bdist_wheel --plat-name=macosx_10_9_x86_64 -d ../dist
- cd ..
artifacts:
paths:
- dist
expire_in: 1 day
mftests :
image : immc/migflow-valid:v0.12
stage : test
script:
- pip3 install --user dist/*manylinux1_x86_64*
- cd validation
- python3 mbtests.py
mfdeploy-test :
image : immc/migflow-build:v0.6
stage : deploy
rules :
- if: '$CI_COMMIT_TAG =~ /^w-.*$/'
when: always
script:
- twine upload --repository testpypi dist/*
mfdeploy :
image : immc/migflow-build:v0.6
stage : deploy
rules :
- if: '$CI_COMMIT_TAG =~ /^v-.*$/'
when: always
script:
- twine upload dist/*
.vscode/c_cpp_properties.json deleted 100644 → 0
View file @ 29e0f581
{
"env" : {
"defaultIncludePath": [
"${workspaceFolder}",
"${workspaceFolder}/hxt",
"${workspaceFolder}/src",
"${workspaceFolder}/scontact"
],
"myCompilerPath": "/usr/bin/clang"
},
"configurations": [
{
"name": "Linux",
"intelliSenseMode": "clang-x64",
"includePath": [ "${defaultIncludePath}"],
"defines": [ "HXT_HAVE_PETSC","DIMENSION=2" ],
"compilerPath": "/usr/bin/clang",
"cStandard": "c11",
"cppStandard": "c++17",
"browse": {
"path": [ "${workspaceFolder}" ],
"limitSymbolsToIncludedHeaders": true,
"databaseFilename": ""
}
}
],
"version": 4
}
CMakeLists.txt
View file @ e90ed53f
......@@ -21,7 +21,6 @@
cmake_minimum_required(VERSION 3.1 FATAL_ERROR)
Project(MigFlow C)
option(ENABLE_PETSC "Enable PETSc" OFF)
## Set a default build type if none was specified
if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
......@@ -36,18 +35,12 @@ list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/migflow)
if(ENABLE_PETSC)
find_package(PETSc)
if(PETSC_FOUND)
add_definitions("-DHXT_HAVE_PETSC ${PETSC_CFLAGS}")
link_directories("${PETSC_LIBDIR}")
list(APPEND EXTERNAL_LIBRARIES ${PETSC_LIBS})
list(APPEND EXTERNAL_INCLUDES ${PETSC_INC})
endif(PETSC_FOUND)
endif(ENABLE_PETSC)
add_subdirectory(python)
add_subdirectory(tools)
add_subdirectory(scontact)
add_subdirectory(fluid)
configure_file(README.md ${PROJECT_BINARY_DIR}/README.md COPYONLY)
configure_file(LICENSE.txt ${PROJECT_BINARY_DIR}/LICENSE.txt COPYONLY)
configure_file(AUTHORS.txt ${PROJECT_BINARY_DIR}/AUTHORS.txt COPYONLY)
configure_file(COPYING.LESSER.txt ${PROJECT_BINARY_DIR}/COPYING.LESSER.txt COPYONLY)
configure_file(COPYING.txt ${PROJECT_BINARY_DIR}/COPYING.txt COPYONLY)
LICENSE.txt
View file @ e90ed53f
......@@ -27,37 +27,3 @@ The list of the contributors to the development of MigFlow is given in the AUTHO
Website of MigFlow: https://git.immc.ucl.ac.be/fluidparticles/MigFlow
Contact: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
This program has dependencies on other libraries and Website depedencies which are or not under GPL and that are commonly distributed with this program core libraries.
The other dependencies are as follows:
PETSC BSD 2-clause https://www.mcs.anl.gov/petsc/
-----------------------------------------------------------------------------
######
PETSC
######
Copyright (c) 1991-2014, UChicago Argonne, LLC and the PETSc Development Team
All rights reserved.
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice, this
list of conditions and the following disclaimer in the documentation and/or
other materials provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-----------------------------------------------------------------------------
doc/mixed.tex 0 → 100644
View file @ e90ed53f
\documentclass[a4]{paper}
\usepackage[utf8]{inputenc}
\usepackage[french]{babel}
\usepackage[vmargin=2cm]{geometry}
\usepackage{amsmath}
\newcommand{\sumop}[1]{\left[#1\right]_{\Sigma}}
\newcommand{\sumopp}[1]{\left[#1\right]_{\Sigma^{n+1}}}
\newcommand{\dpartial}[2]{\frac{\partial #1}{\partial #2}}
\begin{document}
\subsection*{Continu}
\[
\sumop{a} = \sum_p{a_p}\delta_{x_p}
\]
\begin{align*}
\nabla\cdot(u+v) &= 0\\
\text{Particules}&\\
M_p\frac{\partial v_p}{\partial t} &= -V_p\left.\frac{dP}{dy}\right|_{x_p} + M_p g -D_p+R_p\\
\text{Modèle A}&\\
\rho_f \frac{\partial u}{\partial t} &= -\rho_f \nabla \cdot \frac{uu}{c}-(1-\sumop{V})\frac{dP}{dy}+(1-\sumop{V})g\rho_f + \sumop{D}\\
\text{Modèle B}&\\
\rho_f \frac{\partial u}{\partial t} &= -\rho_f \nabla \cdot \frac{uu}{c}-\frac{dP}{dy}+\sumop{V\frac{dP}{dy}}+(1-\sumop{V})g\rho_f + \sumop{D}
\end{align*}
$\text{Modèle B} + \sumop{\text{Particules}}$
\[
\rho_f \frac{\partial u}{\partial t} + \sumop{M\dpartial vt}= -\rho_f \nabla \cdot \frac{uu}{c}-\frac{dP}{dy}+(1-\sumop{V})g\rho_f +\sumop{Mg}+ \sumop{R}
\]
Ce qui est différent de ce qu'on fait pour le moment car on ne reconstruit pas un $c$ continu.
\subsection*{Discrétisation temporelle}
\begin{align*}
M_p\frac{v_p^{n+1}-v_p^n}{\Delta t} =& -V_p\left.\frac{dP^n}{dy}\right|_{x_p} + M_p g -D^n_p+R^{n+1}_p\\
\nabla\cdot(u^{n+1}+v^{n+1}) =&~ 0 \\
\rho_f\frac{u^{n+1}-u^{n}}{\Delta t} + \sumopp{M\frac{v^{n+1}-v^{n}}{\Delta t}}=&
-\rho_f \nabla \cdot \frac{u^{n}u^{n+1}}{c}-\left(\frac{dP}{dy}^{n+1}-g\rho_f\right) \\
&+ \sumopp{Vg(\rho_p-\rho_f)+R^{n+1}}
\end{align*}
On peut ré-écrire la dernière équation :
\begin{align*}
\rho_f\frac{u^{n+1}-u^{n}}{\Delta t} =&
-\rho_f \nabla \cdot \frac{u^{n}u^{n+1}}{c}-\left(\frac{dP}{dy}^{n+1} -g\rho_f\right)
+ \sumopp{V\left(\frac{dP}{dy}^n-g\rho_f\right)+D^n}
\end{align*}
Version certe peu intuitive mais équivalente à la ligne précédente.
L'équation des particules peut faire apparaître le même terme :
\begin{align*}
M_p\frac{v_p^{n+1}-v_p^n}{\Delta t} =& -V_p\left(\left.\frac{dP^n}{dy}\right|_{x_p} -g\rho_f\right) -D^n_p+ g V_p(\rho_p-\rho_f)+R^{n+1}_p\\
\end{align*}
\subsection*{Questions}
\begin{itemize}
\item laquelle on code ? Normalement, c'est vraiment équivalent.
\item $\Sigma^{n+1}$ ou $\Sigma^{n}$ ? Aucune idée. On commence avec $n+1$.
\item quid de l'équation d'incompressibilité ? Pour le moment on garde la version implémentée avec le $\dpartial ct$.
\item est-ce qu'on peut se débarasser complètement de $c$ ? On est embêté avec les $\delta$ au dénominateur.
\item $D^{n+1}$ à la place de $D^n$ ? à priori non, on veut garder l'équivalence. Utiliser un $D$ implicite dans chaque phase ?
\end{itemize}
%
%\begin{align*}
%&\int_\Omega \sum_p M_p\frac{\partial v_p}{\partial t} \delta(x,x_p) \phi_n(x) d\Omega =\\
%&\int_\Omega \sum_p \left[-V_p\left.\frac{dP}{dy}\right|_{x_p} + M_p g -D_p+R_p\right] \delta(x,x_p) \phi_n(x) d\Omega\\
%\end{align*}
%\begin{align*}
%\end{align*}
%Ce qui est sur la feuille :
%\begin{align*}
%\rho_s \frac{v^{n+1}-v^{n}}{\Delta t} &= -(1-c)^{n+1}\frac{dP}{dy}^n + (1-c)^{n+1} g\rho_s -D^{n}+R^{n+1}\\
%u^{n+1}&=-v^{n+1}\\
%\rho_f \frac{u^{n+1}-u^{n}}{\Delta t} + \rho_s\frac{v^{n+1}-v^{n}}{\Delta t} &=
%-\rho_f \nabla \cdot \frac{u^nu^{n+1}}{c}
%\frac{dP}{dy}^{n+1}+g\bar{\rho} + R^{n+1}
%\end{align*}
%Détail :
%Tous les termes correspondent à leur version intégrée vue par le fluide (c'est comme ça qu'on les calcule), sauf le $\frac{dP}{dy}$.
\end{document}
docker/builder/Dockerfile 0 → 100644
View file @ e90ed53f
FROM ubuntu:20.04
env DEBIAN_FRONTEND=noninteractive
#linux
RUN apt update && apt install -y git python3 make cmake gcc
#windows
RUN apt update && apt install -y mingw-w64 mingw-w64-tools
run echo "\
SET(CMAKE_SYSTEM_NAME Windows)\n\
SET(CMAKE_C_COMPILER x86_64-w64-mingw32-gcc)\n\
SET(CMAKE_CXX_COMPILER x86_64-w64-mingw32-g++)\n\
SET(CMAKE_RC_COMPILER x86_64-w64-mingw32-windres)\n\
SET(CMAKE_Fortran_COMPILER x86_64-w64-mingw32-gfortran)\n\
SET(CMAKE_FIND_ROOT_PATH /mingw64)\n\
SET(CMAKE_FIND_ROOT_PATH_MODE_PROGRAM NEVER)\n\
SET(CMAKE_FIND_ROOT_PATH_MODE_LIBRARY ONLY)\n\
SET(CMAKE_FIND_ROOT_PATH_MODE_INCLUDE ONLY)\n"\
> /cmake-mingw
RUN apt update && apt install -y git python3-setuptools python3-wheel twine
### osxcross
RUN apt-get update && apt-get install -y clang patch libssl-dev liblzma-dev libxml2-dev llvm-dev uuid-dev zlib1g-dev
RUN useradd -ms /bin/bash validator
RUN mkdir -p /osxcross && chown validator /osxcross
USER validator
RUN git clone https://github.com/tpoechtrager/osxcross.git /osxcross
COPY MacOSX10.15.sdk.tar.xz /osxcross/tarballs/MacOSX10.15.sdk.tar.xz
RUN cd /osxcross && SDK_VERSION=10.15 UNATTENDED=1 ./build.sh
ENV PATH=/osxcross/target/bin/:$PATH
ENV OSXCROSS_HOST=x86_64-apple-darwin19
VOLUME ["/etc/gitlab-runner"]
WORKDIR /home/validator
docker/Dockerfile → docker/valid/Dockerfile
View file @ e90ed53f
FROM ubuntu:20.04
RUN apt update
RUN DEBIAN_FRONTEND=noninteractive apt install -y git python3 petsc-dev python3-numpy make gmsh cmake g++ python3-distutils
ENV DEBIAN_FRONTEND=noninteractive
RUN apt update && apt install -y git python3 python3-numpy gmsh python3-scipy python3-petsc4py python3-pip
VOLUME ["/etc/gitlab-runner"]
......
fluid/CMakeLists.txt
View file @ e90ed53f
......@@ -25,21 +25,12 @@ set(SRC
fluid_problem_concentration.c
mesh.c
mesh_find.c
hxt_linear_system.c
hxt_linear_system_lu.c
hxt_message.c
../tools/quadtree.c
)
if(PETSC_FOUND)
add_definitions("-DHXT_HAVE_PETSC")
list(APPEND SRC hxt_linear_system_petsc.c)
endif(PETSC_FOUND)
include_directories(. ${CMAKE_SOURCE_DIR}/tools ${EXTERNAL_INCLUDES})
add_library(mbfluid2 SHARED ${SRC})
target_link_libraries(mbfluid2 ${EXTERNAL_LIBRARIES} mbtools2)
target_compile_definitions(mbfluid2 PUBLIC "-DDIMENSION=2")
add_library(mbfluid3 SHARED ${SRC})
target_link_libraries(mbfluid3 ${EXTERNAL_LIBRARIES} mbtools3)
target_compile_definitions(mbfluid3 PUBLIC "-DDIMENSION=3")
fluid/fluid_problem.c
View file @ e90ed53f
......@@ -77,7 +77,7 @@ inline static double node_volume_from_detj(double detj) {
int fluid_problem_n_fields(const FluidProblem *p) {return D+1;}
static void node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
void fluid_problem_node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
const Mesh *mesh = problem->mesh;
const double *porosity = problem->porosity;
......@@ -208,24 +208,22 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
for (int id = 0; id < D; ++id) {
f0[U+id] = c*(pid<0?0:data[pid]-p)*n[id];
f00[(U+id)*n_fields+P] += c*(pid<0?0:-n[id]);
}
if (!(vid<0)) {
double c_du_o_c[D][D];
double sigma = (1+D)/(D*h)*mu*N_N;
for (int iD = 0; iD < D; ++iD) {
for (int jD = 0; jD < D; ++jD)
c_du_o_c[iD][jD] = ds[(U+iD)*D+jD] -u[iD]/c*dc[jD];
}
for (int id = 0; id < D; ++id) {
f0[U+id] += vid<0?0:sigma*(u[id]-data[vid+id]+s_c*n[id]);
f00[(U+id)*n_fields+U+id] += (vid<0?0:sigma);
for (int jd = 0; jd <D; ++jd) {
f0[U+id] -= mu*(c_du_o_c[id][jd]+c_du_o_c[jd][id])*n[jd];
f00[(U+id)*n_fields+U+id] += (vid<0?0:n[id]*n[jd]*sigma) + mu/c*dc[jd]*n[jd];
f00[(U+id)*n_fields+U+jd] += mu/c*dc[id]*n[jd];
f01[(U+id)*n_fields*D+(U+id)*D+jd] -= mu*n[jd];
f01[(U+id)*n_fields*D+(U+jd)*D+id] -= mu*n[jd];
}
}
double c_du_o_c[D][D];
double sigma = problem->ip_factor*(1+D)/(D*h)*mu*N_N;
for (int iD = 0; iD < D; ++iD) {
for (int jD = 0; jD < D; ++jD)
c_du_o_c[iD][jD] = ds[(U+iD)*D+jD] -u[iD]/c*dc[jD];
}
for (int id = 0; id < D; ++id) {
f0[U+id] += vid<0?0:sigma*(u[id]-data[vid+id]+s_c*n[id]);
f00[(U+id)*n_fields+U+id] += (vid<0?0:sigma);
for (int jd = 0; jd <D; ++jd) {
f0[U+id] -= mu*(c_du_o_c[id][jd]+c_du_o_c[jd][id])*n[jd];
f00[(U+id)*n_fields+U+id] += (vid<0?0:n[id]*n[jd]*sigma) + mu/c*dc[jd]*n[jd];
f00[(U+id)*n_fields+U+jd] += mu/c*dc[id]*n[jd];
f01[(U+id)*n_fields*D+(U+id)*D+jd] -= mu*n[jd];
f01[(U+id)*n_fields*D+(U+jd)*D+id] -= mu*n[jd];
}
}
if (problem->advection) {
......@@ -242,7 +240,7 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
}
static double pow2(double a) {return a*a;}
static void compute_stab_parameters(FluidProblem *problem, double dt) {
void fluid_problem_compute_stab_parameters(FluidProblem *problem, double dt) {
const Mesh *mesh = problem->mesh;
problem->element_size = realloc(problem->element_size,sizeof(double)*mesh->n_elements);
......@@ -1117,17 +1115,14 @@ static void fluid_problem_dg_to_cg_grad(FluidProblem *problem, const double *dg,
}
}
static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt)
void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
HXTLinearSystem *lsys = problem->linear_system;
const Mesh *mesh = problem->mesh;
int n_fields = fluid_problem_n_fields(problem);
size_t local_size = N_SF*n_fields;
size_t N_E = mesh->n_elements;
double *all_local_vector = (double*) malloc(N_E*local_size*sizeof(double));
double *all_local_matrix = (double*) malloc(N_E*local_size*local_size*sizeof(double));
for (size_t i = 0; i < N_E*local_size; ++i)
all_local_vector[i] = 0;
for (size_t i = 0; i < N_E*local_size*local_size; ++i)
......@@ -1151,7 +1146,7 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
forced_row[problem->boundaries[iphys].nodes[i]*n_fields+bnd->row] = bnd->field;
}
}
node_force_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
fluid_problem_node_force_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
compute_weak_boundary_conditions(problem, solution_old, dt, all_local_vector, all_local_matrix);
fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
if (problem->n_fluids == 2){
......@@ -1178,15 +1173,11 @@ static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, co
local_vector[inode+ifield*N_SF] = 0;
}
}
hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
}
free(forced_row);
free(all_local_matrix);
free(all_local_vector);
}
static void apply_boundary_conditions(FluidProblem *problem)
void fluid_problem_apply_boundary_conditions(FluidProblem *problem)
{
const Mesh *mesh = problem->mesh;
int n_fields = fluid_problem_n_fields(problem);
......@@ -1217,53 +1208,6 @@ static void apply_boundary_conditions(FluidProblem *problem)
}
}
void fluid_problem_implicit_euler(FluidProblem *problem, double dt, double check_residual_norm)
{
const Mesh *mesh = problem->mesh;
int n_fields = fluid_problem_n_fields(problem);
apply_boundary_conditions(problem);
double *solution_old = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
double *rhs = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
for (int i=0;i<D*problem->mesh->n_boundaries;i++)
problem->boundary_force[i] = 0;
for (int i=0; i<mesh->n_nodes*n_fields; ++i)
solution_old[i] = problem->solution[i];
double *dx = (double*)malloc(sizeof(double)*mesh->n_nodes*n_fields);
double norm0 = 0, norm = 0;
int i;
compute_stab_parameters(problem,dt);
hxtLinearSystemZeroMatrix(problem->linear_system);
for (int i=0; i<mesh->n_nodes*n_fields; ++i){
rhs[i] = 0;
}
fluid_problem_assemble_system(problem, rhs, solution_old, dt);
if (check_residual_norm > 0) {
hxtLinearSystemGetRhsNorm(problem->linear_system,rhs,&norm0);
}
hxtLinearSystemSolve(problem->linear_system, rhs, dx);
for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
problem->solution[j] -= dx[j];
}
if (check_residual_norm > 0) {
hxtLinearSystemZeroMatrix(problem->linear_system);
for (int i=0; i<mesh->n_nodes*n_fields; ++i){
rhs[i] = 0;
}
fluid_problem_assemble_system(problem, rhs, solution_old, dt);
hxtLinearSystemGetRhsNorm(problem->linear_system,rhs,&norm);
printf("residual norm : %g => %g\n", norm0, norm);
if(norm > check_residual_norm) {
printf("wrong derivative or linear system precision\n");
exit(1);
}
}
node_force_volume(problem,solution_old,dt,NULL,NULL);
free(dx);
free(rhs);
free(solution_old);
}
static void fluid_problem_compute_node_volume(FluidProblem *problem) {
free(problem->node_volume);
Mesh *mesh = problem->mesh;
......@@ -1304,23 +1248,16 @@ static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity,
}
}
FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, const char *petsc_solver_type, int drag_in_stab, double drag_coeff_factor, int temporal, int advection)
FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int temporal, int advection)
{
if (n_fluids != 1 && n_fluids != 2) {
printf("error : only 1 or 2 fluids are supported\n");
}
static int initialized = 0;
if (!initialized) {
hxtInitializeLinearSystems(0, NULL);
initialized = 1;
#ifdef HXT_HAVE_PETSC
hxtPETScInsertOptions(petsc_solver_type, "fluid");
#endif
}
FluidProblem *problem = (FluidProblem*)malloc(sizeof(FluidProblem));
problem->n_fluids = n_fluids;
problem->sigma = sigma;
problem->drag_coeff_factor = drag_coeff_factor;
problem->ip_factor = ip_factor;
problem->stab_param = 0;
problem->drag_in_stab = drag_in_stab;
problem->reduced_gravity = 0;
......@@ -1355,7 +1292,6 @@ FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho
problem->oldporosity = NULL;
problem->solution = NULL;
problem->concentration = NULL;
problem->linear_system = NULL;
problem->n_particles = 0;
problem->particle_uvw = NULL;
problem->particle_element_id = NULL;
......@@ -1389,8 +1325,6 @@ void fluid_problem_free(FluidProblem *problem)
free(problem->bulk_force);
free(problem->oldporosity);
free(problem->node_volume);
if(problem->linear_system)
hxtLinearSystemDelete(&problem->linear_system);
free(problem->particle_uvw);
free(problem->particle_element_id);
free(problem->particle_position);
......@@ -1689,13 +1623,6 @@ static void fluid_problem_set_mesh(FluidProblem *problem, Mesh *mesh) {
for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
problem->solution[i] = 0.;
}
if (problem->linear_system)
hxtLinearSystemDelete(&problem->linear_system);
#ifdef HXT_HAVE_PETSC
hxtLinearSystemCreatePETSc(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements,"fluid");
#else
hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
#endif
fluid_problem_compute_node_volume(problem);
if (problem->boundaries) mesh_boundaries_free(problem->boundaries,problem->mesh->n_boundaries);
problem->boundaries = mesh_boundaries_create(problem->mesh);
......@@ -1799,16 +1726,8 @@ void fluid_problem_after_import(FluidProblem *problem)
mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
fluid_problem_compute_node_volume(problem);
compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->n_particles, problem->particle_position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw);
if(problem->linear_system)
hxtLinearSystemDelete(&problem->linear_system);
#ifdef HXT_HAVE_PETSC
hxtLinearSystemCreatePETSc(&problem->linear_system,problem->mesh->n_elements,N_N,n_fields,problem->mesh->elements,"fluid");
#else
hxtLinearSystemCreateLU(&problem->linear_system,problem->mesh->n_elements,N_N,n_fields,problem->mesh->elements);
#endif
}
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact, long int *elid, int init, int reload)
{
......@@ -2055,3 +1974,8 @@ void weak_boundary_values(const Mesh *mesh, const MeshBoundary *bnd, const WeakB
free(x);
}
}
void fluid_problem_reset_boundary_force(FluidProblem *problem) {
for (int i=0;i<D*problem->mesh->n_boundaries;i++)
problem->boundary_force[i] = 0;
}
fluid/fluid_problem.h
View file @ e90ed53f
......@@ -26,7 +26,6 @@
#include "mesh.h"
#include "mesh_find.h"
#include "hxt_linear_system.h"
#if DIMENSION==2
#define N_N 3
......@@ -60,6 +59,7 @@ struct FluidProblem {
double *g;
double coeffStab;
double drag_coeff_factor;
double ip_factor;
double sigma;
double kinetic_energy;
Mesh *mesh;
......@@ -77,7 +77,6 @@ struct FluidProblem {
StrongBoundary *strong_boundaries;
int n_weak_boundaries;
WeakBoundary *weak_boundaries;
HXTLinearSystem *linear_system;
double *boundary_force;
int n_particles;
double *particle_uvw;
......@@ -105,12 +104,11 @@ struct FluidProblem {
// complete force on the particle (including gravity)
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
void fluid_problem_implicit_euler(FluidProblem *problem, double dt, double