3 % computes a numerical diffusive flux
for a diffusion problem
5 % This
function computes a numerical diffusive flux matrix
for a diffusion
6 % problem. Dirichlet boundary treatment is performed, but neuman values are
set
7 % to
'NaN'. These must be handled by the calling
function.
9 % The analytical term inspiring
this flux looks like
10 % `` \nabla \cdot \left(d(u) \nabla l(u) \right). ``
11 % Here, the functions `d(u)` and `l(u)` differ in the numerical implementation.
12 % The latter is evaluated in cell centers during flux computation, whereas the
13 % diffusion functional `d(u)` is evaluated on the cell edges by averaging the
14 % solutions value there.
16 % The implemented flux then can be written
for inner boundaries as
17 % `` g_{ij}(u,v) = d(\{u,v\}) \frac{l(v) - l(u)}{|d_{ij}|} |e_{ij}| ``
18 % and
for Dirichlet boundaries as
19 % `` g_{dir,ij}(u) = d(u_{dir}(x_{ij}))
20 % \frac{l(u_{dir}(x_{ij})) - l(u)}{|d_{ij}|} |e_{ij}|, ``
21 % where `\{u,v\}` refers to the average of `u` and `v`.
24 % P: Dof vector of finite volume
function `p`, i.e. pressure argument
for
27 % required fields of model:
28 % diffusivity_ptr:
function pointer to diffusivity
function `d`
29 % laplacian_ptr:
function pointer to non-linear function `l`
31 % generated fields of num_flux_mat:
32 % epsilon: diffusion coefficient bound
33 % G: the matrix of numerical flux values across the edges.
34 % boundary treatment is performed as follows:
35 % in dirichlet-boundary edges, the neighbour-values are simply set
36 % to this value and the flux is computed.
39 % Bernard Haasdonk 10.4.2006
41 grid = model_data.grid;
44 num_flux_mat.epsilon = 0;
45 num_flux_mat.G = zeros(model_data.gn_edges, 1);
47 % determine index sets of boundary types
48 real_NB_ind = find(grid.NBI>0 & grid.NBI < repmat((1:grid.nelements)
', 1, grid.nneigh));
49 % NBIind = grid.NBI(real_NB_ind);
50 % INBind = grid.INB(real_NB_ind);
51 % NB_real_NB_ind = sub2ind(size(grid.NBI),NBIind,INBind);
53 dir_NB_ind =find(grid.NBI == -1);
54 neu_NB_ind =find(grid.NBI == -2);
56 if model.verbose >= 29
57 disp(['found
',num2str(length(real_NB_ind)),' non-boundary edges.
'])
58 disp(['found
',num2str(length(dir_NB_ind)),' dirichlet-boundary edges.
'])
59 disp(['found
',num2str(length(neu_NB_ind)),' neumann-boundary edges.
'])
62 tmpS = repmat(S, 1, grid.nneigh);
64 real_nb_ind = grid.NBI > 0;
65 neiS(real_nb_ind) = S(grid.NBI(real_nb_ind));
66 %edgeU = 0.5 * (tmpU + neiU);
67 %edgeU = sqrt(tmpU.*neiU);
68 %edgeU = min(tmpU, neiU);
69 %edgeU = 2*(tmpU.*neiU)./(tmpU+neiU);
71 difftmp = model.two_phase.total_mobility([grid.ESX(:), grid.ESY(:)], tmpS(:), model);
72 diffnei = model.two_phase.total_mobility([grid.ESX(:), grid.ESY(:)], neiS(:), model);
74 meanptr = model.mean_ptr;
76 diff.K = meanptr([difftmp.K, diffnei.K]);
77 diff.epsilon = meanptr(max([difftmp.epsilon, diffnei.epsilon],0));
78 if length(diff.K) == 1
79 K = ones(size(grid.ESX))*diff.K;
81 K = reshape(diff.K,size(grid.ESX));
83 nfaces = size(grid.NBI,2);
85 PP = repmat(P,1,nfaces);
86 % matrix with the neighbouring values products
87 PNB = nan * ones(size(PP));
88 PNB(real_NB_ind) = PP(grid.NBI(real_NB_ind));
90 % set real neighbour values
91 orientation = grid.NX(real_NB_ind) + grid.NY(real_NB_ind);
92 num_flux_mat.G(1:model_data.gn_inner_edges) = ...
93 K(real_NB_ind).*grid.EL(real_NB_ind).*(grid.DS(real_NB_ind).^(-1).* ...
94 (PNB(real_NB_ind)- ...
98 % determine dirichlet-boundary values as required by convective
100 if ~isempty(dir_NB_ind)
101 Xdir = grid.ESX(dir_NB_ind);
102 Ydir = grid.ESY(dir_NB_ind);
103 Sdir = model.dirichlet_values_s_ptr([Xdir,Ydir],model);
104 Pdir = model.dirichlet_values_p_ptr([Xdir,Ydir],model);
105 Kdir = model.two_phase.toal_mobility([Xdir,Ydir], Sdir, model);
107 % determine distances circumcenter to edge-midpoint for
108 % gradient approximation
109 %[dir_NB_i,dir_NB_j] = ind2sub(size(UU),dir_NB_ind);
110 %Ddir = sqrt((grid.ESX(dir_NB_ind)-grid.CX(dir_NB_i)).^2 + ...
111 % (grid.ESY(dir_NB_ind)-grid.CY(dir_NB_i)).^2);
112 % set dirichlet neighbour values
113 %num_flux_mat.G(dir_NB_ind) = ...
114 % grid.EL(dir_NB_ind).*(Ddir.^(-1).* Kdir .*(UU(dir_NB_ind)-Udir));
115 orientation = grid.NX(dir_NB_ind) + grid.NY(dir_NB_ind);
116 num_flux_mat.G(model_data.gn_inner_edges+1:end) = ...
117 grid.EL(dir_NB_ind).*((1.0*grid.DS(dir_NB_ind)).^(-1).* Kdir ...
118 .*(PP(dir_NB_ind)-Pdir)) ...
122 % set diffusivity bound
123 num_flux_mat.epsilon = diff.epsilon;