# jfhbrook/anisotropy

as far as I can tell sims are now running. Just wanna parallelize now

jfhbrook committed Jul 12, 2010
Showing with 150 additions and 43 deletions.
1. +11 −0 fea/PLAN
2. +31 −0 fea/controller.m
3. +74 −0 fea/mesher.m
4. +13 −43 fea/{solveit.m → solver.m}
5. +21 −0 fea/symmetric_fromcell.m
 @@ -0,0 +1,11 @@ 1) Run mesher, generate mesh mats. SAVE ORIGINAL ANGLE AND PARAMS! 2) Distribute mesh mats 3) Run each mesh mat over the entire range of k's. SAVE Ks! --> Also parameterize metadata 5) Merge data into collections: a) Separate data structures for different parameter collections! b) Group in 3 dimensions ([angle, kxy,kz]) * (We already should have [kxy,kz]--just stack them horizontally? use cat(dim,v1,v2) and permute(a,[i,j,k]) c) struct(metadata, metadata, data, A(angle,kzy,kx)=cell{[time], [Temp]}
 @@ -0,0 +1,31 @@ %controller.m %does the controlling %angles = linspace(0,90,10); angles = [0]; ks = linspace(0.2,0.4,8); [kxy,kz] = meshgrid(ks,ks); flreport('off'); % Some parameters we won't want to iterate through params=struct('rsnow', 0.4, ... 'rneedle', 0.00025, ... 'lneedle', 0.1, ... 'density_snow', 200, ... 'density_needle', 8000, ... 'cp_snow', 2050, ... 'cp_needle', 460, ... 'q_needle', 635500, ... 'k_needle', 160, ... 'time', [colon(0,0.1,1) colon(2,2,10) colon(20,10,100)]); saveroot='./'; for angle=angles, mesh = mesher(angle,params); solutions = arrayfun(@(x,y) solver(x,y,mesh,params), kxy,kz, 'UniformOutput', false); save([saveroot 'solution-' date '-' num2str(angle)],solutions,angle,ks,params); end
 @@ -0,0 +1,74 @@ % COMSOL Multiphysics Model M-file % Generated in part by COMSOL 3.5a (COMSOL 3.5.0.608, \$Date: 2009/05/11 07:38:49 \$) % the rest of it modified by Joshua Holbrook function fem=mesher(angle,params) % mesh_generate(angle) % generates a mesh for the given angle. fprintf ['meshing for angle=' num2str(angle) '...\n']; flclear fem % COMSOL version clear vrsn vrsn.name = 'COMSOL 3.5'; vrsn.ext = 'a'; vrsn.major = 0; vrsn.build = 608; vrsn.rcs = '\$Name: v35ap \$'; vrsn.date = '\$Date: 2009/05/11 07:38:49 \$'; fem.version = vrsn; % Geometry g1=sphere3(num2str(params.rsnow),'pos',{'0','0','0'},'axis',{'0','0','1'},'rot','0'); g2=cylinder3(num2str(params.rneedle),num2str(params.lneedle),'pos',{num2str(-params.lneedle/2),'0','0'},'axis',{'1','0','0'},'rot','15'); parr={point3(0,0,0)}; g3=geomcoerce('point',parr); parr={point3(params.rsnow,0,0)}; g4=geomcoerce('point',parr); parr={point3(0,params.rsnow,0)}; g5=geomcoerce('point',parr); parr={point3(0,0,params.rsnow)}; g6=geomcoerce('point',parr); parr={point3(-params.rsnow,0,0)}; g7=geomcoerce('point',parr); parr={point3(0,-params.rsnow,0)}; g8=geomcoerce('point',parr); parr={point3(0,0,-params.rsnow)}; g9=geomcoerce('point',parr); % Analyzed Geometry (?) clear p s p.objs={g3,g4,g5,g6,g7,g8,g9}; p.name={'ORIGIN','PT1','PT2','PT3','PT4','PT5','PT6'}; p.tags={'g3','g4','g5','g6','g7','g8','g9'}; s.objs={g1,g2}; s.name={'SNOW','NEEDLE'}; s.tags={'g1','g2'}; fem.draw=struct('p',p,'s',s); fem.geom=geomcsg(fem); % ODE Settings clear ode clear units; units.basesystem = 'SI'; ode.units = units; fem.ode=ode; % Initialize mesh fem.mesh=meshinit(fem, ... 'hauto',5, ... 'hgradsub',[2,1.1], ... 'hmaxsub',[2,0.0005]); % Refine mesh fem.mesh=meshrefine(fem, ... 'mcase',0, ... 'rmethod','longest'); fem=multiphysics(fem); end
 @@ -1,45 +1,11 @@ % COMSOL Multiphysics Model M-file % Generated by COMSOL 3.5a (COMSOL 3.5.0.608, \$Date: 2009/05/11 07:38:49 \$) function solution_generate(meshfile,kmatrix) %solution_generate(meshmat,kmatrix) function answer=solver(kxy,kz,fem,params) %solver(kxy,kz,mesh,params) %uses comsol to pump out a solution using a given mesh-mat and a k-matrix in comsol format. flreport('off'); %params saveroot='./'; %It might be a good idea to take this approach. %flbinaryfile='solve.mphm'; load meshfile; %{ % Geometry clear draw g9=flbinary('g9','draw',flbinaryfile); g7=flbinary('g7','draw',flbinaryfile); g8=flbinary('g8','draw',flbinaryfile); g6=flbinary('g6','draw',flbinaryfile); g5=flbinary('g5','draw',flbinaryfile); g4=flbinary('g4','draw',flbinaryfile); g3=flbinary('g3','draw',flbinaryfile); draw.p.objs = {g9,g7,g8,g6,g5,g4,g3}; draw.p.name = {'B_6','B_4','B_5','B_3','B_2','B_1','ORIGIN'}; draw.p.tags = {'g9','g7','g8','g6','g5','g4','g3'}; g2=flbinary('g2','draw',flbinaryfile); g1=flbinary('g1','draw',flbinaryfile); draw.s.objs = {g2,g1}; draw.s.name = {'NEEDLE','SNOW'}; draw.s.tags = {'g2','g1'}; fem.draw = draw; fem.geom = geomcsg(fem); fem.mesh = flbinary('m1','mesh',flbinaryfile); %} fprintf ['solving for kxy=' num2str(kxy) ' and kz=' num2str(kz) '...\n']; % Application mode 1 clear appl appl.mode.class = 'GeneralHeat'; @@ -54,12 +20,16 @@ function solution_generate(meshfile,kmatrix) appl.bnd = bnd; clear equ equ.sdtype = 'gls'; equ.rho = {200,8000}; % densities equ.rho = {params.density_snow,params.density_needle}; equ.init = 0; equ.shape = 2; equ.C = {2050,460}; equ.Q = {0,635500}; equ.k = {{0.2;0.1;0.1},160}; % Heat capacities equ.C = {params.cp_snow,params.cp_needle}; % Wattage equ.Q = {0,params.q_needle}; % Heat conductivities equ.k = {symmetric_tocell(diag([kxy,kxy,kz])),params.k_needle}; equ.ind = [1,2]; appl.equ = equ; fem.appl{1} = appl; @@ -114,7 +84,7 @@ function solution_generate(meshfile,kmatrix) 'solcomp',{'T'}, ... 'outcomp',{'T'}, ... 'blocksize','auto', ... 'tlist',[colon(0,0.1,1) colon(2,2,10) colon(20,10,100)], ... 'tlist', params.time, ... 'estrat',1, ... 'tout','tlist', ... 'linsolver','gmres', ... @@ -162,5 +132,5 @@ function solution_generate(meshfile,kmatrix) 'solnum','end'); save([saveroot meshfile(1:size(a,2)-4) '-matrix-' cell2str(kmatrix) '-' date '.mat'], 'fem'); answer={T_thermistor,T_surf_avg}; end
 @@ -0,0 +1,21 @@ function A=symmetric_fromcell(a) % inverts cell form of symmetric matrix. % dependencies: triangle.m, triangle_inv.m n=triangle_inv(size(a,1)); % should check for cell dimensions. This needs to be tested. assert(mod(n,1)==0, 'cells must be nx1!'); %makes sure n for nxn matrix implied by t is whole assert(size(a,2)==1, 'cells must be nx1 where n is a triangular number!'); %makes sure other dimension is 1 % The heavy lifting. A=[]; for m = 1:n, i=(m^2-m+2)/2:(m^2+m)/2; %slice of cells to wrap around A A=[A,init([a{i}])'; a{i}]; end function n=triangle_inv(t) n = (sqrt(8*t+1)-1)/2; end