Commit 848f2916 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

wip

parent a848e94f
Pipeline #8560 failed with stages
in 1 minute and 32 seconds
......@@ -173,7 +173,7 @@ void fluid_problem_node_force_volume(FluidProblem *problem, const double *soluti
free(sold);
}
static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double *ds, const double *sold, const double *dsold, const double c, const double *dc, const double rho, double mu, const double dt, int eid, const double *data, double *f00, double *f01)
static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double *ds, const double *sold, const double *dsold, const double c, const double *dc, const double rho, double mu, int eid, const double *data, double *f00, double *f01)
{
const int vid = wbnd->vid;
const int pid = wbnd->pid;
......@@ -450,6 +450,118 @@ static void fluid_problem_f(const FluidProblem *problem, const double *sol, doub
}
}
void fluid_problem_weak_boundary_drag(FluidProblem *problem, const char *tag, double *force)
{
const double *solution_old = problem->solution;
const Mesh *mesh = problem->mesh;
const double *solution = problem->solution;
const int n_fields = fluid_problem_n_fields(problem);
const size_t local_size = N_SF*n_fields;
double *s = malloc(sizeof(double)*n_fields);
double *sold = malloc(sizeof(double)*n_fields);
double *ds = malloc(sizeof(double)*n_fields*D);
double *dsold = malloc(sizeof(double)*n_fields*D);
double *f0 = malloc(sizeof(double)*n_fields);
double *f00 = malloc(sizeof(double)*n_fields*n_fields);
double *f01 = malloc(sizeof(double)*n_fields*n_fields*D);
for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
if(strcmp(wbnd->tag,tag) != 0)
continue;
int bndid= -1;
for (int i = 0; i < mesh->n_boundaries; ++i) {
if (strcmp(mesh->boundary_names[i],wbnd->tag) == 0){
bndid = i;
}
}
if (bndid == -1) {
printf("Mesh has no boundary with name \"%s\".", wbnd->tag);
}
MeshBoundary *bnd = &problem->boundaries[bndid];
int n_value = weak_boundary_n_values(wbnd);
double *data = malloc(n_value*bnd->n_interfaces*N_LQP*sizeof(double));
weak_boundary_values(mesh, bnd, wbnd, data);
for (int dim = 0; dim < D; ++dim) {
force[dim] = 0;
}
for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
const int iel = interface[0];
const int icl = interface[1];
const int *cl = elbnd[icl];
int nodes[N_LSF];
double *x[N_LSF];
const int *el = &mesh->elements[iel*N_N];
for (int i = 0; i < D; ++i){
nodes[i] = el[cl[i]];
x[i] = &mesh->x[nodes[i]*3];
}
double dxidx[D][D], dphi[N_N][D];
mesh_dxidx(mesh,iel,dxidx);
double n[D],detbnd;
get_normal_and_det(x,n,&detbnd);
grad_shape_functions(dxidx, dphi);
double dc[D] = {0};
for (int j = 0; j < n_fields*D; ++j) {
ds[j] = 0;
}
for (int i = 0; i < N_SF; ++i) {
for (int j = 0; j < n_fields; ++j) {
double dof = solution[el[i]*n_fields+j];
double dofold = solution_old[el[i]*n_fields+j];
for (int k = 0; k < D; ++k){
ds[j*D+k] += dphi[i][k]*dof;
dsold[j*D+k] += dphi[i][k]*dofold;
}
}
for (int k=0; k<D; ++k){
dc[k] += problem->porosity[el[i]]*dphi[i][k];
}
}
for (int iqp = 0; iqp < N_LQP; ++iqp) {
double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
double dp[D]={0};
double phi[DIMENSION];
l_shape_functions(LQP[iqp],phi);
for (int i = 0; i < n_fields; ++i) {
s[i] = 0;
sold[i] = 0;
f0[i] = 0;
}
double c = 0;
double a = 0;
for (int i = 0; i < DIMENSION; ++i) {
c += problem->porosity[nodes[i]]*phi[i];
if (problem->n_fluids == 2) {
a += problem->concentration[iel*N_N+cl[i]]*phi[i];
}
for (int j = 0; j < n_fields; ++j) {
double dof = solution[nodes[i]*n_fields+j];
double dofold = solution_old[nodes[i]*n_fields+j];
s[j] += phi[i]*dof;
sold[j] += phi[i]*dofold;
}
}
double rho, mu;
fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
const double jw = LQW[iqp]*detbnd;
f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,c,dc,rho,mu,iel,dataqp,f00,f01);
for (int dim = 0; dim < D; ++dim) {
force[dim] += f0[dim]*jw;
}
}
}
free(data);
}
free(s);
free(ds);
free(sold);
free(dsold);
free(f0);
free(f00);
}
static void compute_weak_boundary_conditions(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
const Mesh *mesh = problem->mesh;
......@@ -464,8 +576,6 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
double *f00 = malloc(sizeof(double)*n_fields*n_fields);
double *f01 = malloc(sizeof(double)*n_fields*n_fields*D);
double i_bnd = 0;
double iv_bnd_up = 0;
double iv_bnd_down = 0;
for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
int bndid= -1;
......@@ -555,15 +665,13 @@ static void compute_weak_boundary_conditions(FluidProblem *problem, const double
for (int i = 0; i < D; ++i){
if (wbnd->vid<0) {
i_bnd -= a*s[U+i]*n[i]*jw*dt;
iv_bnd_up -= s[U+i]*n[i]*jw;
}
else {
iv_bnd_down -= dataqp[wbnd->vid+i]*n[i]*jw;
i_bnd -= a*dataqp[wbnd->vid+i]*n[i]*jw*dt;
}
}
}
f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,c,dc,rho,mu,dt,iel,dataqp,f00,f01);
f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,c,dc,rho,mu,iel,dataqp,f00,f01);
for (int ifield = 0; ifield < n_fields; ++ifield) {
for (int iphi = 0; iphi < D; ++iphi){
if (ifield<D){
......@@ -945,7 +1053,7 @@ static void fluid_problem_surface_tension_bnd(FluidProblem *problem, double *a_c
}
}
static void fluid_problem_surface_tension(FluidProblem *problem, const double *solution_old, double *grad_a_cg, double *all_local_vector){
static void fluid_problem_surface_tension(FluidProblem *problem, const double *solution_old, double *grad_a_cg, double *all_local_vector) {
if (problem->n_fluids == 1) return;
const Mesh *mesh = problem->mesh;
double sigma = problem->sigma;
......@@ -1064,7 +1172,7 @@ static void fluid_problem_dg_to_cg_grad(FluidProblem *problem, const double *dg,
}
}
for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd) {
WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
int bndid= -1;
for (int i = 0; i < mesh->n_boundaries; ++i) {
......@@ -1231,9 +1339,7 @@ void fluid_problem_immersed_boundaries(const FluidProblem *problem, int n_partic
size_t local_size = N_SF*n_fields;
CSR *csr = mesh_tree_particle_to_mesh_csr_element(problem->mesh_tree,
n_particles, particle_position, particle_r);
printf("particle_r : %g\n",particle_r[0]);
for (int p = 0; p < n_particles; ++p) {
printf("%i triangles \n",csr->row_index[p+1]-csr->row_index[p]);
for (int ie = csr->row_index[p]; ie < csr->row_index[p+1]; ++ie) {
int iel = csr->col_index[ie];
const int *el = problem->mesh->elements + iel*N_N;
......@@ -1254,6 +1360,37 @@ void fluid_problem_immersed_boundaries(const FluidProblem *problem, int n_partic
csr_free(csr);
}
void fluid_problem_immersed_forces(const FluidProblem *problem, int n_particles, const double *particle_position, const double *particle_r, double tau, const double *solution, const double *solution_old, double *forces)
{
int n_fields = fluid_problem_n_fields(problem);
size_t local_size = N_SF*n_fields;
CSR *csr = mesh_tree_particle_to_mesh_csr_element(problem->mesh_tree,
n_particles, particle_position, particle_r);
for (int p = 0; p < n_particles; ++p) {
for (int id = 0; id < D; ++id) {
forces[p*DIMENSION+id] = 0;
}
for (int ie = csr->row_index[p]; ie < csr->row_index[p+1]; ++ie) {
int iel = csr->col_index[ie];
const int *el = problem->mesh->elements + iel*N_N;
const double a = csr->v[ie*(N_N+1)];
const double *w = csr->v+ie*(N_N+1)+1;
double dxidx[D][D];
const double detj = mesh_dxidx(problem->mesh, iel, dxidx);
double dphi[N_N][D];
grad_shape_functions(dxidx, dphi);
for (int in = 0; in < N_N; ++in) {
for (int id = 0; id < D; ++id) {
double c = tau*a*w[in];///problem->node_volume[el[in]];
forces[p*DIMENSION+id] += c*solution_old[el[in]*n_fields+U+id]
+a*dphi[in][id]*solution_old[el[in]*n_fields+P];
}
}
}
}
csr_free(csr);
}
void fluid_problem_porosity2(const FluidProblem *problem, double *porosity)
{
const Mesh *mesh = problem->mesh;
......
......@@ -337,7 +337,6 @@ class FluidProblem :
reduced_graviy -- if True simulations are run with a reduced gravity (not to use in two fluids simulations)
stab_param -- if zero pspg/supg/lsic stabilisation is computed otherwise the value is used as a coefficient for a pressure laplacian stabilisation term
"""
print(particles_position.shape, particles_r.shape)
self._lib.fluid_problem_set_reduced_gravity(self._fp,c_int(reduced_gravity))
self._lib.fluid_problem_set_stab_param(self._fp,c_double(stab_param))
self._lib.fluid_problem_apply_boundary_conditions(self._fp)
......@@ -371,6 +370,9 @@ class FluidProblem :
if norm > check_residual_norm :
raise ValueError("wrong derivative or linear system precision")
self._lib.fluid_problem_node_force_volume(self._fp,_np2c(solution_old),c_double(dt),None,None)
self._immersed_forces = np.zeros_like(particles_position)
if particles_position.shape[0] > 0 :
self._lib.fluid_problem_immersed_forces(self._fp, c_int(particles_position.shape[0]),_np2c(particles_position),_np2c(particles_r), c_double(tau), _np2c(self.solution()),_np2c(solution_old), _np2c(self._immersed_forces))
def compute_cfl(self):
"""Compute the CFL number divided by the time step
......@@ -515,6 +517,11 @@ class FluidProblem :
"""Give access to the volume force at fluid nodes."""
return self._get_matrix("bulk_force", self.n_nodes(), self._dim)
def get_weak_boundary_forces(self, tag) :
f = np.zeros([self.dimension()])
self._lib.fluid_problem_weak_boundary_drag(self._fp, tag.encode(), _np2c(f))
return f
def node_volume(self) :
"""Return the volume associated with each node"""
return self._get_matrix("node_volume",self.n_nodes(),1)
......
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