Inhomogeneous Dirichlet BC's implemented into test file advec_diffusion1 and fokkerplanck2_6p1_withLHS.

PDE/advec_diffusion1.m

@@ -0,0 +1,154 @@
+function pde = advec_diffusion1
+% Example PDE using the 1D Advection-Diffusion Equation. This example PDE is
+% time dependent (although not all the terms are time dependent). This
+% implies the need for an initial condition. 
+% PDE:
+% 
+% df/dt = -df/dx + d^2 f/dx^2 + (cos(x) - sin(x))exp(-t)
+%
+% Domain is [0, pi/4]
+% Inhomogeneous Dirichlet boundary condition 
+% Code is expected to be added to the 1d Advection Equation 
+%
+% Diffusion terms are dealt with via LDG, i.e., splitting into two first
+% order equations:
+%
+% d^2 f / dx^2 becomes
+%
+% dq/dx with free (homogeneous Neumann BC)
+%
+% and
+%
+% q=df/dx with Dirichlet BCs specified by analytic solution
+%
+% Run with
+%
+% explicit
+% asgard(advec_diffusion1);
+%
+% implicit
+% asgard(advec_diffusion1,'implicit',true);
+
+pde.CFL = 0.01;
+
+%% Setup the dimensions
+% 
+% Here we setup a 1D problem (x,y)
+
+soln_x = @(x) cos(x) + sin(x);
+soln_t = @(t) exp(-t);
+
+BCFunc = @(x) soln_x(x);
+BCFunc_t = @(t) soln_t(t);
+
+% Domain is (a,b)
+
+% The function is defined for the plane
+% x = a and x = b
+BCL_fList = { ...
+    @(x,p,t) BCFunc(x), ... % replace x by a
+    @(t,p) BCFunc_t(t)
+    };
+
+BCR_fList = { ...
+    @(x,p,t) BCFunc(x), ... % replace x by b
+    @(t,p) BCFunc_t(t)
+    };
+%Incorrect BC to check whether each term requires BC
+%Each term DOES require correct BC
+%BCIncorrect_fList= { ...
+%    @(x,p,t) x.*0 + 2, ... % replace x by b
+%    @(t,p) BCFunc_t(t)
+%    };
+
+dim_x.name = 'x';
+dim_x.domainMin = 0;
+dim_x.domainMax = pi/4;
+dim_x.init_cond_fn = @(x,p,t) soln_x(x)*soln_t(t);
+
+
+%%
+% Add dimensions to the pde object
+% Note that the order of the dimensions must be consistent with this across
+% the remainder of this PDE.
+
+pde.dimensions = {dim_x};
+num_dims = numel(pde.dimensions);
+
+%% Setup the terms of the PDE
+%
+% Here we have 2 terms, with each term having nDims (x) operators.
+
+%%
+% Setup the -d_dx term
+
+g1 = @(x,p,t,dat) x.*0-1;
+pterm1 = GRAD(num_dims,g1,-1,'D','D', BCL_fList, BCR_fList);
+
+term1_x = TERM_1D({pterm1});
+term1   = TERM_ND(num_dims,{term1_x});
+
+%% 
+% Setup the d^2_dx^2 term
+
+% term1
+%
+% eq1 :  df/dt   == d/dx g1(x) q(x,y)   [grad,g1(x)=1, BCL=N, BCR=N]
+% eq2 :   q(x,y) == d/dx g2(x) f(x,y)   [grad,g2(x)=1, BCL=D, BCR=D]
+%
+% coeff_mat = mat1 * mat2
+
+g1 = @(x,p,t,dat) x.*0+1;
+g2 = @(x,p,t,dat) x.*0+1;
+
+pterm1 = GRAD(num_dims,g1,+1,'N','N');
+pterm2 = GRAD(num_dims,g2,-1,'D','D', BCL_fList, BCR_fList);
+
+term2_x = TERM_1D({pterm1,pterm2});
+term2   = TERM_ND(num_dims,{term2_x});
+
+%%
+% Add terms to the pde object
+
+ pde.terms = {term1, term2};
+
+%% Construct some parameters and add to pde object.
+%  These might be used within the various functions below.
+
+params.parameter1 = 0;
+params.parameter2 = 1;
+
+pde.params = params;
+
+%%
+% Sources
+
+s1x = @(x,p,t) cos(x) - sin(x);
+s1t = @(t,p) exp(-t);
+source1 = {s1x, s1t};
+
+pde.sources = {source1};
+
+%% Define the analytic solution (optional).
+% This requires nDims+time function handles.
+
+pde.analytic_solutions_1D = { ...
+    @(x,p,t) soln_x(x), ...
+    @(t,p) soln_t(t) 
+    };
+
+    function dt=set_dt(pde,CFL)
+
+        dims = pde.dimensions;
+
+        % for Diffusion equation: dt = C * dx^2
+
+        lev = dims{1}.lev;
+        dx = 1/2^lev;
+        dt = CFL*dx^2;
+
+    end
+
+pde.set_dt = @set_dt;
+
+end

