[CIG-MC] Re: CitcomCU in 2D?

Magali Billen billen at geology.ucdavis.edu
Tue Nov 6 14:29:08 PST 2007


Hello Eh and Tobias,

I've been using CitcomCU in 2D using the conjugate gradient solver.

However, I found that there was a line missing from the function  
solve_del2_u
when choosing the conjugate gradient solver.  In the version I have  
(1.0.2),
the parameter called valid was not set if you choose the conjugate  
gradient solver.
This parameter is returned by the call to this function and is needed  
for the program
to run correctly.

I added the line below marked by /* MIB, FEB 07 */. This should be  
line 416 in
General_matrix_functions.c

      if(!(E->control.NMULTIGRID || E->control.EMULTIGRID))
	{
	  valid = (residual < acc) ? 0 : 1;  /* MIB, FEB 07 */
	  cycles = E->control.v_steps_low;
	  time = CPU_time0();
	  residual = conj_grad(E, D1, r, Au, acc, &cycles, high_lev);
	  for(i = 0; i < neq; i++)
	    {
	      d0[i] += D1[i];
	      D1[i] = 0.0;
	    }

- Magali


On Nov 6, 2007, at 2:04 PM, Eh Tan wrote:

> Hi Tobias,
>
> You can run "2D" model by using only 2 nodes in y direction. A small
> (but non-zero) dimeny might be required to ensure the element has an
> aspect ratio close to 1. You also have to use the conjugate gradient
> solver, not the multigrid solver. The multigrid solver has special
> requirement for the number of nodes in each direction.
>
> A working input file is attached. I only confirmed that the solver is
> working and converging, but didn't check the output, so the model may
> not make sense. Customize the input file for your need, and let me  
> know
> if you have any problem running it.
>
> Cheers,
>
> -- 
> Eh Tan
> Staff Scientist
> Computational Infrastructure for Geodynamics
> 2750 E. Washington Blvd. Suite 210
> Pasadena, CA 91107
> (626) 395-1693
> http://www.geodynamics.org
>
>
>
>
> Tobias Hoeink wrote:
>> Hello from Houston,
>>
>> Is there a way to run CitcomCU (cartesian geometry) in 2  
>> dimensions? I
>> have looked into the source a little and it looks like there is a
>> 'cart2d' option. However, it appears disabled.
>>
>> Would it work to use 'cart3d' and to play with the geometry
>> parameters, e.g. by setting dimeny=0.0 or the like?
>>
>> I would appreciate your help.
>>
>> Thanks,
>> Tobias
>>
>>
>>
>>
>>
>>
>> Tobias Höink
>> Keith-Wiess Postdoctoral Fellow
>>
>> Department of Earth Science, MS-126
>> Rice University
>> 6100 Main St, Houston, Texas 77005
>>
>> Mailing Address:
>> P.O. Box 1892, Houston, Texas 77251-1892
>>
>> Tel. (713) 348 4497
>> Fax. (713) 348 5214
>>
>> Tobias.Hoeink at rice.edu <mailto:Tobias.Hoeink at rice.edu>
>> http://terra.rice.edu/department/staff/tobias.hoeink/
>>
>>
>>
>>
>
> # Start up file for citcom using getpar.
>
>
> # 1. Input and Output Files Information
>
> 	datafile="CASE1/caseA"	# directory and file head for output files.
> 	use_scratch="local"	# output files under current directory.
> #	use_scratch="szhong"	# output files onto local disk to each  
> compute node.
>
> 	oldfile="CASE1/caseA"	# directory and file head for restart files.
> 	restart=0		# check Convection.c for different options
> 	restart_timesteps=20000	# timestep to restart
>
> 	stokes_flow_only=0	# 1: only solve velocity once; 0: time- 
> dependent problem.
> 	maxstep=20000 		# max time steps
> 	storage_spacing=10	# write data every ...
>
> # 2. Geometry, Ra numbers, Internal heating, Thermochemical/Purely  
> thermal convection
>
> 	Solver=cgrad node_assemble=1	 # conjugate gradient: unsure if it  
> still works
>
> 	rayleigh=10.97394e5	# Rayleigh number
> 	rayleigh_comp=1e6	# Compositional Rayleigh number. relevant for  
> composition=1
> 	composition=0		# 0: purely thermal convection; 1: thermochem
> 	Q0=0			# Dimensionless internal heating rate
> 	Q0_enriched=0		# Q0 for C=1 layer, when relevant
>
> 	markers_per_ele=15	# number of particles per element
> 	comp_depth=0.605	# initial depth of compositional interface
>
> 	visc_heating=0		# 1: visc heating on; or 0: visc heating off.
> 	adi_heating=0		# 1: adiabatic heating on; or 0: adiabatic heating  
> off.
> 	
> # 3. Grid And Multiprocessor Information
>
> 	nprocx=1		# number of processors in x or theta dir	
> 	nprocz=1		# number of processors in z or r dir
> 	nprocy=1		# number of processors in y or fi dir
>
> 	nodex=15  nodez=15  nodey=2 	# only relevant with conj-grad
>
> 	mgunitx=14		# multigrid base level in x or t
> 	mgunitz=14		# multigrid base level in z or r
> 	mgunity=1		# multigrid base level in y or f
> 	levels=1		# and how many times it gets doubled
>
> # 4. Coordinate Information
>
> 	Geometry=cart3d 	# cart3d or Rsphere
>
> #CART:  irrelevant if Geometry=Rsphere.
> 	dimenx=1.0		# box size in x direction
> 	dimenz=1.0		# box size in z direction
> 	dimeny=0.01		# box size in y direction
>
> 	z_grid_layers=4		# minus 1 is number of layers with uniform grid  
> in z.
> 	zz=0.0,0.1,0.9,1.0	#    starting and ending z coodinates
> 	nz=1,3,13,15		#    starting and ending node in z direction
>
> 	x_grid_layers=2
> 	xx=0,1
> 	nx=1,15
>
> 	y_grid_layers=2
> 	yy=0,1
> 	ny=1,2
>
> 	z_lmantle=0.76655052		# 2870 km-670 km
> 	z_410=0.857143			# 2870 km-410 km
> 	z_lith=0.9651568		# 2870-100 km
>
>
> # Regional SPHERICAL: irrelevant if Geometry=cart3d
>
> 	radius_inner=0.55	radius_outer=1.0 	# inner and outer radius
> 	theta_north=73.5	theta_south=106.5	# colatitude in deg
> 	fi_west=0		fi_east=36.0 		# longitude in deg
>
> 	r_grid_layers=4
> 	rr=0.55,0.59,0.96,1.0
> 	nr=1,6,44,49
>
> 	t_grid_layers=2
> 	tt=73.5,106.5
> 	nt=1,49
>
> 	f_grid_layers=2
> 	ff=0,36
> 	nf=1,49
>
> 	r_lmantle=0.89482	# (Ro-670 km)/Ro
> 	r_410=0.9356358
> 	r_lith=0.984301		# (Ro-100 km)/Ro
>
> # 5. Rheology
>
> 	rheol=0	# 0,1,2,... diff. option for diff. rheology. Check Visco....c
> 	TDEPV=off		# on/off for temperature-dependent viscosity
> 	VISC_UPDATE=off		# on/off for updating viscosity or not.
> 	update_every_steps=2 	# update viscosity every n timesteps
>
> 	num_mat=4 # number of material group with possibly different rheology
> 	visc0=1.0e0,1.0e0,1.0e0,1.0e0	# pre-exponential constant
> 	viscE=6.9077553,6.9077553,6.9077553,6.9077553	#activation energy
>
> 	viscT=273,273,273,273		#surface temperature:
> 	viscZ=5e-6,5e-6,5e-6,5e-6	#activation volume
>
> 	SDEPV=off			# on/off for non-Newtonian
> 	sdepv_misfit=0.010		# accuracy for non-Newtonian iteration
> 	sdepv_expt=1,1,1,1		# for each mat group, stress exponent
> 	sdepv_trns=1.e0,1.e0,1.e0,1.e0	# transition stress
>
> 	VMIN=on visc_min=5.0e-2		# viscosity lower cutoff
> 	VMAX=on visc_max=2.0e04		# viscosity upper cutoff
>
> 	visc_smooth_cycles=1		# how viscosity is smoothed in multigrid
> 	Viscosity=system		# always
>
> # 6. DIMENSIONAL INFORMATION and Depth-dependence
>
> 	layerd=2870000.0 #meter
> 	radius=6370000.0
> 	ReferenceT=3800.0
> 	refvisc=1.0e20
> 	density=3300.0
> 	thermdiff=1.0e-6
> 	gravacc=9.8
> 	thermexp=5e-5
> 	cp=1250
> 	wdensity=0.0
>
> 	visc_factor=1.0
> 	thermexp_factor=1.0
> 	thermdiff_factor=1.00
> 	dissipation_number=2.601
> 	surf_temp=0.078947
>
> # 7. phase changes: to turn off any of the phase changes, let Ra_XXX=0
>
> 	Ra_410=0.0 #kg/m^3
> 	Ra_670=0.0
> 	clapeyron410=3.0e6 #Pa K-1
> 	clapeyron670=-3.0e6
>
> 	width410=3.5e4 #meter
> 	width670=3.5e4
>
>
> # 8. BOUNDARY CONDITIONS and Initial perturbations
> 	
> 	topvbc=0		# velocity boundary conditions top and bottom
> 		topvbxval=0.0
> 		topvbyval=0.0
> 	botvbc=0
> 		botvbxval=0.0
> 		botvbyval=0.0
> 					#
> 	toptbc=1 bottbc=1		# temperature bc's top and bottom
> 	toptbcval=0.0 bottbcval=1.0	#
>
> 	periodicx=off 		#
> 	periodicy=off		#
> 	flowthroughx=off	#
> 	flowthroughy=off	#
>
> 	num_perturbations=1	#  N, Number of perturbations
> 	perturbmag=0.001	#  A list of N magnitudes
> 	perturbk=1.0		#  A list of N wavenumbers (/PI)
> 	perturbl=6.0		#  A list of N wavenumbers (/PI)
> 	perturbm=0.0		#  A list of N wavenumbers (/PI)
>
>
> # 9. SOLVER RELATED MATTERS
>
> 	Problem=convection	# always, almost
> 	aug_lagr=on
> 	aug_number=1.0e3
> 	precond=on
> 	orthogonal=off
>
> 	maxsub=1
>
> 	viterations=2		# Uzawa iteration loops.
> 	vhighstep=3		# Smoothing passes at highest level (finest grid).
>
> 	piterations=375		# Uzawa iteration loops.
> 	accuracy=1.0e-4		# Desired accuracy of Uzawa algorithm.
> 	tole_compressibility=1e-7
>
> # Tuning of energy equation
>
> 	adv_sub_iterations=2
> 	finetunedt=0.75
>
> 	ll_max=20
> 	nlong=180
> 	nlati=90
>
> # Data input and program debugging
>
> 	DESCRIBE=off		#
> 	BEGINNER=off		#
> 	VERBOSE=off		#
> 	verbose=off		#
> 	COMPRESS=off		#
> 	see_convergence=1
>
> # vim:ts=8:sw=8
> _______________________________________________
> CIG-MC mailing list
> CIG-MC at geodynamics.org
> http://geodynamics.org/cgi-bin/mailman/listinfo/cig-mc

---------------------------------------------
Assistant Professor
U.C. Davis, Geology
(530) 754-5696
billen at geology.ucdavis.edu
----------------------------------------------
"Life is not measured by the breaths you take, but by the moments  
that take your breath away."


-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://geodynamics.org/pipermail/cig-mc/attachments/20071106/0bfb79be/attachment.html


More information about the CIG-MC mailing list