Commit 867f77d5 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

sub-iterate advection

parent b1b2c161
Pipeline #4768 passed with stage
in 23 seconds
......@@ -1548,6 +1548,11 @@ void fluid_problem_boundary_name(const FluidProblem *p, int i, char **phys_name)
*phys_name = p->mesh->boundary_names[i];
}
double *fluid_problem_element_size(const FluidProblem *p)
{
return p->element_size;
}
/*void fluid_problem_private_set_elements(FluidProblem *p, int n_nodes, double *x, int n_elements, int *elements)
{
Mesh *m = malloc(sizeof(Mesh));
......
......@@ -113,6 +113,7 @@ int *fluid_problem_elements(const FluidProblem *p);
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_set_wall_boundary(FluidProblem *p, const char *tag, BoundaryCallback *pressurecb);
int fluid_problem_private_n_physical(const FluidProblem *p);
......
......@@ -274,7 +274,15 @@ class FluidProblem :
"""
nit = self._lib.fluid_problem_implicit_euler(self._fp, c_double(dt), c_double(newton_atol), c_double(newton_rtol), c_int(newton_max_it), c_int(reduced_gravity), c_double(stab_param))
if self._n_fluids == 2 :
self._lib.fluid_problem_advance_concentration(self._fp,c_double(dt))
nv = np.linalg.norm(self.solution()[:,:self._dim],axis=1)
nvmax = np.max(nv[self.elements()],axis=1,keepdims=True)
h = self._get_matrix("element_size",self.n_elements(),1)
cfl = dt*nvmax / (0.2*h)
nsub = int(np.max(cfl)+1)
if (nsub != 1) :
print("sub-iterating advection for cfl : %i sub-iterations"%nsub)
for i in range(nsub) :
self._lib.fluid_problem_advance_concentration(self._fp,c_double(dt/nsub))
return nit
def set_particles(self, mass, volume, position, velocity, reload = False):
......
......@@ -60,8 +60,8 @@ N = 10000
rhop = 1500
# numerical parameters
tEnd = 10 # final time
dt = .001 # time step
outf = 10 # number of iterations between output files
dt = .01 # time step
outf = 1 # number of iterations between output files
def outerBndV(x) :
a = np.where(np.abs(x[:,0]-0.111) < 0.01, 0.4*max(np.sin(t*np.pi*2./1),0),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