PDE/fokkerplanck2_6p1.m

@@ -110,7 +110,7 @@
 
 %% Setup the dimensions 
 
-dim_p.domainMin = 0;
+dim_p.domainMin = 0.01;
 dim_p.domainMax = +10;
 dim_p.init_cond_fn = @(x,p,t) f0_p(x);
 

PDE/fokkerplanck2_6p1_withLHS.m

@@ -15,7 +15,7 @@
 % asgard(fokkerplanck2_6p1_withLHS)
 %
 % implicit
-% asgard(fokkerplanck2_6p1_withLHS,'implicit',true,'num_steps',20,'CFL',1.0,'deg',3,'lev',4)
+% asgard(fokkerplanck2_6p1_withLHS, 'implicit', true,'CFL', 1.0, 'num_steps', 50, 'lev', 4, 'deg', 4)
 
 pde.CFL = 0.01;
 
@@ -32,7 +32,7 @@
 tau = 10^5;
 gamma = @(p)sqrt(1+(delta*p).^2);
 vx = @(p)1/vT*(p./gamma(p));
-p_min = 0.005;
+p_min = 0.01;
 
 Ca = @(p)nuEE*vT^2*(psi(vx(p))./vx(p));
 
@@ -43,7 +43,6 @@
     function ret = phi(x)
         ret = erf(x);
     end
-
     function ret = psi(x,t)     
         dphi_dx = 2./sqrt(pi) * exp(-x.^2);
         ret = 1./(2*x.^2) .* (phi(x) - x.*dphi_dx);   
@@ -115,23 +114,23 @@
             case 0.005
                 Q = 0.49859;    % for p_min = 0.005
             case 0.01              
-                Q = 0.497179;   % for p_min = 0.01
+                Q = 0.502844;  % for p_min = 0.01
             case 0.1
                 Q = 0.471814;   % for p_min = 0.1
             case 0.2
                 Q = 0.443769;   % for p_pim = 0.2
             case 0.5
                 Q = 0.361837;   % for p_min = 0.5
             case 1.0
-                Q = 0.23975;     % for p_min = 1.0
+                Q = 0.23975;    % for p_min = 1.0
         end
         ret = Q * 4/sqrt(pi) * exp(-x.^2);
     end
 
 %% Setup the dimensions 
 
 dim_p.domainMin = p_min;
-dim_p.domainMax = +10;
+dim_p.domainMax = 5;
 dim_p.init_cond_fn = @(x,p,t) f0_p(x);
 
 dim_z.domainMin = -1;
@@ -143,12 +142,38 @@
 % Note that the order of the dimensions must be consistent with this across
 % the remainder of this PDE.
 
+
+
 pde.dimensions = {dim_p,dim_z};
 num_dims = numel(pde.dimensions);
 
