test_2d_resistor

PURPOSE ^

Create 2D model of a cylindrical resistor

SYNOPSIS ^

function test_2d_resistor(opt)

DESCRIPTION ^

 Create 2D model of a cylindrical resistor
 $Id: test_2d_resistor.m 5399 2017-04-12 17:16:01Z aadler $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function test_2d_resistor(opt)
0002 % Create 2D model of a cylindrical resistor
0003 % $Id: test_2d_resistor.m 5399 2017-04-12 17:16:01Z aadler $
0004 
0005 if nargin>0 && strcmp(opt,'UNIT_TEST'); do_unit_test; return; end
0006 
0007 resistor_test;
0008 
0009 function img = this_mdl(conduc, z_contact, current);
0010    nn= 12;     % number of nodes
0011    ww=3;       % width = 4
0012    scale = .35;
0013    mdl=mk_grid_model([],3+scale*(1:ww), scale*(1:nn/ww));
0014    mdl= rmfield(mdl,'coarse2fine'); % don't calc this.
0015 
0016    mdl.gnd_node = 1;
0017    elec_nodes= [1:ww];
0018    elec(1).nodes= elec_nodes;      elec(1).z_contact= z_contact;
0019    elec(2).nodes= nn-elec_nodes+1; elec(2).z_contact= z_contact;
0020    stim.stim_pattern= [-1;1]*current;
0021    stim.meas_pattern= [-1,1];
0022    mdl.stimulation= stim;
0023    mdl.electrode= elec;
0024    n_el = size(mdl.elems,1);
0025    img= eidors_obj('image','2D rectangle', ...
0026          'elem_data', ones(n_el,1) * conduc );
0027    img.fwd_model = mdl;
0028    img.fwd_model.normalize_measurements = 0;
0029 
0030 function V = analytic_soln(img);
0031    % analytical solution
0032    nodes = img.fwd_model.nodes;
0033    wid_len= max(nodes) - min(nodes);
0034    conduc =  mean(img.elem_data);
0035    Block_R = wid_len(2) / wid_len(1) / conduc;
0036    % Contact R reflects z_contact / width. There is no need to scale
0037    %  by the scale, since this is already reflected in the size of the
0038    %  FEM as created by the grid. This is different to the test_3d_resistor,
0039    %  where the FEM is created first, and then scaled, so that the ww
0040    %  and hh need to be scaled by the scale parameter.
0041    z_contact = sum([img.fwd_model.electrode(:).z_contact]);
0042    current = max(img.fwd_model.stimulation(1).stim_pattern(:));
0043    Contact_R = z_contact/wid_len(1);
0044    R = Block_R + Contact_R;
0045 
0046    V= current*R;
0047    fprintf('Solver %s: %f\n', 'analytic', V);
0048    fprintf('Solver %s: %f\n', 'analytic (no z_contact)', V - 2*Contact_R*current);
0049 
0050 function V = first_order_solver( img );
0051    % AA_SOLVER
0052    img.fwd_model.solve = @fwd_solve_1st_order;
0053    img.fwd_model.system_mat = @system_mat_1st_order;
0054    fsol= fwd_solve(img);
0055    fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0056    V = fsol.meas;
0057 
0058 function vals = RC_test;
0059 
0060    current= 4;  % Amps
0061    conduc=  .4 + 2*pi*j*10; % conductivity in Ohm-meters
0062    z_contact= 1e-1;
0063    img = this_mdl( conduc, z_contact, current);
0064 
0065    vals.analytic = analytic_soln(img);
0066    vals.aa_solver = first_order_solver(img);
0067 
0068 function vals = resistor_test;
0069 
0070    current= 4;  % Amps
0071    conduc=  .4; % conductivity in Ohm-meters
0072    z_contact= 1e-1;
0073    img = this_mdl( conduc, z_contact, current);
0074 
0075    show_fem(img.fwd_model);
0076 
0077    vals.analytic = analytic_soln(img);
0078 
0079 % AA_SOLVER
0080    vals.aa_solver = first_order_solver(img);
0081 
0082 warn_state = warning('query','EIDORS:Deprecated');
0083 warning('off','EIDORS:Deprecated');
0084 
0085 
0086 % NP_SOLVER
0087 img.fwd_model.solve = @np_fwd_solve;
0088 img.fwd_model.system_mat = @np_calc_system_mat;
0089 fsol= fwd_solve(img);
0090 fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0091 vals.np_solver = fsol.meas;
0092 
0093 % NOW CALCULATE THE ANALYTICAL JACOBIAN
0094 img.fwd_model.solve = @np_fwd_solve;
0095 img.fwd_model.jacobian = @np_calc_jacobian;
0096 
0097 img.elem_data= ones(num_elems(img),1) * conduc ;
0098 Jnp= calc_jacobian(img);
0099 
0100 img.fwd_model.jacobian = @jacobian_perturb;
0101 Jp1= calc_jacobian(img);
0102 
0103 img.elem_data= ones(num_elems(img),1) * conduc ;
0104 Jp2= zeros(size(Jnp));
0105 delta= 1e-8;
0106 for i=1:num_elems(img)
0107    fsol_h= fwd_solve(img);
0108    img.elem_data(i) = img.elem_data(i) + delta;
0109    fsol_i= fwd_solve(img);
0110    Jp2(:,i) = (fsol_i.meas - fsol_h.meas)/delta; 
0111 end
0112    
0113 mdl.solve = @fwd_solve_1st_order;
0114 mdl.jacobian = @jacobian_adjoint;
0115 Jaa= calc_jacobian(img);
0116 
0117 
0118 fprintf('Jacobians: Cols by Jaa, Jnp, Jp1, Jp2:\n')
0119 JJ = [Jaa;Jnp;Jp1;Jp2];
0120 disp(JJ(:,1:6))
0121 
0122 
0123 
0124 
0125 % TEST RESISTOR THAT IS NOT RECTANGULAR
0126 nn= 12;     % number of nodes
0127 ww=2;       % width = 4
0128 if ~exist('conduc');conduc=  .4;end  % conductivity in Ohm-meters
0129 current= 4;  % Amps
0130 z_contact= 1e-2;
0131 scale = .35;
0132 mdl=mk_grid_model([],3+scale*(1:ww), scale*(1:nn/ww));
0133 xdelta= .05*mdl.nodes(1:ww:nn,2);
0134 mdl.nodes(1:ww:nn,1) = mdl.nodes(1:ww:nn,1) + xdelta;
0135 mdl.nodes(2:ww:nn,1) = mdl.nodes(2:ww:nn,1) - xdelta;
0136 mdl.gnd_node = 1;
0137 elec_nodes= [1:ww];
0138 elec(1).nodes= elec_nodes;      elec(1).z_contact= z_contact;
0139 elec(2).nodes= nn-elec_nodes+1; elec(2).z_contact= z_contact;
0140 stim.stim_pattern= [-1;1]*current;
0141 stim.meas_pattern= [-1,1];
0142 mdl.stimulation= stim;
0143 mdl.electrode= elec;
0144 show_fem(mdl);
0145 n_el = size(mdl.elems,1);
0146 img= eidors_obj('image','2D rectangle', 'fwd_model',mdl, ...
0147       'elem_data', ones(n_el,1) * conduc );
0148 % AA_SOLVER
0149 mdl.solve = @fwd_solve_1st_order;
0150 mdl.system_mat = @system_mat_1st_order;
0151 img.fwd_model = mdl;
0152 fsol= fwd_solve(img);
0153 fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0154 % NP_SOLVER
0155 mdl.solve = @np_fwd_solve;
0156 mdl.system_mat = @np_calc_system_mat;
0157 fsol= fwd_solve(img);
0158 fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0159 
0160 % Analytic
0161 e1= mdl.nodes(mdl.electrode(1).nodes,:);
0162 d1= sqrt( sum(( [1,-1]*e1 ).^2 ));
0163 e2= mdl.nodes(mdl.electrode(2).nodes,:);
0164 d2= sqrt( sum(( [1,-1]*e2 ).^2 ));
0165 hig= max(mdl.nodes(:,2)) - min(mdl.nodes(:,2));
0166 % R = int_h (1/d) dh
0167 % d= a+b*h; a= d1; b= (d2-d1)/hig
0168 % R = int 1/(a+b*h)*dh = log(a+bh)/b
0169 % R = [log(d2) - log(d1)]*hig/(d2-d1)
0170 % doesn't account for arc in conductivity pattern
0171 R = (log(d2)-log(d1))*hig/(d2-d1)/conduc + 2*z_contact/scale;
0172 V= current*R;
0173 fprintf('Solver %s: %f\n', 'analytic', V);
0174 fprintf('Analytic is not expected to be same in last case\n');
0175 
0176 warning(warn_state.state,'EIDORS:Deprecated');
0177 
0178 
0179 function do_unit_test
0180   vals = resistor_test;
0181 
0182   unit_test_cmp('Analytic vs AA', vals.analytic, vals.aa_solver, 1e-10);
0183   unit_test_cmp('Analytic vs NP', vals.analytic, vals.np_solver, 1e-10);
0184 
0185   vals = RC_test;
0186 
0187   unit_test_cmp('Analytic vs AA', vals.analytic, vals.aa_solver, 1e-10);
0188
Generated on Fri 30-Dec-2022 19:44:54 by m2html © 2005