+%Inhomogeneous BC's set to analytical solution
+BCL_fList = { ...
+    @(x,p,t) soln_p(x,t), ... %replace x with p_min
+    @(x,p,t) soln_z(x,t), ...
+    @(t,p) 1
+    };
+
+BCR_fList = { ...
+    @(x,p,t) soln_p(x,t), ... %replace x with 10
+    @(x,p,t) soln_z(x,t), ...
+    @(t,p) 1
+    };
+
+%BCL_rList = { ...
+%    @(x,p,t) -2*x.*soln_p(x,t), ... %replace x with p_min
+%    @(x,p,t) soln_z(x,t), ...
+%    @(t,p) 1
+%    };
+%BCR_rList = { ...
+%    @(x,p,t) -2*x.*soln_p(x,t), ... %replace x with 10
+%    @(x,p,t) soln_z(x,t), ...
+%    @(t,p) 1
+%    };
+
 %% Setup the terms of the PDE
 
-%%
+%%r
 % termLHS == p^2 d/dt f(p,z,t)
 
 g1 = @(x,p,t,dat) x.^2;
@@ -167,11 +192,13 @@
 % termC1 == d/dp g2(p) r(p)   [grad, g2(p) = p^2*Ca, BCL=D,BCR=N]        
 %   r(p) == d/dp g3(p) f(p)   [grad, g3(p) = 1,      BCL=N,BCR=D]
 
+% Inhomogenous Dirichlet BC's for f(p_max)
+
 g2 = @(x,p,t,dat) x.^2.*Ca(x);
 g3 = @(x,p,t,dat) x.*0+1; 
 
 pterm2  = GRAD(num_dims,g2,+1,'D','N');
-pterm3  = GRAD(num_dims,g3,-1,'N','D');
+pterm3  = GRAD(num_dims,g3,-1,'N','D', BCL_fList, BCR_fList);
 term1_p = TERM_1D({pterm2,pterm3});
 termC1  = TERM_ND(num_dims,{term1_p,[]});
 
@@ -184,7 +211,7 @@
 
 g2 = @(x,p,t,dat) x.^2.*Cf(x);
 
-pterm2  = GRAD(num_dims,g2,+1,'N','D');
+pterm2  = GRAD(num_dims,g2, +1,'N','D', BCL_fList, BCR_fList);
 term2_p = TERM_1D({pterm2});
 termC2   = TERM_ND(num_dims,{term2_p,[]});
 
@@ -202,8 +229,8 @@
 pterm1  = MASS(g1);
 term3_p = TERM_1D({pterm1});
 
-g2 = @(x,p,t,dat) (1-x.^2);
-g3 = @(x,p,t,dat) x.*0+1;
+g2 = @(x,p,t,dat) 1-x.^2;
+g3 = @(x,p,t,dat) x.*0 + 1;
 pterm1  = GRAD(num_dims,g2,+1,'D','D');
 pterm2  = GRAD(num_dims,g3,-1,'N','N');
 term3_z = TERM_1D({pterm1,pterm2});
@@ -213,6 +240,7 @@
 %%
 % Add terms to the pde object
 
+
 pde.terms = {termC1,termC2,termC3};
 
 %% Construct some parameters and add to pde object.

PDE/fokkerplanck2_complete.m

@@ -58,7 +58,7 @@
     case '6p1b'
         delta = 0.042;
         Z = 1;
-        E = 0.25;
+        E = 0.1;
         tau = 10^5;
     case '6p1c' 
         delta = 0.042;
@@ -143,8 +143,8 @@
 
 %% Setup the dimensions 
 
-dim_p.domainMin = 0.1;
-dim_p.domainMax = +10;
+dim_p.domainMin = 0.01;
+dim_p.domainMax = +5;
 dim_p.init_cond_fn = @(x,p,t) f0_p(x);
 
 dim_z.domainMin = -1;

apply_A.m

@@ -31,9 +31,11 @@
 implicit = opts.implicit;
 
 totalDOF = nWork * elementDOF;
+A = 0;
 if opts.implicit
     A = zeros(totalDOF,totalDOF); % Only filled if implicit
 end
+ALHS = 0;
 if num_terms_LHS > 0
     ALHS = zeros(totalDOF,totalDOF); % Only filled if non-identity LHS mass matrix
 end